Study on Deposition Characteristics of Microparticles in Terminal Pulmonary Acini by IB–LBM

As an indicator of health risk, the deposition of microparticles in terminal pulmonary acini is of great significance in the medical field. To control particulate pollution and optimize aerosol delivery, it is necessary to perform an in-depth study of the microparticle deposition in terminal pulmonary acini; however, little research has been done on this topic. This paper proposes a respiratory movement model of terminal pulmonary acini using an immersed boundary–lattice Boltzmann method. In addition, we explored the effect of gravity direction, respiratory rate, microparticle diameter, and other parameters on the microparticles deposition process and distribution, under the airflow in the acinar wall. It was found that the deposition of microparticles is sensitive to gravity direction, and the growth of the respiratory rate increases the rate of microparticle migration and deposition. It was observed that the gravity effect is enhanced by increasing the diameter of microparticles, causing a high deposition and dispersion rate. The study reveals the dynamic correlation between the respiration process and the movement of microparticles, which is of reference value to figure out the pathogenicity mechanism of inhalable particles and to optimize the aerosol delivery.


Introduction
There are 274-790 million (480 million on average) alveoli in human lungs that constitute a large surface area of 70-75 m 2 for gas exchange between the inhaled air and the blood [1]. An acinus is defined as the blind-ending parenchymal unit beginning with a transitional bronchiole. Within an acinus, all airways (alveolar ducts and respiratory bronchioles) have alveoli attached to their walls, which are the primary sites of oxygen and carbon dioxide exchange [2]. The acini are composed of three distinct structures: respiratory bronchioles, alveolar ducts, and alveolar sacs [3]. Respiratory bronchioles split into alveolar ducts, which, further, branch into alveolar sacs, terminally [4]. Alveolar sacs are of great importance to pulmonary function [5].
Some respiratory diseases affect pulmonary acini, where microparticles are inhaled and deposited, causing infection, inflammation, and even tumor formation [6]. Moreover, certain viruses and bacteria are responsible for lung diseases by invading pulmonary acini directly, such as inhalational anthrax, which results from the deposition of spore-bearing microparticles of 1 to 5 µm into the pulmonary acini [7]. Knowing the deposition of microparticles in terminal pulmonary acini helps predict and prevent respiratory diseases.
Injection is the most common method of drug delivery. Substantial progress has been made in recent years regarding nanoparticle and microparticle delivery in blood circulation. It is found that ellipsoidal particles with submicron shapes are more liable to adhere to the wall surface [8,9]. Drug particles need to overcome a series of complex biological obstacles to reach biological targets in blood circulation. By comparison, pulmonary aerosol delivery is more efficient. As an emerging method, pulmonary aerosol delivery is significant to treat respiratory diseases. It is an effective and non-invasive pathway for therapeutic delivery, especially for drugs that cannot be delivered orally due to degradation in the digestive system [10]. Alveolar sacs are the functional targets for inhalation drug deposition and absorption, where, for example, inhaled insulin and some other drugs are absorbed. The medicines used for distal lungs, such as pentamidine depositing in upper airways, may cause significant side effects [11]. Thus, the size and the way of inhalation of atomized droplets should be well designed.
Based on the reasons above, it is of great scientific and clinical significance to study the deposition of microparticles in pulmonary acini. In some studies, a hemispherical shell structure and a cylindrical pipe were used to simulate the alveolus and bronchiole, respectively, and a preset airflow was designed to drive microparticles. The application of this rigid wall structure helps to study the movement and deposition of microparticles in pulmonary acini [12][13][14]. Compared to the rigid wall, the simulation results of using alveolus with rhythmic wall motion model is more similar to the physiological situation, and the deposition characteristics of microparticles are more accurate, accordingly. Tsuda et al. [15] simulated a two-dimensional axisymmetric alveolar structure, expanding and contracting sinusoidally, and the ratio of alveolar flow to ductal flow was found to significantly influence the flow patterns in the alveolus, resulting in changes to the characteristics of alveolus movement, such as alteration in the area of circulatory flow within the alveolus during respiration and the dispersion of microparticles [16]. In general, the movement of the alveolar wall significantly affects the airflow inside the alveolus [17][18][19]. Meanwhile, the change in flow patterns will alter the movement and distribution of microparticles [20,21]. In addition, much attention has been paid to flow patterns and the movement of microparticles in the bronchial-alveolar system. Tsuda et al. [22] studied the multi-alveolar structure in pulmonary acini, and they found that the deposition of microparticles is strongly affected by the alveolar structural characteristics. Additionally, the complexity of airway structure in pulmonary acini impacts the deposition of microparticles [23,24].
The deposition of microparticles smaller than 20 µm has been generally studied [25], and the effect of Brownian movement to the microparticles, whose diameters are larger than 1 µm, are negligible. The movement and migration of microparticles is mainly controlled by aerodynamics and gravity [26]. Acinar flow shows the characteristics of chaotic flow [27]. Considering the inertial effect of microparticles, it is necessary to propose a model that can fully predict the movement.
So far, a number of studies have been conducted on the deposition of microparticles in the respiratory bronchioles and alveolar ducts. However, the studies of deposition in terminal pulmonary acini are still scarce, and the corresponding rules of deposition need to be further explored. In addition, the airflow in terminal acini is driven entirely by the expansion and contraction of the pulmonary acinar wall, in which the motion of microparticles is different from that in alveolar ducts. The dynamic transport process of microparticles in a closed sac is worth studying.
In order to study the movement and deposition of microparticles in terminal pulmonary acini, we built a respiratory movement model of terminal pulmonary acini using immersed boundary-lattice Boltzmann method (IB-LBM) [28][29][30][31][32][33][34], with an immersed boundary (IB) method describing the pulmonary acinar wall and the morphology of the microparticles. The lattice Boltzmann (LB) method has an advantage in dealing with flows in complex geometries [35,36]. We used an elastic mesh structure that is independent of the flow field to simulate the solid lung tissue, so that the elastic support of the terminal pulmonary acinus wall boundary was set. By changing the internal and external pressure of terminal pulmonary acini, we controlled the movement of the terminal wall to simulate the expansion and contraction process of terminal pulmonary acini. In addition, periodic expansion and contraction of the terminal wall was repeated to simulate multiple respiratory cycles. By building the model, we studied the movement and deposition rules of the microparticles in air flow under different respiratory rates, the diameter of microparticles, and other factors. This provides a deductive research method for evaluating the effects of haze and flying dust on the function of terminal lung acini. From the perspective of medical treatment, it provides a theoretical reference for improving the design of atomized drug administration.
The rest of the paper is organized as follows: Section 2 set up the model. Section 3 verifies the correctness of the model. In Section 4, the deposition characteristics of microparticles are studied in detail. The final summary and conclusion are given in Section 5.

Alveolar Model
In this work, first, we built a model of terminal pulmonary acini that responded to the actual respiratory process (Figure 1). When pulmonary acini are under low pressure, the alveolar ducts are narrowed, and folds characterize the alveolar septums. When pulmonary acini are under high pressure, the alveolar ducts are wide, and the alveolar septums are stretched [37]. Our model can simulate this physiological process. and other factors. This provides a deductive research method for evaluating the effects of haze and flying dust on the function of terminal lung acini. From the perspective of medical treatment, it provides a theoretical reference for improving the design of atomized drug administration. The rest of the paper is organized as follows: Section 2 set up the model. Section 3 verifies the correctness of the model. In Section 4, the deposition characteristics of microparticles are studied in detail. The final summary and conclusion are given in Section 5.

Alveolar Model
In this work, first, we built a model of terminal pulmonary acini that responded to the actual respiratory process (Figure 1). When pulmonary acini are under low pressure, the alveolar ducts are narrowed, and folds characterize the alveolar septums. When pulmonary acini are under high pressure, the alveolar ducts are wide, and the alveolar septums are stretched [37]. Our model can simulate this physiological process. Each healthy human alveolus can be considered a small sphere, approximately 200 μm to 450 μm in diameter. Pulmonary acini are composed of many complex connections of alveoli and pulmonary ducts. In this study, part of the terminal pulmonary acini is selected for numerical modeling, which is the representative area of pulmonary acini. The diameter of the alveolar model is set as 400 μm, and the terminal pulmonary acini is composed of an alveolar duct of 200 μm diameter and surrounding alveoli. The alveolar membrane in this model is a continuous movable wall, whose position is from 0 to 1, as indicated. The respiratory flow curve of a healthy person approximates a sinusoidal curve [38]. The left side shows inlet and outlet of airflow, and the boundary condition is set as = 0.0015 sin 200000 where is the inlet density and t is the time step. Air fluid density in terminal pulmonary acini is a Newtonian fluid with density = 1 and kinematic viscosity = 0.4.
Navier-Stokes equation: The elastic structure of the alveolus is simulated by a series of unstructured grids [39,40]. Each healthy human alveolus can be considered a small sphere, approximately 200 µm to 450 µm in diameter. Pulmonary acini are composed of many complex connections of alveoli and pulmonary ducts. In this study, part of the terminal pulmonary acini is selected for numerical modeling, which is the representative area of pulmonary acini. The diameter of the alveolar model is set as 400 µm, and the terminal pulmonary acini is composed of an alveolar duct of 200 µm diameter and surrounding alveoli. The alveolar membrane in this model is a continuous movable wall, whose position is from 0 to 1, as indicated. The respiratory flow curve of a healthy person approximates a sinusoidal curve [38]. The left side shows inlet and outlet of airflow, and the boundary condition is set as where ρ 1 is the inlet density and t is the time step. Air fluid density ρ 0 in terminal pulmonary acini is a Newtonian fluid with density ρ 0 = 1 and kinematic viscosity nu = 0.4. Fluid dynamics are governed by Equations (2) and (3). Continuity equation: Navier-Stokes equation: The elastic structure of the alveolus is simulated by a series of unstructured grids [39,40]. These grids are directly connected to the alveolar boundary. To simulate the elastic support of the pulmonary tissue, the outer edges of these grids are set to be fixed, which move as the alveolar boundary deforms and generates elastic force at the alveolar boundary. These grids, which move as elastic force and bending force change, are set independently of the flow field. Elastic force is calculated by the equation [41]: where K S1 is the elasticity coefficient of pulmonary tissue, ∆x is the elastic elongation between adjacent grid points, and → n is the unit direction vector. Bending force is computed by the equation: where K b is the bending stiffness and ∂X (s,t) /∂s is the unit tangent vector. The elastic force of the alveolar membrane point is calculated by a neo-Hookean elastic membrane force equation [42]: where K S2 is the elasticity coefficient of the elastic membrane and λ is the stretch ratio of the distance between membrane points. The alveolar movement is generated by changing the inlet density boundary condition, simulating the pressure change in pulmonary acini. The position of the alveolar membrane point changes by the flow field pressure, and the elastic force and bending force change accordingly. Eventually, the alveolus stops moving when reaching a new equilibrium.

Microparticles Model
We studied the movement of micron-sized particles so that gravity is taken into account and the effect of Brown movement is negligible.
Solid microparticles are simulated by the node-spring model, which consists of a central supporting point and some boundary points ( Figure 2). These grids are directly connected to the alveolar boundary. To simulate the elastic support of the pulmonary tissue, the outer edges of these grids are set to be fixed, which move as the alveolar boundary deforms and generates elastic force at the alveolar boundary. These grids, which move as elastic force and bending force change, are set independently of the flow field. Elastic force is calculated by the equation [41]: where K is the elasticity coefficient of pulmonary tissue, ∆ is the elastic elongation between adjacent grid points, and is the unit direction vector. Bending force is computed by the equation: where K is the bending stiffness and ( , ) / is the unit tangent vector. The elastic force of the alveolar membrane point is calculated by a neo-Hookean elastic membrane force equation [42]: where K is the elasticity coefficient of the elastic membrane and is the stretch ratio of the distance between membrane points.
The alveolar movement is generated by changing the inlet density boundary condition, simulating the pressure change in pulmonary acini. The position of the alveolar membrane point changes by the flow field pressure, and the elastic force and bending force change accordingly. Eventually, the alveolus stops moving when reaching a new equilibrium.

Microparticles Model
We studied the movement of micron-sized particles so that gravity is taken into account and the effect of Brown movement is negligible.
Solid microparticles are simulated by the node-spring model, which consists of a central supporting point and some boundary points ( Figure 2). The elastic force between the adjacent grid points, bending force, and elastic force generated by the central support point are applied to each boundary point. Elastic force obeys Hooke's Law, calculated by the equation: where K S3 is tensile rigidity and ∂X (s,t) /∂s is the unit tangent vector. Bending force is calculated by Equation (5).
In addition to elastic force calculated by Equation (7), a central support point is acted on by gravity and inertia force. Gravity G is calculated by the equation: where m is microparticle quality and g is gravity coefficient. Inertia force F k is calculated by the equation: where m s is inertia coefficient, which controls particle inertia. X (s,t) is the center point coordinate vector. When a microparticle approaches the alveolar boundary, if the shortest distance between the center and boundary points is less than three times the microparticle radius, the adhesion mechanism of the alveolar boundary will be activated. Adhesion is simulated by a virtual spring, whose original length is twice the microparticle radius, which connects the central point to the nearest boundary point. Adhesion force F c is calculated by the equation: where k c is the elastic coefficient, l 0 is the expected riveting point distance, l is the actual distance, and R is microparticle radius. When microparticles are close to each other and collide, the repulsive force F r is set according to the momentum conservation rule. The activation distance is set to be twice the microparticle radius, and the limit proximity distance is not less than the microparticle radius. The repulsive force F r is calculated by the equation: where k r is the coefficient of repulsive force.

IB-LBM
When simulating particles in a flow, the lattice Boltzmann and Langevin dynamics (LB-LD) method can implement multi-scale particles movement simulation to study complex flows [43,44]. The IB-LBM can implement deformable boundaries and simulate the interaction between particles and the interaction between particles and the flow field. Additionally, at the micron scale, Brownian motion can be ignored, so, in this study, the IB-LBM is more suitable. The multi-relaxation time (MRT) model has a wider selection of parameters than the lattice Bhatnagar-Gross-Krook (LBGK) model, which is beneficial to propose a whole lung model in the future.
In this work, the discrete lattice Boltzmann method of the multi-relaxation time model is utilized to solve the Navier-Stokes equations. We can find the lattice Boltzmann equation: where g i (x, t) is the distribution function of velocity e i at position x and moment t, ∆t is the time step, S is collision matrix, and g eq i is the equilibrium distribution function. Here, the 9 discrete velocities in 2 dimensions (D2Q9) model is used, and the nine velocities e i are given by the equation: where ∆x is the lattice spacing. The settings of e i enable a strategy to control the migratory directions within a time step. In Equation (12), g eq i is calculated by the equation: where u is fluid velocity and ω i is the weight defined by ω 0 = 4/9, ω i = 1/9 for i = 1 − 4, and ω i = 1/36 for i = 5 − 8. c s = ∆x/ √ 3∆t is the sound speed. According to Lallermand et al. [45], Equation (12) is transformed into the equation: where M is the transformation matrix.Ŝ is the collision matrix of the moment space. In this model, the value ofŜ is given by the equation [46]: where τ is the relaxation time, calculated by the equation [47]: where υ is the kinematic viscosity in Navier-Stokes equations. Once the density distribution functions are known, the macroscopic fluid density, velocity, and pressure can be calculated by equations: where f is the vector of the body force density. The non-equilibrium extrapolation method is used to obtain the particle distribution function on the boundaries, which is marked with g i,b , and calculated by the equation: where g i,n represents the particle distribution function of the neighboring grid and g eq i,b and g eq i,n are, respectively, the equilibrium forms of g i,b and g i,n .
The interaction between the fluid and particles is achieved by the immersed boundary method. In this method, Lagrangian force F is spread onto the collocated grid points near the boundary, which is calculated by the equation: where D(x − X) is the Dirac delta function, calculated by the equation [48]: where n is the dimension and δ(x i − X i ) is calculated by the Equation: Then, the velocity U of the moving boundary X can be calculated by the equation: The position is updated by the equation:

Grid Independence Validation
First, different support grids are set in the same alveolar model to compare the changes in alveolar expansion. By comparing the variation of alveolar expansion, the influence of the supporting grid on alveolar deformation can be identified, which provides a basis for the selection of the scale of the supporting grid. The outer boundary that supports the grid is a free flow boundary but a fixed border for the supporting grid with no displacement of nodes. The verification steps are as follows: 1.
In the initial state, the pressure difference between the alveolar inlet and the free flow boundary is 0. Calculate static pressure P 0 using the equation: where ρ 0 = 1.00 is initial density. The alveoli and supporting mesh are free of stress and deformation; 2.
Change the alveolar inlet pressure by the equation: where ρ in = 1.02 is inlet density. We can find the pressure difference ∆P by the equation: 3.
Record the volume of the alveolar expansion when the alveolar expansion reaches a steady state. Expansion ratio D r is calculated by the equation: where V is the volume of the alveoli after expansion. V 0 is the volume of the alveoli without expansion. We use grid proportion M r to represent the variation of support grid, which is calculated by the equation: where n w is the number of grid points and n m is the number of alveolar membrane points.
In this validation, we set M r as 64, 81, 100, and 121, respectively, with the spatial distribution of the corresponding supporting grid shown in Figure 3a where is the number of grid points and is the number of alveolar membrane points.
In this validation, we set as 64, 81, 100, and 121, respectively, with the spatial distribution of the corresponding supporting grid shown in Figure 3a

points.
In this validation, we set r as 64, 81, 100, and 121, respectively, with the spatial distribution of the corresponding supporting grid shown in Figure 3a Changes in r are shown in Figure 4. The expansion ratio r of r = 100 increases by 0.0065 over the expansion ratio r of r = 81. The expansion ratio r of r = 121 increases by 0.0035 over the expansion ratio r of r = 100. The results indicate that the influence of the supporting grid size on the alveolar deformation tends to converge with the increase in the grid proportion. The alveolar deformation, therefore, is considered not related to the size of the external supporting grid when the grid reaches a certain number.

Validation of Microparticle Inertia Expression Model
The present numerical method is validated by simulating the inertial migration of microparticles in a microchannel (  The central supporting point is initially located at 0.5l from the upper boundary and moves due to gravity , which is proportional to the volume of the microparticle (the volume of the microparticle is represented by microparticle area in the two-dimensional case). When the microparticle moves to l, it continues to travel for some distance without gravity . Microparticles with greater inertia move farther. The test results are shown in Table 1 and Table 2. Table 1 indicates the sliding distances of microparticles with different diameters of D under the same inertia conditions. The larger the microparticle size, the greater the inertia force, so that the microparticle travels farther. Table 2 shows the sliding distances of the same microparticle under different inertial conditions. The sliding distance of the 5 μm particle, which is 11.37 μm, is the longest of the three sizes. Therefore, the 5 μm particles have the maximum inertia force. Moreover, as is 100, the sliding distance is the longest, at 16.08 μm. So, the larger the is, the greater the inertia force is, which makes the microparticle move farther. The above results reveal that the microparticle model can effectively simulate the inertia of the microparticle in the movement process, which meets the requirements for the simulation of the inertia of solid microparticles in the lung.  The central supporting point is initially located at 0.5l from the upper boundary and moves due to gravity G, which is proportional to the volume of the microparticle (the volume of the microparticle is represented by microparticle area in the two-dimensional case). When the microparticle moves to l, it continues to travel for some distance without gravity G. Microparticles with greater inertia move farther.
The test results are shown in Tables 1 and 2. Table 1 indicates the sliding distances of microparticles with different diameters of D under the same inertia conditions. The larger the microparticle size, the greater the inertia force, so that the microparticle travels farther. Table 2 shows the sliding distances of the same microparticle under different inertial conditions. The sliding distance of the 5 µm particle, which is 11.37 µm, is the longest of the three sizes. Therefore, the 5 µm particles have the maximum inertia force. Moreover, as m s is 100, the sliding distance is the longest, at 16.08 µm. So, the larger the m s is, the greater the inertia force is, which makes the microparticle move farther. The above results reveal that the microparticle model can effectively simulate the inertia of the microparticle in the movement process, which meets the requirements for the simulation of the inertia of solid microparticles in the lung.

Respiratory Deformation and Microparticles Movement in Pulmonary Acini
In the continuous respiratory cycle, pulmonary acini show the periodic movement of expansion and contraction, and microparticles entering pulmonary acini move and migrate by a variety of influencing factors. These factors include deformation amplitude, deformation velocity, microparticle mass, gravity direction, and mucus on the pulmonary acini surface. The Reynolds number (Re) of gas in pulmonary acini is very small (Re << 1) so that the airflow can be considered an isothermal incompressible laminar flow. We find the flow in pulmonary acini is relatively smooth with no obvious backflow area by simulation.
The deformation of pulmonary acini during one respiratory cycle is shown in Figure 6. In one respiratory cycle T, the times of inhalation and exhalation are each T/2. t ∈ (nT, nT + T/2) is the inspiratory cycle. When t = nT + T/2, pulmonary acini expand to their maximum volume at a constant positive inlet pressure. t ∈ n + 1 2 T, (n + 1)T is the expiration cycle; when t = (n + 1)T, pulmonary acini shrink to their minimum at the constant negative inlet pressure.

Respiratory Deformation and Microparticles Movement in Pulmonary Acini
In the continuous respiratory cycle, pulmonary acini show the periodic movement of expansion and contraction, and microparticles entering pulmonary acini move and migrate by a variety of influencing factors. These factors include deformation amplitude, deformation velocity, microparticle mass, gravity direction, and mucus on the pulmonary acini surface. The Reynolds number (Re) of gas in pulmonary acini is very small (Re << 1) so that the airflow can be considered an isothermal incompressible laminar flow. We find the flow in pulmonary acini is relatively smooth with no obvious backflow area by simulation.
The deformation of pulmonary acini during one respiratory cycle is shown in Figure 6. In one respiratory cycle T, the times of inhalation and exhalation are each T/2. ∈ ( T， T + T/2) is the inspiratory cycle. When = T + T/2, pulmonary acini expand to their maximum volume at a constant positive inlet pressure. ∈ + T，( + 1)T is the expiration cycle; when = ( + 1)T, pulmonary acini shrink to their minimum at the constant negative inlet pressure. To observe the movement and migration of particles in the process of pulmonary acini expansion and contraction, 160 microparticles of diameter 4 μm are set at the entrance (Figure 7, t = 0), which are subjected to airflow and gravity.

Respiratory Deformation and Microparticles Movement in Pulmonary Acini
In the continuous respiratory cycle, pulmonary acini show the periodic movement of expansion and contraction, and microparticles entering pulmonary acini move and migrate by a variety of influencing factors. These factors include deformation amplitude, deformation velocity, microparticle mass, gravity direction, and mucus on the pulmonary acini surface. The Reynolds number (Re) of gas in pulmonary acini is very small (Re << 1) so that the airflow can be considered an isothermal incompressible laminar flow. We find the flow in pulmonary acini is relatively smooth with no obvious backflow area by simulation.
The deformation of pulmonary acini during one respiratory cycle is shown in Figure 6. In one respiratory cycle T, the times of inhalation and exhalation are each T/2. ∈ ( T， T + T/2) is the inspiratory cycle. When = T + T/2, pulmonary acini expand to their maximum volume at a constant positive inlet pressure. ∈ + T，( + 1)T is the expiration cycle; when = ( + 1)T, pulmonary acini shrink to their minimum at the constant negative inlet pressure. To observe the movement and migration of particles in the process of pulmonary acini expansion and contraction, 160 microparticles of diameter 4 μm are set at the entrance (Figure 7, t = 0), which are subjected to airflow and gravity.  Figure 7 shows the movement and distribution of microparticles in ten respiratory cycles. The groups of microparticles are observed to spread gradually within acini over time, and some microparticles deposit on the acini membrane, which occurs with a certain regu-larity. This indicates that the model we constructed can effectively simulate the movement and distribution of multiple microparticles in pulmonary acini in a proper setting.

Movement and Distribution of Microparticles in Pulmonary Acini
When a microparticle approaches the inner surface of pulmonary acini, it will be captured by mucus, so that its position will not be affected by airflow or gravity. In order to include the trapping effect of mucilage on microparticles, a fixing mechanism is introduced to rivet the microparticles close to the surface of the mucus. A specific method is as follows: When the shortest distance between the point of microparticle boundary and the point of pulmonary acinar boundary is smaller than the diameter of the microparticle, the center of the microparticle is set to connect with the corresponding border point of pulmonary acini by a spring with the original length of the microparticle diameter and spring with the original length of 0.5 µm longer than the length between the center of the microparticle and the corresponding border point of pulmonary acini.
This method better simulates the phenomenon of microparticle trapping on the inner surface. The elastic force is calculated by the equations: where K s,m is the corresponding elastic control coefficient, which keeps the change in the distance between the center of the microparticle and the related point of pulmonary acini boundary within 0.1D. In addition, in order to indicate the number of microparticles captured by mucus, and describe the distribution of microparticles deposition, we counted the number of microparticles near the surface of pulmonary acini according to the following method.
We set the static distance between points on the surface of pulmonary acini as ∆s (the distance between adjacent points when no expansion or contraction occurs), and counted the number of microparticles in the distance band with a width of 30 ∆s inside the boundary of the pulmonary acini (the position of microparticles is represented by the center point). According to the point order of pulmonary acinar boundary, by starting from the left end, we normalized and located the boundary with a set of relative positions from 0 to 1, where 0 refers to the left end and 1 is the right end. The distribution of microparticle deposition on the inner surface of pulmonary acini can be obtained by sampling the number of microparticles with a statistical interval of 0.02.

Effect of Respiratory Depth on Deposition and Distribution of Microparticles
The depth of respiration describes the volume change of the alveolus during a single respiratory cycle. We change the inlet pressure to adjust the alveolar expansion (contraction) rate and define the alveolar deformation rate by the equation: where ∆s p is the change of alveolar volume (alveolar area in the two-dimensional (2D) case), which is calculated by the equation: where s p,t is the alveolar volume at t time and s p is the equilibrium volume (area in 2D case) of the same internal and external alveolar pressure. The increase in the alveolar deformation rate represents the enhancement of respiratory depth. We set the variation range of alveolar deformation rate during normal respiration as (0~11.97%), and deep respiration as (0~29.43%).
We use microparticles with a diameter D of 4 µm to simulate the deposition progress, of which the distribution state of initial position at t = 0 is shown in Figure 7. In order to study the movement and distribution of microparticles in different gravity conditions, we use θ to represent the direction of gravity, which is set as 0 • , 30 • , 60 • and 90 • . Different gravity directions simulate various body position changes, which lead to direction changes between acinar and gravity. After 10 respiratory cycles, most of the microparticles are close to the alveolar wall, and we calculate the microparticle distribution and deposition (Figure 8).
where , is the alveolar volume at t time and is the equilibrium volume (area in 2D case) of the same internal and external alveolar pressure. The increase in the alveolar deformation rate represents the enhancement of respiratory depth. We set the variation range of alveolar deformation rate during normal respiration as (0~11.97%), and deep respiration as (0~29.43%).
We use microparticles with a diameter D of 4 μm to simulate the deposition progress, of which the distribution state of initial position at t = 0 is shown in Figure 7. In order to study the movement and distribution of microparticles in different gravity conditions, we use θ to represent the direction of gravity, which is set as 0°, 30°, 60° and 90°. Different gravity directions simulate various body position changes, which lead to direction changes between acinar and gravity. After 10 respiratory cycles, most of the microparticles are close to the alveolar wall, and we calculate the microparticle distribution and deposition ( Figure 8). In Figure 8a-d, Case 1 represents the process of the natural sedimentation of the microparticles, when no alveolar movement (no respiration) occurs, which is shown as the basic reference. As we can see in Figure 8a, more microparticles distribute at the inlet of alveolus when not breathing. Case 2 represents the normal respiratory process, and more microparticles distribute in the middle alveolus compared with Case 1. Case 3 reveals the deep respiratory process, in which more microparticles distribute to the distal alveolus. In Figure 8a-d, Case 1 represents the process of the natural sedimentation of the microparticles, when no alveolar movement (no respiration) occurs, which is shown as the basic reference. As we can see in Figure 8a, more microparticles distribute at the inlet of alveolus when not breathing. Case 2 represents the normal respiratory process, and more microparticles distribute in the middle alveolus compared with Case 1. Case 3 reveals the deep respiratory process, in which more microparticles distribute to the distal alveolus. We find that, under the same gravity condition, the higher the alveolar deformation rate is, the more microparticles distribute in the distal alveolar wall. It suggests that deep respiration is conducive to the movement and deposition of microparticles in the deeper region of pulmonary acini. Figure 8b-d shows the effect of gravity. In Case 3, microparticles distribute mostly in the deep area of pulmonary acini, followed by Case 1 and Case 2. Therefore, the results of this part of the study suggest that deep breath facilitates the migration and deposition of microparticles to the deeper area in pulmonary acini.

Effect of Microparticle Diameter on Deposition and Distribution of Microparticles
The diameter of microparticles affects their movement and distribution by affecting gravity and its morphology. We use microparticles with diameters D of 2, 3, 4, and 5 µm to simulate the deposition progress, of which the distribution state of initial position at t = 0 is shown in Figure 7. We observe the effect of different microparticle diameters on their deposition during ten deep respiration cycles. Similarly, we use θ to represent the gravity direction, which is set as 0 • , 30 • , 60 • and 90 • . We calculate the microparticle distribution and deposition in Figure 9.
is, the more microparticles distribute in the distal alveolar wall. It suggests that deep respiration is conducive to the movement and deposition of microparticles in the deeper region of pulmonary acini. Figure 8b-d shows the effect of gravity. In Case 3, microparticles distribute mostly in the deep area of pulmonary acini, followed by Case 1 and Case 2. Therefore, the results of this part of the study suggest that deep breath facilitates the migration and deposition of microparticles to the deeper area in pulmonary acini.

Effect of Microparticle Diameter on Deposition and Distribution of Microparticles
The diameter of microparticles affects their movement and distribution by affecting gravity and its morphology. We use microparticles with diameters D of 2, 3, 4, and 5 μm to simulate the deposition progress, of which the distribution state of initial position at t = 0 is shown in Figure 7. We observe the effect of different microparticle diameters on their deposition during ten deep respiration cycles. Similarly, we use θ to represent the gravity direction, which is set as 0°, 30°, 60°and 90°. We calculate the microparticle distribution and deposition in Figure 9. The distribution of microparticles with different diameters in acini is shown in Figure 9a when gravity direction θ is 0°. We can observe that fewer microparticles of 2 μm distribute in the middle alveolus. The distribution of 3 μm microparticles is similar to those of 2 μm, yet the number of microparticles is larger. Some microparticles of 4 μm and 5 μm distribute near the distal wall of pulmonary acini. In general, larger particles are more likely to The distribution of microparticles with different diameters in acini is shown in Figure 9a when gravity direction θ is 0 • . We can observe that fewer microparticles of 2 µm distribute in the middle alveolus. The distribution of 3 µm microparticles is similar to those of 2 µm, yet the number of microparticles is larger. Some microparticles of 4 µm and 5 µm distribute near the distal wall of pulmonary acini. In general, larger particles are more likely to deposit in the alveolar wall at the same respiratory cycle and depth. The results under four gravity conditions (Figure 9a-d) reveal that the direction of gravity is the main reason that affects the deposition location of microparticles, which is more obviously observed in larger microparticles. Our results, therefore, indicate that the pulmonary spatial orientation and the diameter of microparticles strongly impact the efficacy of atomization drug delivery.

Effect of Respiratory Rate on Deposition and Distribution of Microparticles
The respiratory rate describes the rate of volume change of the alveolus during a single respiratory cycle. We set the normal rate in the deep respiratory process as U r . In the slow respiratory process, the respiratory rate is 0.5U r , while the respiratory rate is 2U r in the rapid process. The distribution state of the initial position at t = 0 is shown in Figure 7. We use microparticles with a diameter D of 4 µm to simulate the deposition progress, and we observe the effect of different respiratory rates and different gravity directions on microparticle deposition in a deep breath during ten respiratory cycles. The results are shown in Figure 10.
four gravity conditions (Figure 9a-d) reveal that the direction of gravity is the main reason that affects the deposition location of microparticles, which is more obviously observed in larger microparticles. Our results, therefore, indicate that the pulmonary spatial orientation and the diameter of microparticles strongly impact the efficacy of atomization drug delivery.

Effect of Respiratory Rate on Deposition and Distribution of Microparticles
The respiratory rate describes the rate of volume change of the alveolus during a single respiratory cycle. We set the normal rate in the deep respiratory process as . In the slow respiratory process, the respiratory rate is 0.5 , while the respiratory rate is 2 in the rapid process. The distribution state of the initial position at t = 0 is shown in Figure  7. We use microparticles with a diameter D of 4 μm to simulate the deposition progress, and we observe the effect of different respiratory rates and different gravity directions on microparticle deposition in a deep breath during ten respiratory cycles. The results are shown in Figure 10. The distribution of microparticles in acini during different respiratory rates is shown in Figure 10a when gravity direction θ is set to 0°. It can be observed that, within the normal range of respiratory rates, the rate has no obvious effect on the microparticle deposition or distribution in pulmonary acini, and the microparticle deposition pattern almost remains unchanged. The experimental results under four gravity directions show that the The distribution of microparticles in acini during different respiratory rates is shown in Figure 10a when gravity direction θ is set to 0 • . It can be observed that, within the normal range of respiratory rates, the rate has no obvious effect on the microparticle deposition or distribution in pulmonary acini, and the microparticle deposition pattern almost remains unchanged. The experimental results under four gravity directions show that the distribution of microparticles in acini is approximately the same. Therefore, we believe that respiration rate has relatively little influence on distribution and deposition in pulmonary acini.
In low Reynolds number cases, the increase in respiratory rate does not affect the flow in pulmonary acini obviously, but only changes the flow velocity. This change will cause some effect on the inertia of particles. However, for particles of 2-5 µm, the influence is limited enough.

Conclusions and Discussion
In this study, we used IB-LBM to simulate the movement and deposition of solid (liquid or hydrogel) particles numerically in pulmonary acini. We proposed a method based on the elastic mesh support outside the flow field to propose the elastic lung tissue to support the alveolus, and to simulate the alveolar movement during respiration by changing the pressure at the alveolar entrance to promote alveolus expansion and contraction. On this basis, we studied the movement and distribution law of microparticle populations with different microparticle diameters, respiratory depths, respiration rates, and gravity directions. The main conclusions are as follows: (1) Large breathing depth is conducive to the movement and deposition of microparticles in pulmonary acini. (2) The microparticles with a larger diameter are easier to deposit after entering the alveoli. (3) The direction of gravity has a powerful effect on the deposition location of the microparticle. (4) In the same respiratory depth and gravity direction, the respiratory rate has little effect on the distribution of microparticles.
The research conclusions above provide a certain reference for designing aerosol delivery and the reduction in particle pollution exposure. For example, in the design of aerosol delivery, to achieve the best efficacy of an aerosol, the number of drug particles should be maximized in alveoli. To evenly distribute the microparticles, deep breath and change in body position should be recommended when administrating aerosol. In addition, to limit the impacts of haze pollution on health, accelerating respiratory rate is more effective than deep respiration for inhibiting the deposition of microparticles in the lungs.
There are also some limitations to this study. First, we used a two-dimensional model to propose pulmonary acini and microparticle models. However, the distribution and deposition of microparticles are three-dimensional, so what the model reveals is from a limited perspective. Second, this model is limited by the grid resolution; thus, it does not provide the possibility of studying microparticles with tiny diameters (less than 1 µm), which are liable to enter the pulmonary acini. Therefore, to deepen the corresponding research, the current model needs to be further improved.