**4.2. Lattice Boltzmann method (LBM) for aerosol transport in alveoli**

A portion of the inhaled aerosol will enter the respiratory zone (G17–G23) of human lung, which is mainly the alveolar region. In this case, the difficulty in using CFPD model is that the point-mass and point-force assumptions do not hold and particle-particle interactions are not negligible. Lattice Boltzmann method (LBM) can be employed for the two-way fluid-particle interaction and many-particle interaction in alveoli, due to its advantage in resolving particles that are comparative or larger than the mesh cells in size.

## *4.2.1. Background and literature review*

visible via CT/MRI imaging. Accordingly, two methods can be further developed: (1) coupling methods with single-path or multiple-path lung airway trees and (2) virtual lung model with idealized airways for higher generations covering the entire lung conducting zone (G0–G17). This is crucial to create a full-scale lung aerosol dynamics simulation featuring full breathing

Coupling methods combine coupled boundary conditions with truncated lung airway trees to represent the whole conducting zone. Realization of coupling methods can be performed in

**1.** Divide the whole lung airway trees into segments [40] or truncate the airway tree with

**3.** The boundary conditions applied to each truncated airway ends will be obtained from the previous simulation of the next larger section (inhalation) or the next smaller section

Kleinstreuer and Zhang [40] used the TBU of Weibel type A and constructed the symmetric and in-plane tracheobronchial airways from G0 to G15, and used coupling method to simulate airflow and particle deposition in the lung airways with steady-state inhalation flow rate *Q*in = 30 L/min. Walters and Luke [42] created an 8-generation airway tree with random connection angles following Weibel's lung morphology. They applied coupled pressure boundary conditions for inhalation cases, and validated the coupling method for inhalation cases. They extended the application of this coupling method on particle transport in lung airways [43]. Tian et al. [41] developed the stochastic individual path (SIP) and multiple stochastic individual path (MSIP) approaches to represent the whole-lung geometry. Compared to the coupling method applied on detached airway geometries [41], this model is able to simulate transient

Since modeling the whole lung is computationally expensive, using the coupling method is potentially advantageous for fast-solution purposes. However, accuracy of the coupling methods still needs to be tested compared with using virtual lung model containing the entire

To encompass the entire conducting zone, the main objective of virtual lung model reconstruction is to develop a stochastic algorithm that generates 3D asymmetric human lower lung airways, which can then be connected to the subject-specific upper lung airways [46–48]. To establish such a virtual lung model, the upper airway (i.e., mouth and nasal airways, throat, and trachea) will be exactly reconstructed using image processing software that converts CT/ MRI images. Conducting airways, from the trachea down to the level of the terminal bronchioles, will be generated via a stochastic algorithm [46–48]. Specifically, by generating

breathing scenarios. Recently, similar models have also been developed [44, 45].

*4.1.1. Coupling methods with single-path or multiple-path lung airway trees*

single path or multiple paths left [41–43];

**2.** Simulate each segment in series or in parallel;

conducting zone without any truncations.

*4.1.2. Virtual lung model for the entire conducting zone*

cycles.

the following steps:

70 Aerosols - Science and Case Studies

(exhalation).

Lattice Boltzmann method (LBM) solves the Lattice Boltzmann equation (LBE) to obtain the probability density distribution, which is then used to obtain the pressure and velocity field through integration. Details of the LBE, discrete speeds, discretization, and boundary conditions are given in [49] and avoided here for brevity. LBM consists of two main steps: (1) propagation: particle distribution moves along the lattice bonds to neighboring lattice nodes and (2) collision: particles on the same node interact and adjust velocities to conserve momentum and mass. Communication between nodes in the lattice is required only in the propagation step requiring message passing or update of the shared memory location.

LBM has been applied to the airflow and particle transport simulations in the human airways due to the advantages mentioned above compared with CFPD models. Many studies focused on the flow structures in nasal cavities. For example, Finck et al. [50] first employed LBM for steady simulations of laminar nasal flow. The model predictions agree well with those of finitevolume based model using structured grids. However, the LBM has high flexibility in fast grid generation and easy boundary-condition implementation for such complex geometries compared to conventional Navier-Stokes solvers. Hörschler et al. [51] compared steady and unsteady flows in real nasal cavity, and concluded that the transient effects are only important at Re < 1500. Henn et al. [52] coupled Lagrangian particle tracking method to LBM to simulate aerosol dynamics in a human nasal cavity. They concluded that transient simulations are necessary to study particle dynamics in nasal cavities. Other studies modeled airflow in the nasal cavities with diagnosed pathologies [53].

There are also existing researches focusing on air flow simulation in human upper airway using LBM. Ball et al. [54] showed that LBM yields higher fidelity than the RANS models for laminarto-turbulent transition flow regime predictions inside an idealized human extrathoracic airway, as compared to experimental results. Schroeder et al. [55] presented LB simulation results of pulsatile flow at two breathing frequencies in a human respiratory system down to G6. They showed that inspiratory flow is characterized by flow separation and secondary flow, while expiratory flow is characterized by more homogeneous structure but higher vorticity. This is in contrast with the study by Eite et al. [56], which showed that less vorticity is generated at exhalation than at inhalation in a nasal cavity. Krause [53] investigated respirations in the upper part of the human lungs using LBM coupled with a model for the rest of lung, thereby constructing a two-scale whole lung model. This latter model assumed that the lower part of the bronchial tree can be described by a symmetrical dyadic tree of pipes, where Poiseuille type flow holds. Imai et al. [57] developed a LBM with GPU acceleration using an adaptive meshing method. The model was used to study the effect of breath holding on the deposition of microparticles in pulmonary airways during inhalation. The subject-specific airway model contains thirteen generations of airways. The study found that breath holding can effectively increase particle deposition; yet the effect varies for different particle sizes. Wang and Elghobashi [58] analyzed flow structures in human upper airways including the nasal cavity, the pharynx, the larynx, and the upper part of the trachea for both healthy and diseased subjects. The laminar-transitional-turbulent flow in the airway was captured in inspiration, and laminar flow was found during expiration. The authors proposed that the maximum positive value of the local time-averaged second derivative of pressure can be used to locate the obstructions due to disease in the upper airway. It should be noted that most studies claimed that LBM provides better accuracy for air flow simulation in respiratory systems. However, Chen and Gutmark [59] did a comparison study of airflow in an idealized human extrathoracic airway. They claimed that LES is more accurate in predicting root-mean-square flow field compared with Reynolds stress model and LBM.

In contrast to many existing LBM simulations on airflow patterns in human respiratory systems, there are only a few studies using LBM for lung aerosol dynamics. Henn et al. [52] used LBM coupled with Lagrangian tracking method to study microparticle transport in human nasal cavities. Trunk et al. [60] developed an Euler-Euler LBM method to investigate the particle dynamics in an idealized bifurcation geometry representing the split from trachea to main bronchi. Particle capturing BC and numerical stabilization were discussed for lung aerosol dynamics simulations. Li and Kleinstreuer [61] analyzed transient airflow in twodimensional alveoli and bifurcating alveolar ducts. The results highlighted the importance of alveolar geometries and shapes, as well as the effect of the moving alveolus walls. Implications of these effects on particle deposition in the alveolar sac were discussed. However, simulation of particle suspension has not been performed.

## *4.2.2. LBM vs. CFPD models*

necessary to study particle dynamics in nasal cavities. Other studies modeled airflow in the

There are also existing researches focusing on air flow simulation in human upper airway using LBM. Ball et al. [54] showed that LBM yields higher fidelity than the RANS models for laminarto-turbulent transition flow regime predictions inside an idealized human extrathoracic airway, as compared to experimental results. Schroeder et al. [55] presented LB simulation results of pulsatile flow at two breathing frequencies in a human respiratory system down to G6. They showed that inspiratory flow is characterized by flow separation and secondary flow, while expiratory flow is characterized by more homogeneous structure but higher vorticity. This is in contrast with the study by Eite et al. [56], which showed that less vorticity is generated at exhalation than at inhalation in a nasal cavity. Krause [53] investigated respirations in the upper part of the human lungs using LBM coupled with a model for the rest of lung, thereby constructing a two-scale whole lung model. This latter model assumed that the lower part of the bronchial tree can be described by a symmetrical dyadic tree of pipes, where Poiseuille type flow holds. Imai et al. [57] developed a LBM with GPU acceleration using an adaptive meshing method. The model was used to study the effect of breath holding on the deposition of microparticles in pulmonary airways during inhalation. The subject-specific airway model contains thirteen generations of airways. The study found that breath holding can effectively increase particle deposition; yet the effect varies for different particle sizes. Wang and Elghobashi [58] analyzed flow structures in human upper airways including the nasal cavity, the pharynx, the larynx, and the upper part of the trachea for both healthy and diseased subjects. The laminar-transitional-turbulent flow in the airway was captured in inspiration, and laminar flow was found during expiration. The authors proposed that the maximum positive value of the local time-averaged second derivative of pressure can be used to locate the obstructions due to disease in the upper airway. It should be noted that most studies claimed that LBM provides better accuracy for air flow simulation in respiratory systems. However, Chen and Gutmark [59] did a comparison study of airflow in an idealized human extrathoracic airway. They claimed that LES is more accurate in predicting root-mean-square flow field compared

In contrast to many existing LBM simulations on airflow patterns in human respiratory systems, there are only a few studies using LBM for lung aerosol dynamics. Henn et al. [52] used LBM coupled with Lagrangian tracking method to study microparticle transport in human nasal cavities. Trunk et al. [60] developed an Euler-Euler LBM method to investigate the particle dynamics in an idealized bifurcation geometry representing the split from trachea to main bronchi. Particle capturing BC and numerical stabilization were discussed for lung aerosol dynamics simulations. Li and Kleinstreuer [61] analyzed transient airflow in twodimensional alveoli and bifurcating alveolar ducts. The results highlighted the importance of alveolar geometries and shapes, as well as the effect of the moving alveolus walls. Implications of these effects on particle deposition in the alveolar sac were discussed. However, simulation

nasal cavities with diagnosed pathologies [53].

72 Aerosols - Science and Case Studies

with Reynolds stress model and LBM.

of particle suspension has not been performed.

As compared with the CFPD models, LBM has advantages in the following aspects [54]: the convection operator is linear in LBM; the conservation laws are exactly preserved in LBM; the calculations are arithmetic operations in LBM; and the boundary conditions (BCs) like no-slip or moving boundaries are easy to implement. In addition, the entropic lattice Boltzmann method is attractive in modeling transitional and fully turbulent flows for high-Re simulations [62]. It should be noted that when the particle size is much larger than the lattice size, LBM can readily account for the two-way fluid-particle interaction and many body interaction by accurately calculating the hydrodynamic force and torque exerted on the particle by the fluids [63–65]. Also, Brownian motion of particles arises spontaneously from stochastic fluctuations in the fluid stress tensor, avoiding application of a random force on particles. Therefore, LBM is suitable for the simulation of aerosol dynamics, especially in the alveoli sacs. Furthermore, LBM can also be used to simulate the mucociliary clearance process when coupled with immersed boundary method or finite difference method [66, 67].

In conclusion, LBM is a rapidly developing method in lung aerosol dynamics simulation. Many researches have been done concerning airflow dynamics. On the other hand, LBM simulation including particle transport and deposition in alveoli is still largely unexplored. Much more efforts are demanded to couple particle transport with airflow simulation using LBM, with both sub-grid-size particles, and macro-size particles.

## **4.3. Physiologically based pharmacokinetic (PBPK) models for aerosol translocation into systemic regions**

### *4.3.1. Background and literature review*

Carried by ultrafine aerosols, toxicants and therapeutic components can easily translocate from lung deposition sites into the systemic circulation, and accumulate in extrapulmonary organs. Therefore, after obtaining deposited dose of environmental pollutants or pulmonary drugs in the human respiratory systems, the judgment of potential adverse effects or curing effectiveness depends on agents' toxicokinetic/pharmacokinetic information with more insights into their in-vivo behavior, and then toxicology, i.e., absorption, distribution, metabolism, and elimination. In order to provide comprehensive regional data for risk assessment, it is necessary to link a whole-body PBPK model to CFPD models, which extends exposure and deposition modeling to health endpoints. The integration of preclinical data generated by CFPD-PBPK models is also a key to optimize the drug delivery processes for clinical practice.

Physiologically based pharmacokinetic (PBPK) models are based on differential equations that are designed to simulate the absorption, distribution, metabolism, and excretion (ADME) in compartments which are connected biological regions, usually in whole animals or humans [11, 66]. Realistic descriptions of physiology and relevant pathways of metabolism can be incorporated into PBPK models. The history and details of PBPK models have been well discussed in [67, 68].

Focusing on the development of hybrid CFPD-PBPK models, the respiratory system serves as one compartment, but has significant influence on the whole body as it is the main barrier between inhaled aerosols and systemic circulation (see **Figure 11**). As direct inputs to the PBPK model, regional airway-tissue doses must be accurately calculated from the lung deposited doses acquired from the CFPD simulations. Therefore, the aim of incorporating CFPD models with PBPK models is to develop a new boundary condition for site-specific airway-tissue dosimetry for PBPK modeling that includes essential mechanisms, such as clearance, tissue reaction, tissue metabolism, diffusion, and blood perfusion. Assuming a linear first-order constitutive relationship, the material conservation law for compartment *i* of the airway tissue can be expressed as:

$$\left\|\mathbf{\hat{c}c}\_{i} \mid \mathbf{\hat{c}t} = \sum\_{j=1}^{N} (k\_{\boldsymbol{y}}\mathbf{c}\_{j} - k\_{\boldsymbol{y}}\mathbf{c}\_{i}) + S\_{i}(t) \right\| \tag{21}$$

where *ci* is the amount of pathogen per unit volume accumulating in compartment *i*, *kij* is the transfer rate from compartment *j* to *i*, and *Sij* is a general source term representing exterior inputs, reactions, clearance, etc. The transfer rate, *kij*, can be determined from existing experimental data, or in vivo deposition/retention and mass transfer data obtained from laboratory animals. Other biological material properties and site-specific boundary conditions for the PBPK model will be obtained from the open literature.

Efforts have been made to link CFPD and PBPK/PBTK models for simulations of inhaled drugs or toxicants, from deposition to systemic regions of interests [69–75]. For example, Rigg et al. [69] established a numerical model for particle dissolution, absorption, and clearance (DAC) in the nasal cavity. Their DAC model can estimate the tissue dose for particles after deposition, which are the more accurate inputs than deposited dose for PBPK/PBTK calculation. Using subject-specific respiratory system geometries of different species, Corley et al. [70] compared the parameter difference between humans and animals and adapted the nasal CFD-PBPK model [71] to perform interspecies variability analyses and sensitivity analyses for vapor exposure. Anatomic, physiological, and biochemical parameter values were also listed [70]. Parameters such as airway geometry, inhalation flow rates, vapor concentration, partition coefficient between air and tissue, tissue thickness, and metabolism rate were identified as the significant factors that influence local vapor uptakes. Recently, they extended the CFD-PBPK model to transient simulation for toxic vapors [72]. Furthermore, PBPK models with different compartmental divisions have been established to estimate pulmonary drug delivery [73, 74]. Parameter values were also provided in the context.

#### *4.3.2. Comparisons among CFPD-PBPK, CFPD, and PBPK models*

Compared with PBPK models, CFPD-PBPK models have the advantage in considering the subject variabilities induced by interindividual morphologies of human respiratory systems, which can be reflected by the different local deposition patterns provided by CFPD calculations. Additionally, CFPD-PBPK models extend exposure and deposition modeling of CFPD to health endpoints, i.e., delivered dosage estimations at site of interests in the whole human body.

Although research efforts have been made to develop hybrid CFPD-PBPK models as discussed in Section 4.3.1, none of the existing models display the realism needed for accurate dosimetry calculations, because they ignored the influence of essential physical and chemical mechanisms such as hygroscopic characteristics on the dispersion and deposition of droplets in the respiratory system, and used incomplete human lung airway geometries.

Therefore, future directions to develop more accurate and computationally efficient CFPD-PBPK models are to make more realistic compartment divisions, mechanism selections, and parameter value determinations to reflect the underlying aerosol dynamics in systemic regions [76]. For example, the compartmental division of cellular-level PBPK model should be based on specific anatomy and mass transfer characteristics of different lung regions. It is also important to collaborate closely with experimentalists and clinical doctors to find reliable data sets reflecting responses and pharmacokinetics for evaluating the performance of the CFPD-PBPK model [77]. Moreover, the determination of human physiological and biochemical parameters requires strong collaborations. These parameters include breathing rate, skin surface area, cardiac output, tissue blood flow rates, tissue volumes, partition coefficients, and biochemical rate constants [68]. Subject-variability must be considered too, i.e., physiological and biochemical parameter differences due to the differences in interindividual response to the biological effects of inhaled aerosols. It will not only contribute to the evaluation of the confidence in PBPK models on the lack of precise knowledge of the parameter values, but also contribute to the uncertainty analyses for the extrapolations between animals and humans. The error source must be identified and reduced to enhance the accuracy of the numerical results [75].

#### **4.4. Others**

Focusing on the development of hybrid CFPD-PBPK models, the respiratory system serves as one compartment, but has significant influence on the whole body as it is the main barrier between inhaled aerosols and systemic circulation (see **Figure 11**). As direct inputs to the PBPK model, regional airway-tissue doses must be accurately calculated from the lung deposited doses acquired from the CFPD simulations. Therefore, the aim of incorporating CFPD models with PBPK models is to develop a new boundary condition for site-specific airway-tissue dosimetry for PBPK modeling that includes essential mechanisms, such as clearance, tissue reaction, tissue metabolism, diffusion, and blood perfusion. Assuming a linear first-order constitutive relationship, the material conservation law for compartment *i* of the airway tissue

1

= ¶ ¶= - + å *N*

*j*

PBPK model will be obtained from the open literature.

Parameter values were also provided in the context.

*4.3.2. Comparisons among CFPD-PBPK, CFPD, and PBPK models*

/ ( ) ()

where *ci* is the amount of pathogen per unit volume accumulating in compartment *i*, *kij* is the transfer rate from compartment *j* to *i*, and *Sij* is a general source term representing exterior inputs, reactions, clearance, etc. The transfer rate, *kij*, can be determined from existing experimental data, or in vivo deposition/retention and mass transfer data obtained from laboratory animals. Other biological material properties and site-specific boundary conditions for the

Efforts have been made to link CFPD and PBPK/PBTK models for simulations of inhaled drugs or toxicants, from deposition to systemic regions of interests [69–75]. For example, Rigg et al. [69] established a numerical model for particle dissolution, absorption, and clearance (DAC) in the nasal cavity. Their DAC model can estimate the tissue dose for particles after deposition, which are the more accurate inputs than deposited dose for PBPK/PBTK calculation. Using subject-specific respiratory system geometries of different species, Corley et al. [70] compared the parameter difference between humans and animals and adapted the nasal CFD-PBPK model [71] to perform interspecies variability analyses and sensitivity analyses for vapor exposure. Anatomic, physiological, and biochemical parameter values were also listed [70]. Parameters such as airway geometry, inhalation flow rates, vapor concentration, partition coefficient between air and tissue, tissue thickness, and metabolism rate were identified as the significant factors that influence local vapor uptakes. Recently, they extended the CFD-PBPK model to transient simulation for toxic vapors [72]. Furthermore, PBPK models with different compartmental divisions have been established to estimate pulmonary drug delivery [73, 74].

Compared with PBPK models, CFPD-PBPK models have the advantage in considering the subject variabilities induced by interindividual morphologies of human respiratory systems, which can be reflected by the different local deposition patterns provided by CFPD calculations. Additionally, CFPD-PBPK models extend exposure and deposition modeling of CFPD

*c t kc k c S t* (21)

*i ij j ji i i*

can be expressed as:

74 Aerosols - Science and Case Studies

Other numerical methods that can be integrated or employed to improve or replace current CFPD models are direct numerical simulations (DNSs) for turbulence [78], complete numerical simulations (CNSs) for multiphase flow [79], discrete element methods (DEMs) [80], etc. They will be able to extend the CFPD model capabilities by encompassing more underlying physics (e.g., particle-particle interactions) and less simplifications (e.g., point-mass and point-force assumptions).
