Nanoparticle Diffusion in Respiratory Mucus Influenced by Mucociliary Clearance: A Review of Mathematical Modeling

Abstract

Background:

Inhalation and deposition of particles in human airways have attracted considerable attention due to importance of particulate pollutants, transmission of infectious diseases, and therapeutic delivery of drugs at targeted areas. We summarize current state-of-the art research in particle deposition on airway surface liquid (ASL) influenced by mucociliary clearance (MCC) by identifying areas that need further investigation.

Methodology:

We aim to review focus on governing and constitutive equations describing MCC geometry followed by description of mathematical modeling of ciliary forces, mucus rheology properties, and numerical approaches to solve modified time-dependent Navier–Stokes equations. We also review mathematical modeling of particle deposition in ASL influenced by MCC, particle transport in ASL in terms of Eulerian and Lagrangian approaches, and discuss the corresponding mass transport issues in this layer. Whenever required, numerical predictions are contrasted with the pertinent experimental data.

Results:

Results indicate that mean mucus and periciliary liquid velocities are strongly influenced by mucus rheological characteristics as well as ciliary abnormalities. However, most of the currently available literature on mucus fiber spacing, ciliary beat frequency, and particle surface chemistry is based on particle deposition on ASL by considering a fixed value of ASL velocity. The effects of real ASL flow regimes on particle deposition in this layer are limited. In addition, no other study is available on modeling nonhomogeneous and viscoelastic characteristics of mucus layer on ASL drug delivery.

Conclusion:

Simplification of assumptions on governing equations of drug delivery in ASL influenced by MCC leads to imposing some limitations on numerical results.

Introduction

Mucus

Mucus is a slippery hydrogel fluid substance produced by mucosal tissues of human body, including airways,^1,2 gastrointestinal tract,^3,4 genitourinary tract,^5,6 and ocular membrane.⁷ Airway surface fluid (ASL), which contains mucus, is a thin fluid layer that covers the conducting airways of the respiratory tract and is essential for tissue lubrication, homeostasis of the hydrated epithelial layer, and nutrition of the underlying epithelium. ASL is also a first-line defense against pathogens and foreign substances.^8,9 Mucus is mainly composed of cross-linked and entangled mucin fibers and is secreted by goblet cells and submucosal glands.^10–12 Its branched chains create a three-dimensional (3D), covalently linked mesh-like network (Fig. 1).^13,14

FIG. 1.

Respiratory mucus sample image using scanning electron microscopy shows that a mucus sample has a meshwork architecture (Scale bar shows 500 nm). Reproduced with permission from Schuster et al.¹⁴

Mucosal layer is an important protective barrier against inhaling pathogens and other particles that enter this hydrogel matrix.¹⁵ Mucus is mainly composed of water (∼95%), mucins (∼0.2%–5.0%), globular proteins (∼0.5%), salts (∼0.5%–1.0%), lipids (1%–2%), deoxyribonucleic acid (DNA), epithelial cells, and cellular debris, which can make a dense viscoelastic layer over epithelial cells as a selective barrier to drugs and other particles.¹⁶

Pore size (mesh network size), viscoelasticity, pH, ionic strength, and charge are the main physicochemical key parameters influencing drug delivery across mucus barrier.¹⁶ For instance, mucus viscoelasticity, which significantly affects the mucus filtering mechanism, creates a mesh network filter that reduces particle penetration and diffusion rates.^17,18 Complex non-Newtonian and viscoelastic mucus behavior results from viscous (flow) and elastic (resistance to deformation) components.⁹ Studies show that mucus behaves as a shear thinning fluid, and at low shear rates, viscosity of the mucus is 100–10,000 times higher than water. In contrast, at higher shear rates, mucus has a viscosity comparable to water.¹⁹

In addition to its protective properties, mucosal surface shows promise in delivery of aerosolized medicines in diseases such as chronic obstructive pulmonary disease (COPD), asthma, cystic fibrosis (CF), and lung cancer,^20–23 even though passage of nanoparticles of a drug through the barrier can be technically challenging because of its protective nature. There are studies on the effect of particle size and surface chemistry on mobility of nanoparticles in respiratory mucus in healthy state^14,24–28 or diseases such as CF,^29–31 chronic rhinosinusitis,³² COPD,³³ and lung cancer.³⁴

Newby et al. showed that 100 nm and smaller drug-loaded nanoparticles, as well as pure drug powders or drug solutions, can easily diffuse through the mucosal layer by Brownian motion, and mucosal mesh does not have a great impact on the motion of these particles in mucus. However, there are challenges in mucosal drug delivery with larger drug-loaded nanoparticles that are comparable to, or larger than, mucus pore size (∼200 nm) where diffusion is not by simple Brownian motion and is typically slower. In this regard, advanced particle tracking of drug-loaded nanoparticles or computational modeling is necessary to estimate particle trajectories (Fig. 2).³⁵

FIG. 2.

Schematic illustration of 3D mucus meshwork structure and effect of particle size on transport in mucus layer by Brownian and non-Brownian motion. Blocks on the left show epithelium and ciliated cells. Reproduced with permission from Newby et al.³⁵ 3D, three dimensional.

ASL and mucociliary clearance

Tracheobronchial tree of lung consists of a conducting region formed by the trachea and tracheal bronchi and a gas exchange region formed by respiratory bronchioles, alveolar ducts, and finally, alveoli, where gas exchange between the atmosphere and capillaries occurs. Trachea and bronchi are covered by ASL, which is produced by a variety of cells, including ∼50% by ciliated cells, 12% by secretory cells (mostly goblet), 20%–30% by basal cells, and some by undifferentiated cells.^36,37 As the airways continue to branch out into bronchioles and decrease in diameter, there is a decrease in the density of mucus-producing and basal goblet cells, while goblet cells, which produce aqueous secretions, are more common in small airways; however, the density of ciliated cells remains almost constant in various generations of bronchioles.³⁸ Due to shear thinning behavior of ASL, this multiphase fluid is reported to have two main sublayers:

(a) Aqueous periciliary liquid (PCL) layer: it has a depth of about 5–10 μm in upper airways of bronchi, trachea, and nasal cavity,³⁹ and provides an environment for beating of a carpet of motile microstructures called cilia as well as cell surface lubrication.⁴⁰ This layer has a watery characteristic due to ciliary forces and a high shear rate in this region.

(b) Sticky, high-viscosity, gel-like and viscoelastic mucus layer⁴¹: it captures foreign particles such as dust, fungi, bacteria, and other environmental particles. This layer has a depth of 5–15 μm in nasal cavity,^42,43 10–30 μm in trachea, and 2–5 μm in bronchi.^18,44,45 Mucus has a high viscosity and relaxation time at low shear rates,⁴⁶ which leads to movement above PCL resembling solid state.

About 2 × 10¹² of motile cilia,^47,48 each with a length of 5–7 μm,⁴⁹ beat in a coordinated wave-like motion within PCL, and propel inhaled foreign particles trapped in mucosal layer toward the trachea and eventually outside of respiratory tract. This process is known as mucociliary clearance (MCC).⁵⁰

Figure 3 shows a schematic diagram of ASL showing ciliary motion, mucus, and PCL layer, as well as the fate of soluble and insoluble particles after deposition in ASL, which rely heavily on characteristics of noncellular pulmonary barrier; this barrier is composed mostly of an ASL blanket in upper and central segments or bronchi, or a pulmonary surfactant layer in distal branches and alveoli, the latter is devoid of mucus in a healthy subject.⁵¹ The peripheral lung has a large surface area (100 m²) and very thin epithelial layer (0.2–0.7 μm), resulting in excellent sink conditions and a desirable environment for gas exchange and high blood perfusion.^52,53

FIG. 3.

Schematic illustration of ASL and MCC and fate of an aerosol after deposition in ASL blanket in upper and central segments or bronchi (left panel) or pulmonary surfactant layer in distal branches and alveoli (right panel). (1) First contact of deposited particles with ASL. (2) Penetration of soluble or ultrafine nanoparticles into pulmonary epithelium and blood. (3) Clearance of large and undissolved particle. Reproduced with permission from Ruge et al.⁵¹ ASL, airway surface liquid; MCC, mucociliary clearance.

Mathematical modeling indicates that efficiency of the MCC is affected by many parameters such as ciliary abnormalities, which can be related to ciliary beat frequency,^54–61 ciliary density,^{55,59,60,62,63} ciliary length,^{58,59,61,64,65} ciliary beat pattern,^55,59,65,66 ciliary beat orientation,^67–69 missing ciliary regions,^55,59 phase difference between cilia,^55,61 cilia lattice geometry,^55,61 and metachronal waves produced by cilia.⁷⁰ MCC can also be affected by ASL and airway characteristics such as mucus rheology,^46,71–77 mucus and PCL viscosity,^56,64,74,78 depth of PCL and mucosal layers,^{56–58,60,61,64,79} mucus production rate,⁷⁹ stiffness transition between PCL and mucosal layers,^56,65,78 high expiratory airflow rates above mucus,^80–87 mucosal dehydration,^88–90 or even airway pressure gradient.⁹¹ Other physiological factors (e.g., sex, posture, age, sleep, exercise etc.) can also influence MCC by changing properties of mucus, PCL, and even cilia, individually or in combination.⁹² The available voluminous literature on biological fluid characteristics in scores of mucus is reviewed thoroughly by Spagnolie.⁹³

Mathematical modeling of inhaled particle deposition on ASL influenced by MCC is critical to evaluate proper drug delivery on ASL to treat common respiratory diseases. Numerical modeling has been regarded as a key tool to achieve this goal. Our review mainly focuses on recent advances in numerical simulation of effect of various parameters on MCC in health and disease followed by mathematical modeling of particle deposition on ASL, which are influenced by MCC.

Numerical Modeling of MCC

Numerical simulation of MCC has been studied extensively⁹⁴ and various aspects of pulmonary ASL modeling were reported by Levy et al.⁹⁵ However, due to the unsteady nature of this problem, small grid and time spacing, a partially large value of ciliary consecutive beat cycle, as well as high computational cost of simulation of ciliary forces and mucus as a viscoelastic fluid, almost all previous studies simulated geometry of MCC in bronchi, which have minimal ASL thickness. From a numerical point of view, there are three major challenges in modeling such complicated geometry: (1) interaction mechanism between cilia and ASL, (2) considering mucus viscoelasticity by employing a proper constitutive equation and, (3) employing an accurate numerical method to solve governing and constitutive equations. Each component will be reviewed in detail.

Ciliary modeling

Mucus clearance is facilitated by unidirectional transport of mucus over an underlying bed of beating cilia. By studying ciliary beat cycle in rabbit trachea,⁹⁶ Fulford and Blake⁵⁰ represented position of each cilium ( $ξ$ ) in x and y directions at 13 time steps (t) by truncated Fourier series as follows:

where s measures length of arc from base of cilium, σ is the angular beat frequency, and a_n and b_n are the Fourier coefficients, which were previously reported by Fulford and Blake.⁵⁰ They also represented location of an array of cilia by⁵⁰ $ξ^{*} (s, t) = (n a i + m b k) + ξ (s, τ), m, n = 0, \pm 1, \pm 2, \dots$ (2)

where a and b are the spacing of cilia in the x and z directions, respectively, and $τ = k n a + σ t$ . Figure 4a shows two-dimensional (2D) ciliary motion in the x-y plane from the cilium beating cycle, which is represented by Equation (1). Figure 4b also shows configuration of the 13 cilia in x and z directions in a symmetric segment of ASL, which has one cycle phase difference in both x and z directions [Eq. (2)]. As Figure 4a illustrates, each cilium performs two main repetitive beating motion called recovery and effective stroke.⁹⁷

FIG. 4.

(a) Diagram of ciliary beat cycle derived by Fulford and Blake.⁵⁰ (b) Schematic illustration of the cilia motion of a symmetric segment of ASL with 13 cilia in x and z directions with one phase difference in both the x and z directions [Eq. (2)]. Note: PCL (y < 4 μm), mucus layer (4 μm < y < 10 μm). PCL, periciliary liquid.

In effective stroke phase (states 1–4 in Fig. 4a), cilia impose the maximum force during the forward motion, which leads to forward flow of ASL. On the other hand, in the recovery stroke phase (states 5–13 in Fig. 4a), cilia slowly move backward by bending down to the epithelial surface, which minimizes their retarding impact on the ASL flow.⁹⁸ The work performed during the effective stroke is magnitudes higher than amount of work done during the recovery stroke.⁹⁸

Due to the complex nature of ciliary motion, as well as the significant effect of ciliary forces on MCC, modeling of ciliary motion and forces has been challenging in the study of this complicated geometry. Barton and Raynor were seemingly the first to model ciliary forces in MCC problem.⁹⁹ They assumed cilium as a rigid rod that is long in effective stroke and becomes short in recovery stroke. They used simplified resistance/drag coefficients to evaluate the effect of ciliary tips on fluid to obtain mucus velocity profile.⁹⁹ Although they calculated realistic mucus flow rates, erroneous depiction of ciliary motion, ignoring effect of PCL and epithelial surface on velocity profile, and coordination wave on an array of cilia are the most important shortcomings of this study.¹⁰⁰

Blake introduced an “active porous medium” model to study ciliary propulsive force by considering cilia as a porous media in the fluid and reported a simplified formulation for velocity profile.¹⁰¹ Other mathematical approaches to study ciliary forces include slender body theory,¹⁰² traction layer model,^74,103,104 continuous envelope moving boundary model,¹⁰⁵ and hybrid grid-particle method.¹⁰⁶ A detailed review of modeling techniques of cilia-ASL interactions, including continuum ciliary modeling and discrete ciliary modeling with its relative merits and demerits, has been recently reviewed thoroughly by Vanaki et al.⁹⁴

Immersed boundary method (IBM) has been regarded as one of the most efficient methods to model stationary and moving bodies in fluid. Peskin appears to be the first to propose IBM and to solve blood flow around heart valves by using this method.¹⁰⁷ In this technique, fluid is represented on the Eulerian domain and moving or stationary boundary immersed in the fluid is represented on Lagrangian points. Force on the immersed object [ $F (X, t)$ ] is calculated based on the difference between real velocity of immersed object and interpolated Eulerian fluid velocity at the Lagrangian node. Then this force should be interpolated on the Eulerian points (fluid domain) using a Dirac delta function as follows¹⁰⁸: $f (x, t) = \int_{Γ} F (X, t) δ_{d} (x - X) d V,$ (3)

where x corresponds the Eulerian coordinates, X is the Lagrangian coordinate system, $δ_{d} (x - X)$ is a Dirac delta function, and d V is the volume of an element on the immersed object. In this equation, force density of the fluid [ $f (x, t)$ ] is computed from the immersed boundary force density [ $F (X, t)$ ] at the Lagrangian domain ( $Γ$ ). A very accurate interpolation operations based on a regularized Dirac delta function is reported by Roma et al. as follows¹⁰⁹: $δ_{d} (x - X) = d (x_{i j k} - X_{l}) \times d (y_{i j k} - Y_{l}) \times d (z_{i j k} - Z_{l})$ (4a) $d (x_{i j k} - X_{l}) = \frac{1}{h} φ (\frac{x_{i j k} - X_{l}}{h})$ (4b)

φ (r) = \{\begin{matrix} \frac{1}{6} (5 - 3 |r| - \sqrt{- 3 {(1 - |r|)}^{2} + 1}), & 0 . 5 \leq |r| \leq 1.5 \\ \frac{1}{3} (1 + \sqrt{- 3 r^{2} + 1}), & |r| \leq 0.5 \\ 0, & o t h e r w i s e, \end{matrix}

(4c)

where $h = Δ x = Δ y = Δ z$ is the spacing between the grids in x, y, and z directions.

Due to efficacy and simplicity of implementation of IBM as well as its robustness when dealing with moving and flexible boundaries, especially in the case of infinitely thin structures, this method has been regarded as a popular technique to couple interaction between cilia and fluid. Thereby, many recent numerical studies have used this method to simulate ciliary forces.^{54–61,65,75,110–115} Because of high computational cost of ciliary modeling, most previous investigations used one-way interaction of cilia on ASL and the effect of ASL flow on ciliary beat pattern was not taken into calculations. However, Loiseau et al. suggested a 2D model that predicts a phase diagram of mucus transport, which influences ciliary beating direction.⁶⁹ Their results showed that cilia/mucus hydrodynamic interactions control collective dynamics of ciliary beat directions.

Modeling of mucus viscoelasticity

Mucus is perceived as a viscoelastic and nonhomogenous fluid.¹¹⁶ Although viscoelastic properties of mucus play a vital role in its flow,^46,74,75,110 due to numerical challenges, many earlier investigators considered mucus as a Newtonian,^{60,61,65,111,115} or non-Newtonian (shear-thinning)^{64,78,106,117} fluid. There are limited studies on viscoelastic behavior of mucus. Smith et al. are seemingly the first to report a numerical study of viscoelastic behavior of mucus by considering it as a linear Maxwell viscoelastic fluid.⁷⁴

Linear viscoelastic models could be useful to simulate viscoelastic fluids in small deformations (<2%–3%)¹¹⁸ and their results are not usually correct for large deformations and fluid flow, and these models are widely used in solid mechanics and polymer engineering for problems in which the size of deformation is very small.¹¹⁹ Due to not satisfying the principle of frame invariance, equations of linear viscoelasticity are not valid for large deformations. To remedy this issue, a set of frame-invariant coordinates are presented to define time derivatives in frames that deform with the material elements.¹¹⁹ Thereby, some recent studies by using this coordinate simulated mucus as quasilinear viscoelastic models, such as upper-convected Maxwell⁸⁷ or Oldroyd-B model,^{54,56,75,110,114} to study the effect of mucus viscoelasticity.

Although these semilinear models such as the Oldroyd-B model can capture the fundamental viscoelastic property of mucus at large deformation, several drawbacks such as assumption of constant viscosity and first normal stress differences coefficient, as well as inability to model second normal stress differences of fluid in shear¹¹⁹ remain challenging in simulating this fluid. Besides, mucus typically has nonlinear rheology properties at greater strains.⁴⁶ Therefore, the nonlinear characteristics and behaviors of these biological systems as a function of shear rate are important factors in numerical simulation.¹²⁰ To study the effect of nonlinear behavior of mucus, Vasquez et al. reported a comprehensive equation of mucus.⁴⁶ They used a five-mode nonlinear Giesekus model in approximating the experimentally verified relaxation spectrum of mucus. In this model, the elastic contribution of stress tensor ( $τ$ ) is calculated by adding elastic stress tensor of each mode as follows:

Each mode of the Giesekus model can be calculated as follows¹²¹:

where, $λ_{a}$ , $η_{a}^{p}$ , and $α_{a}$ are relaxation time, elastic viscosity contribution, and mobility parameter of each mode, respectively. Total viscosity of mucus ( $η_{M}$ ) is obtained by adding solvent ( $η_{M}^{s}$ ) and elastic viscosity ( $η_{M, a}^{p}$ ) contributions:

In addition, in Equation (6), D is rate of deformation tensor and is upper convected derivative defined as follows:

D = \frac{1}{2} (\nabla u + \nabla u^{T})

(8)

Nonlinear mobility term in the Giesekus model [Eq. (6)] is a key parameter to capture the shear-thinning behavior of mucus, as well as second normal stress difference,¹¹⁹ while these behaviors are ignored in some semilinear models such as Oldroyd-B model. Due to ability of the Giesekus model to capture shear thinning viscosity and first and second normal stress difference in shear,¹²¹ this model can be a desirable model to simulate mucus viscoelasticity. Table 1 presents values of parameters based on the five-mode Giesekus model of mucus in healthy state.⁴⁶ Our group recently used this model to replicate 2D⁵⁸ and 3D^55,57 simulation of MCC geometry.

Table 1.

The Five-Mode Giesekus Model Parameters of Mucus in Healthy State⁴⁶

Mode	Relaxation time, $λ_{a}$ (s)	Elastic part of mucus viscosity, $η_{M, a}^{p}$ (Pa.s)	Giesekus parameter, $α_{a}$
1	0.0089	0.0298	0.2
2	0.0821	0.062	0.3
3	0.466	0.2493	0.5
4	3.1290	1.4215	0.5
5	49.733	42.2034	0.5

Numerical simulation of governing equations

Recent investigations have elucidated that PCL contains a fine matrix network that prevents mucus from penetrating it.⁴⁰ Thereby, in most previous numerical investigations, ASL is modeled as two-phase immiscible fluids with different properties.⁷⁴ For this purpose, a proper boundary condition should be applied at the interface of the PCL and mucus layers. The modified Navier–Stokes equations, including immersed boundary force term (to simulate ciliary and PCL-mucus membrane forces) as well as elastic stress tensor (to estimate viscoelastic behavior of the mucus), have been presented as governing equations of this geometry as follows⁶⁰: $\nabla \cdot u = 0,$ (10a) $\frac{\partial u}{\partial t} + u . \nabla u = - \frac{1}{ρ} \nabla p + υ_{s} \nabla^{2} u + \frac{1}{ρ} \nabla . τ + f,$ (10b)

where $ρ$ is density of ASL, $υ_{s}$ is Newtonian kinematic viscosity, p is static pressure, t is time, u is ASL velocity vector, and $τ$ is elastic stress tensor (due to mucus viscoelasticity). Equations (10a) and (10b) correspond to continuity and momentum equations, respectively. Momentum equation includes forces f acting on fluid due to boundary force of cilia and force on PCL-mucus interface.⁶⁰ Ciliary forces can be calculated by IBM explained in Ciliary Modeling section and elastic force exerted by membrane interface of PCL and mucus layer can be calculated as follows⁶⁰: $F_{M e m} = \frac{1}{ρ} \frac{\partial}{\partial s} [T (s) κ (s)]$ (11)

where $T (s)$ and $κ (s)$ are defined as follows:

T (s) = T_{0} (\frac{\partial X_{M e m}}{\partial s_{0}} - 1)

(12a)

In Equations (12a) and (12b), is a unit tangential vector to PCL-mucus membrane interface. Arc lengths s and s₀ are measured along current and initial configuration of membrane interface, respectively. Scalar T₀ describes stiffness constant, which shows elastic characteristics of flexible boundary. A few studies in 2D numerical simulations considered effect of surface tension at PCL-mucus interface by computing boundary membrane force explained in Equation (11).^56,60,65,75 However, previous experimental and numerical investigations have shown that interface force between mucus and PCL does not have a significant effect on mucus and PCL velocities^{56,60,74,96,122}; therefore, most numerical studies ignored the effect of membrane interface between PCL and mucus layers by assuming a flat interface between these two layers and imposing a continuous fluid velocity and stress across interface in their numerical modeling.

In most previous studies, a symmetric segment of ASL (Fig. 4b) with following boundary conditions has been investigated:

ASL flow is considered periodic in x and z directions.⁷⁴

Bottom boundary is assumed to be a no-slip wall.⁶⁰

Top boundary is perceived to be a slip wall.¹²³

This two-phase model has the advantage of simulating mucus and PCL layers as a unique domain with different properties simultaneously. It can also estimate penetration of cilia into moving mucus layer reported in previous investigations.^{50,96,124–127}

Table 2 represents a summary of 2D and 3D numerical investigations of MCC along with a detailed numerical method, ciliary modeling method, and mucus rheology model. As this table indicates, due to the complex nature of this geometry, bulk of literature is limited to 2D flow regimen. However, 2D MCC overestimates mean mucus velocity since fluid is not allowed to flow around cilia.⁵⁷

Table 2.

Summary of Studies Available on the Two-Dimensional and Three-Dimensional Numerical Simulations of Mucociliary Clearance

Nature of study	Investigator	Cilia modeling method	Numerical method	Mucus model
Two dimensional	Smith et al.^74,169	Fourier expansion	FDM	Viscoelastic: Maxwell
	Dillon et al.¹²²	IBM	FDM	Viscoelastic: Jeffreys
	Lukens et al.¹⁷⁰	IBM	FDM	Newtonian
	Lee et al.⁶⁰ and Jayathilake et al.⁶⁵	IBM	FDM	Newtonian
	Roy et al.¹⁷¹	Fourier expansion	FEM	Viscoelastic: Maxwell
	Kurbatova et al.¹⁰⁴	Porous medium	FDM	Newtonian
	Sedaghat et al.^54,56,75,114	IBM	LBM	Viscoelastic: Oldroyd-B
	Guo and Kanso¹¹⁰	IBM	FVM	Viscoelastic: Oldroyd-B
	George and Pidaparti⁹¹	Elastic model	FEM	Newtonian
	Sedaghat et al.⁵⁸	IBM	FDM	Viscoelastic: 5-mode Giesekus
	Modaresi and Shirani⁸⁷	Velocity profile	FVM	Viscoelastic: Upper-Convected Maxwell
	Modaresi¹¹⁷	Velocity profile	FVM	Non-Newtonian (Power Law and Thixotropic)
Three dimensional	Mitran¹²⁴	Fixed moving Grid	FVM	Viscoelastic: Maxwell
	Jayathilake et al.⁶¹ and Vanaki et al.⁵⁹	IBM	FDM	—
	Chatelin and Poncet^64,106 and Chatelin et al.¹⁷²	Penalization method	FDM	Non-Newtonian (shear-thinning viscosity)
	Chateau et al.^111,112	IBM	LBM	Newtonian
	Quek et al.¹¹⁵	IBM	FDM	Newtonian
	Sedaghat et al.^55,57	IBM	FDM	Viscoelastic: 5-mode Giesekus

The numerical method include FDM, FEM, FVM, and LBM.

FDM, finite difference method; FEM, finite element method; FVM, finite volume method; IBM, immersed boundary method; LBM, lattice Boltzmann method.

Due to the unsteady nature of this problem, a partially large value of ciliary consecutive beat cycle (about 0.1 seconds in healthy state), small grid, time spacing of geometry, and high computational cost of numerical methods to simulate ciliary forces and mucus viscoelasticity, the 3D simulation of this complex geometry, even in Newtonian mucosal fluid, has not been extensively researched. Jayathilake et al., in 3D numerical simulation, studied the effect of various parameters on PCL flow.⁶¹ They used IBM to simulate ciliary forces and finite difference projection method to simulate PCL as a Newtonian fluid.

Their results showed that maximum PCL velocity happens if cilia have phase differences in both stream-wise and span-wise directions. They also showed that PCL velocity is linearly dependent on ciliary beat frequency, while average PCL velocity decreases with increasing PCL height. Furthermore, shortened cilia significantly reduce PCL velocity, while changing PCL viscosity does not substantially affect PCL velocity.

Vanaki et al. did a study with the same approach by studying some ciliary abnormalities found in various chronic respiratory diseases, for example, ciliary beat frequency, ciliary length, ciliary beat pattern, ciliary density, and missing ciliary region on the PCL transport.⁵⁹ Their results showed that PCL surface velocity and flow rate are significantly affected by changes in ciliary beat frequency and ciliary length observed in COPD patients. They also showed that stiff planar and wiper ciliary beat patterns in primary ciliary dyskinesia patients significantly reduce PCL surface velocity and flow rate.

Furthermore, these ciliary abnormalities decrease PCL velocity and flow rate compared to a healthy state. Although these two studies^59,61 have been considered two of the most comprehensive 3D numerical simulations of PCL geometry, they did not consider mucus layer as an integral part of this geometry in their mathematical modeling. Our group has recently studied 3D simulation of MCC geometry by considering PCL as a Newtonian and mucus as a nonlinear Giesekus viscoelastic fluid.⁵⁷

In our study, the projection finite difference method was used to solve governing and constitutive equations of fluid flow. In addition, to simulate ciliary forces more accurately, IBM was employed. Ciliary configuration at t = 0 seconds as well as computational domain of this study are shown in Figure 5. As this figure indicates, cilia have cyclic motion according to truncated Fourier series, as reported by Fulford and Blake⁵⁰ [Eq. (2) and Fig. 4a]. Values for standard parameter set (parameters used in a healthy state) of mucus rheology properties and ASL characteristics employed in this 3D simulation are represented in Tables 1 and 3, respectively. Some numerical simplifications have been used in this modeling, which include the following:

FIG. 5.

Schematic geometry of 3D MCC in study of Sedaghat et al.⁵⁷ Reproduced with permission.

Table 3.

Specifications of Airway Surface Liquid Under Healthy Conditions in Numerical Studies of Sedaghat Et Al.^{54–58,75,114}

Parameter	Standard value	Explored range
Cilia length	6 μm	5–7⁴⁹
Cilia beat frequency	60 rad/s	50–90¹⁷³
Cilia diameter	0.2 μm	0.2–0.3¹⁷⁴
Cilia spacing (in x and z directions)	3 μm	3⁶⁰
PCL deph	6 μm	5–10³⁹
Mucus depth (bronchi)	4 μm	2–5^30,44,45
Mucus and PCL density	1000 kg/m³	1000⁷⁴
PCL viscosity	0.001 Pa.s	0.001⁷⁴
Newtonian part of mucus viscosity	0.001 Pa.s	0.001⁷⁴

PCL, periciliary liquid.

(1)

Air-mucus and PCL-mucus interfaces are considered to be flat.⁷⁴ Furthermore, this force is neglected due to minimal membrane force between PCL and mucus layers.^{56,60,74,96,122}

(2)

The cilia move in evenly spaced filament-like structures with a constant beat frequency.^59,61,112

(3)

One-way interaction of cilia on ASL is considered, which means propulsive impact of cilia on ASL is taken into consideration; however, cilia are made to move in a predefined beat pattern and influence of ASL flow on cilia was ignored.⁶⁹

(4)

Due to limited height of PCL and mucus layer, the effect of ASL gravity force would be negligible compared to the ciliary forces.¹²⁸

(5)

No mass transfer is present at the bottom epithelial cell wall.

Three-dimensional simulation in a healthy state predicts mean mucus and PCL velocities of 44.33 and 20.49 μm/s, respectively, which agree with experimentally determined clearance rates.^129,130 In contrast, our 2D simulation⁵⁸ with the same approach and rheology properties for mucus showed ∼2.5-fold larger clearance rates with respect to 3D models. This is mainly due to neglecting ciliary spacing in the z-direction in 2D simulation, which overestimates mucus and PCL velocities since fluid is not allowed to flow around cilia in the z-direction. 3D numerical results⁵⁷ also showed that for PCL layer, even though flow direction mainly changes along the z direction, PCL has a 3D flow structure because of localized ciliary forces.

On the other hand, in the mucus layer, effects of ciliary motions have been mainly dampened and velocity change in z direction is much less than that in PCL due to significantly higher mucus viscosity and elasticity. Our 3D results⁵⁷ showed that mean mucus and PCL velocities enhance linearly with ciliary beating frequency. However, the effect of this parameter on mucus velocity is more than PCL velocity. Furthermore, with increasing ciliary length, from 5 to 6.5 μm, mean mucus velocity almost doubled. In contrast, PCL velocity was much less affected, resulting from cilia penetration into the mucus layer. Numerical results also indicated that changes in PCL depth and mucus depth have little effect on mucus and PCL velocities.

In another 3D numerical investigation with the same approach, our group has recently studied the effect of cilia abnormalities such as phase difference between cilia, cilia beat pattern, cilia beat frequency, cilia lattice geometry, missing cilia regions, and cilia density on MCC.⁵⁵ Results indicated that some cilia abnormalities, such as cilia density, cilia beat pattern, and cilia beat frequency, have a dominant effect on MCC and 3D structure of the fluid, and some abnormalities, such as missing cilia regions and phase differences between cilia, have a moderate influence on that.

Numerical Simulation of Particle Deposition on ASL

Simulation of particle transport and deposition in various parts of human airways have been extensively studied and bulk of the literature has been reviewed recently.^131,132 However, there needs to be more numerical activity on the fate of particles deposited on ASL by considering the effect of MCC. The protective nature of mucus against foreign particles and efficient inhaled drug delivery makes it challenging to tackle this issue from a computational standpoint.

Experiments showed that proteins and viruses <100 nm in diameter, such as h-lysozyme (3.5 nm), horse myoglobin (3.8 nm), chicken lysozyme (4.1 nm), porcine pepsin (4.5 nm), h IgM Fab (6 nm), horse ferritin (12 nm), and viruses, including polio (28 nm), Norwalk (38 nm), hepatitis B (43 nm), human papillomavirus (55 nm), adenoviruses (60–90 nm), and rotavirus (75 nm), could easily penetrate mucus.^11,12 Influenza virus (100 nm) can also penetrate mucus unless strongly immobilized by adhesive interactions with mucosal components.¹⁴

The concentration distribution of soluble particles or trajectories of insoluble particles in ASL are affected by ASL velocity profile, as well as various forces. Gravity, drag, lift, and Brownian forces are regarded as important forces that influence particle motion. There are two main methods for studying particle deposition in airways: Eulerian and Lagrangian approaches. Eulerian method, which is based on solving classical mass diffusion equation, is appropriate for the simulation of soluble or very small nanoparticles where the Brownian effect of particle is dominant.¹³³ On the other hand, Lagrangian method, which is based on solving particle equation of motion, is more efficient for simulation of insoluble larger particles where particle inertia effect is more important.¹³⁴ In the following sections, numerical studies dealing with simulation of particle transport on ASL influenced by MCC using these two approaches will be reviewed.

Eulerian modeling

Nanoparticle transport based on Eulerian modeling is the most used model to simulate solutes or ultrafine particles; the motion and deposition are dominated by their Brownian forces. In this condition, Lagrangian trajectory analysis becomes inefficient since the particle relaxation times are very small and a large number of particles need to be traced for precise estimation of local particle deposition.¹³⁵ In addition, simplicity, low computational cost, and rapid evaluation of particle deposition rates are some other characteristics of Eulerian simulation compared to Lagrangian modeling.

To simulate particle deposition on ASL based on Eulerian method, by assuming dissolution of particles in mucus and using velocity distribution of the ASL, classical mass diffusion equation at ASL needs to be solved and the concentration of inhaled particles at various parts of ASL versus time is obtained. Equation of mass diffusion of an inhaled particle can be presented as follows: $\frac{\partial c}{\partial t} + c . \nabla u = D_{A S L} \nabla^{2} c$ (13)

where c is ASL particle concentration (mass of particle/volume of ASL) and $D_{A S L}$ is diffusion coefficient of solutes through ASL (mucus and PCL layers). As stated earlier, mucus is a viscoelastic fluid interwoven by a network of flexible mucin fibers. To diffuse through this medium, particles must pass through pores of these fibers. In addition, the movement of solutes through this 3D network is hindered by these fibers, leading to decreased space for diffusion.^136–142

Olmsted et al. reported diffusion coefficient of various proteins and viruses with different sizes in human cervical mucus.¹⁴³ Yang et al. also calculated diffusion coefficient of the swine influenza virus in porcine mucus.¹⁴⁴ Although diffusion characteristics are affected by particles' size and shape, globular or chain-like solutes can adopt themselves in a 3D configuration in fluid medium.¹⁴⁵ Desai et al. reported effective diffusion coefficients of various solutes through aqueous and porcine mucus layers.¹⁴⁶ One of the preferred mathematical models of diffusion coefficients of solute in mucus layer was reported as the Obstruction-Scaling model¹³⁸: $\frac{D_{M}}{D_{0}} = e^{\frac{- π}{4} {(\frac{r_{p} + r_{f}}{r_{g} + r_{f}})}^{2}}$ (14a) $D_{0} = \frac{K_{B} T}{6 π η_{0} r_{p}}$ (14b)

In Equation (14a), $r_{f} = 3.5 n m$ is mucin fiber radius, $r_{g} = 50 n m$ is mucin network's effective mesh fiber spacing,¹⁴⁷ andr_p is the effective radius of solute agent. In addition, D_M is mucus diffusion coefficient and D₀ is corresponding solute diffusion coefficient in water, estimated by Stokes-Einstein equation as Equation (14b). In this equation, $η_{0}$ is water viscosity, K_B is the Boltzmann's constant, and T is the temperature. In addition, Ariane et al. reported a wide range of drug diffusion coefficients in PCL by considering the effects of drug diffusion, ciliary beat frequency, and ciliary flexibility.¹⁴⁸

In the Eulerian approach, mass conservation boundary condition at the interface separating layers should be imposed to obtain the concentration of particles in the multilayer domain. For example, between PCL and mucus layer, this boundary condition is given by the following equation: $D_{M} \frac{\partial c_{M}}{\partial n} = D_{P C L} \frac{\partial c_{P C L}}{\partial n}$ (15)

The same equation can simulate mass conservation between air-mucus¹⁴⁹ or PCL-epithelium^149,150 layers.

There have been a few studies dealing with simulation of particle deposition in ASL solved by Eulerian method. In some of these studies, an attempt has been made to simulate particle transport in both airway and ASL; however, due to numerical complexities arising from this geometry, there are some simplifications in the simulation of ASL flow. Corley et al. studied inhaled acrolein concentration in various parts of rat, monkey, and human respiratory tract by considering steady-state inhalation.¹⁵¹ Their results showed that regional acrolein uptake was influenced by airway geometry, airflow rates, inhaled acrolein concentrations, ASL diffusion coefficient, ASL thickness, and maximum rate of metabolism.

In another study, they extended their previous findings by considering transient breathing patterns for predicting kinetics of three reactive constituents of cigarette smoke: acrolein, acetaldehyde, and formaldehyde.¹⁵² Yoo and Ito used an integrated model of actual shape of human body geometry with a virtual airway of a realistic human respiratory tract to simulate time-dependent inhaled formaldehyde concentration, including adsorption distributions, that is, heterogeneous tissue dosimetry in respiratory tract (Fig. 6).¹⁴⁹ As this figure indicates, they used Eulerian approach to study formaldehyde adsorption/deposition flux or equilibrium concentration in a three-layer model of airway, mucus, and subepithelium. Their results showed that over 50% of inhaled formaldehyde was adsorbed by ASL at nasal cavity.

FIG. 6.

A schematic geometry of numerical model in study of Yoo and Ito,¹⁴⁹ which was done by PBPK-computational fluid dynamics (CFD) hybrid model. Reproduced with permission. CFD, computational fluid dynamics; PBPK, physiologically based pharmacokinetic.

Rygg and Longest¹⁵³ and Rygg et al.^154,155 studied dissolution, absorption, and clearance of drug aerosols in nasal mucus. In their studies, a flat 2D surface of inner nasal wall is considered to study particle deposition in a 10-micron mucus layer. The study also investigated dissolution of three different drugs on nasal mucus layer, as well as penetration of particles into epithelial surface over time. Shang et al. extended these works by considering mucus as a 3D surface-shell model.¹⁵⁶ Considering ASL as a single fluid with an average viscosity, as well as ignoring drug adsorption by epithelium are important limitations of this study. In addition, in studies by Rygg and Longest,¹⁵³ Rygg et al.,^154,155 and Shang et al.,¹⁵⁶ a uniform mass flow rate is defined to impose a fixed axial velocity in the ASL layer.

Recently Chari et al. simulated a realistic 3D nasal cavity airway model of nose of a healthy nonsmoker person.¹⁵⁰ They employed a cutting-edge 3D meshing technology that enables smooth capture of mucus layer at a micrometer scale and comparatively large air-flow domain (Fig. 7). To obtain a more precise velocity profile of ASL results from beating of cilia, a wall-driven velocity assumption at interface between mucus and PCL layer was employed. They studied influence of partition coefficient, particle size, drug solubility, as well as particle distribution on mucus and its uptake in epithelium. They showed that drugs with a higher partition coefficient or solubility can be taken up easier in the epithelium. At the same dose, smaller particles with their relatively large surface areas dissolve faster and are absorbed easier in ASL compared to larger ones.

FIG. 7.

A schematic geometry of nasal geometry and mucosal layer in study of Chari et al.¹⁵⁰ Mucosal layer (black) separates airflow of nasal surface (dark gray) from epithelial layer (light gray). Reproduced with permission.

Ariane et al. used a discrete multiphysics method to simulate diffusive and convective drug delivery in PCL.¹⁴⁸ In this study, ASL and cilia were simulated by smoothed particle hydrodynamics (SPH) and mass-spring models, respectively, in a 2D simulation. They studied effects of ciliary flexibility, ciliary beat frequency, and drug diffusion on PCL drug delivery. Results showed that at a low-frequency regime, cilia hinder drug delivery. In contrast, diffusion is dominant drug delivery mechanism at an intermediate-frequency regime. Cilia play a vital role in mass transfer in ASL at a high-frequency regime due to convection mechanism.

It should be noted that in some of these studies, for simplicity, particles were assumed to dissolve immediately into mucus after depositing on upper surface of ASL layer.^149,157 However, some investigations studied particle size reduction owing to dissolution and later absorption in mucus layer^{148,150,153–155} according to Noyes–Whitney mass diffusion equation: $\frac{d m}{d t} = \frac{A D_{M} (C_{s} - C_{b})}{h}$ (16)

where m is mass of the drug, A is surface area of the particle, D_m is diffusion coefficient of the dissolved drug in the mucus layer, C_s is solubility of the drug, C_b is concentration of the dissolved drug in bulk phase, and h is diffusion layer thickness.

The above numerical studies to study particle dissolution in ASL assume either constant mucus velocity or some simple assumption for ciliary forces. Accurate simulation of governing equation in ASL, to obtain more precise velocity distribution in PCL and mucus layers, which leads to calculation of convective term in diffusion equation, is scarce in literature.

Lagrangian modeling

In this method, individual particles are tracked within ASL flow field. The primary advantage of this approach is that the effect of effective forces such as gravity, drag, inertia, and Brownian forces are directly considered in the motion of the particle. In addition, accurate velocity and location of particles, as well as percentage of particles that deposit on ciliary or epithelial walls can also be estimated by this method.

Second law of Newton governs dynamics of particle motion and for each particle, it can be represented as follows¹³⁴: $\frac{d u_{p}}{d t} = F_{D} + F_{L} + F_{B} + \frac{g (ρ_{p} - ρ_{f})}{ρ_{p}}$ (17a) $\frac{d x_{p}}{d t} = u_{p}$ (17b)

where u _p is particle velocity vector, F _D is drag force per unit mass, F _L is lift force per unit mass, F _B is Brownian force per unit mass, g is acceleration gravity vector, and x _p is particle position vector.

Drag force can be calculated as follows¹⁵⁸: $F_{D} = \frac{3 μ_{a} C_{D} {Re}_{p}}{4 ρ_{p} d_{p}^{2} C_{c}} (u - u_{p})$ (18a) ${Re}_{p} = \frac{ρ_{f} |u - u_{p}| d_{p}}{μ_{f}}$ (18b)

C_{D} = \frac{24}{{Re}_{p}} (1 + 0.15 {Re}_{p}^{0.687})

(18c)

C_{c} = 1 + \frac{2 λ}{d_{p}} [1.257 + 0.4 exp (- 1.1 \frac{d_{p}}{2 λ})]

(18d)

where d_p is particle diameter, $ρ_{p}$ and $ρ_{f}$ are particle and fluid density, $μ_{f}$ is fluid dynamic viscosity, and ${Re}_{p}$ and C_Dare particle Reynolds number and drag coefficient, respectively, which can be calculated as Equations (18b) and (18c). In addition, C_c is Cunningham correction factor, which is given by Equation (18d). In this equation, $λ$ is mean free path along particle trajectories.

Particle mass diffusivity in ASL can be obtained as follows¹⁵⁹: $D_{A S L} = \frac{K T C_{c}}{3 π μ_{f} d_{p}}$ (19)

where K is the Boltzmann constant and T is fluid temperature. Based on the above equation, drag force for a small spherical particle in ASL can be calculated as follows: $F_{D} = \frac{3 K T}{4 ρ_{p} d_{p}^{3} D_{A S L}} (u - u_{p})$ (20)

Diffusion coefficient of a particle in the ASL ( $D_{A S L}$ ) can be estimated by Obstruction-Scaling model [Eq. (14)].⁹²

Shear lift force can be calculated from a generalized expression of Saffman lift force, shown as follows¹⁶⁰: $F_{L} = \frac{2 k ν^{\frac{1}{2}} ρ_{f} D}{ρ_{p} d_{p} {(D D)}^{\frac{1}{4}}} (u - u_{p})$ (21)

where $ν$ is fluid kinematic viscosity, k = 2.594 is constant coefficient of Saffman lift force,¹⁶¹ and D is deformation rate tensor, which is calculated by Equation (8).

Brownian force per unit mass of particles can also be calculated as follows¹⁶²: $F_{B} = ξ \sqrt{\frac{π S_{0}}{Δ t}}$ (22a) $S_{0} = \frac{216 ν K T}{π^{2} ρ_{f} {d_{p}}^{5} {(\frac{ρ_{p}}{ρ_{f}})}^{2} C_{c}}$ (22b)

In Equation (22a), $ξ$ is the zero-mean, unit-variance-independent Gaussian random numbers, and $Δ t$ is the integration time step.

There have been a few studies dealing with simulation of nanoparticle deposition in ASL by Lagrangian trajectory analysis approach. Kirch et al. simulated deposition of deposited particles in ASL and clearance of these particles in a 2D domain.^163,164 Mucus flow velocity was calculated by solving 2D steady-state equations of motion of fluid in Eulerian domain analytically.

By solving time-dependent Lagrangian equation of motion of nanoparticles, velocity components and location of each particle in ASL were estimated. Their results showed that ASL velocity distribution plays a vital role in the fate of nanoparticles in ASL, while sedimentation, diffusion, and impaction seem to have little effect on particle transport. Also, intrinsic plasticity of mucus slabs and collision of such slabs increase particle translocation toward epithelium. Two-dimensional simulation, ignoring the effect of Brownian forces of nanoparticles, and using an inaccurate model to simulate ciliary forces and ASL velocity distribution are main shortcomings of these studies.^163,164

Ernst et al., by considering mucus as a porous structure of Newtonian fluid-filled random-sized cells, calculated subdiffusion of particles in mucus layer (human sputum taken from CF patients), by measuring mean squared displacements of the particles.¹⁶⁵ They considered the effect of Brownian diffusion of particles in permeable membrane domain. The foundation of this strategy is modeling of particle trajectories in the presence of Robin boundaries.¹⁶⁶ They showed that only particles with a diameter of <40 nm can pass through a mucosal layer by passive Brownian motion. This finding contradicts results of previous studies that imply 100 nm and smaller particles can easily diffuse through mucosal layer through Brownian motion regardless of mucus mesh spacing.³⁵

Bajd and Serša¹⁶⁷ modified the method presented by Ernst et al.¹⁶⁵ and used bond-fluctuation model to simulate motion of different-size and chain-like particles through mucus as a porous media with three different 2D restrictive-scaffold geometries (isotropic, anisotropic, and random). Their results showed that steric interactions play a vital role in the motion of chain particles in confining mucus fluid. The main shortcoming of the above two studies is not taking into effect the ASL velocity on particles' motion and velocity.^165,167

Vanaki et al. simulated deposition of nanoparticles in 3D ciliary arrays interacting with PCL.¹⁶⁸ They used IBM to estimate ciliary forces and location versus time and also projection method for solving a 3D Navier–Stokes equation in ASL and a particle tracking method to simulate transport of inhaled particles deposited onto the PCL in healthy and diseased lungs. Although this study was the first that solved complete 3D Navier–Stokes equations in ASL to address velocity and location of particles more accurately, it had two main drawbacks. First, it did not consider mucus layer in simulation, while mucus has been perceived as one of the most important protective barriers against inhaling particles. Second, it used small nanoparticles in simulations and effect of particle Brownian motion was not studied, while Brownian force is dominant force in this situation.

An essential step in simulating more accurate circumstances for both healthy and diseased states is developing a 3D model that can depict interaction of MCC and mucus-deposited inhaled particles. While micron-sized particles seem to be cleared by the mucus effectively, specific nano-sized particles have been shown to penetrate ASL into and beyond PCL.²⁵ These particles have both therapeutic (retention of hazardous particles in lungs for an extended period of time) and toxicological consequences (leveraging of mucus-penetrating nanoparticles for drug delivery to bronchial epithelium), which can only be simulated with a 3D cilia configuration and ASL model. However, 3D simulation studies of mucus and PCL layers by considering effect of ciliary motion and simulation of particle motion by considering impact of effective forces [explained in Eq. (17)–(22)] are indeed scarce.

Conclusion and Outlook

Numerical simulation of MCC as well as particle deposition in human airways have been studied extensively. Current state-of-the-art is only limited to study particle deposition in human airway combined with ASL, which is influenced by MCC. Findings shed further light into drug deposition as well as penetration of harmful particles into human respiratory system.

Studies reveal that three key parameters have influenced numerical simulation of MCC, including ciliary forces, mucus rheology properties, and numerical method, to solve time-dependent modified Navier–Stokes equations. The effect of various parameters on both mucosal and PCL velocities has been reported in previous numerical investigations by overcoming these three issues. In this regard, IBM , one of the most effective cilia modeling techniques many researchers use, is extensively reviewed. The five-mode nonlinear Giesekus model in approximating experimentally verified relaxation spectrum (reported by Vasquez et al.⁴⁶) has been explained as one of complete constitutive equations of mucus used in some recent numerical studies.

Protective nature of mucus against large particles as well as efficient inhaled drug delivery have been controversial research issues. A small percentage of current literature addresses particle deposition in human airways influenced by MCC. Two main approaches have been used to study particle transport in ASL: Eulerian and Lagrangian. Eulerian approach, which is based on solving classical mass diffusion equation, is appropriate for the simulation of dissolved nanoparticles (<100 nm) in mucosal layer. It has some advantages compared to Lagrangian method, which includes simplicity, low computational cost, and rapid evaluation of particle deposition rates.

Most studies have used Eulerian method to simulate particle deposition on ASL. Lagrangian method is based on solving particle equation of motion by considering effect of inertia, drag, lift, gravity, and Brownian forces imposed on each particle. This method is more efficient for simulation of larger particles where the effect of particle inertia is more significant. Tracing exact velocity and location of each particle as well as exact location of deposited particles (i.e., cilia, epithelium wall, etc.) are advantages of this method. However, there have been few studies dealing with particle deposition on ASL using this method.

Reliable results based on Eulerian and Lagrangian approaches for studying distribution of dissolved drug powders or drug solutions in ASL are available up to using fixed value of ASL velocity. The effect of real ASL flow regime on drug distribution in this layer is still not well studied. Nonhomogeneous and viscoelastic behavior of mucus layer as well as consideration of mucus shear-thinning viscosity may also influence particle deposition rate. Considering an accurate velocity, distribution of ASL can lead to a significant effect on particle distribution in this layer. Most of this literature is based on numerical simplifications of governing differential equations in this area and 3D numerical results are scarce.

Footnotes

Authors' Contributions

M.H.S.: conceptualization, data curation, formal analysis, investigation, methodology, project administration, writing original draft, and writing review and editing. M.B.: project administration, writing review and editing, and investigation. O.A.: data curation, investigation, methodology, project administration, supervision, and writing review and editing. All three authors approved the final version.

Author Confirmation Statement

First author (M.H.S.) is from Technical and Vocational University (Tehran, Iran) and the third author (O.A.) is from Shiraz University (Shiraz, Iran); both institutions are primarily involved in education or research.

Author Disclosure Statement

The authors declare they have no conflicting financial interest.

Funding Information

No funding was received for this article.

Reviewed by:

Georg Dietze

Tim Corcoran

References

Livraghi-Butrico

, Grubb

, Wilkinson

, et al. Contribution of mucus concentration and secreted mucins Muc5ac and Muc5b to the pathogenesis of muco-obstructive lung disease. Mucosal Immunol, 2017; 10(2):395–407; doi:10.1038/mi.2016.63

Roy

, Livraghi-Butrico

, Fletcher

, et al. Muc5b is required for airway defence. Nature, 2014; 505(7483):412–416; doi:10.1038/nature12807

Ensign

, Cone

, Hanes

. Oral drug delivery with polymeric nanoparticles: The gastrointestinal mucus barriers. Adv Drug Delivery Rev, 2012; 64(6):557–570; doi: 10.1016/j.addr.2011.12.009

Lundquist

, Artursson

. Oral absorption of peptides and nanoparticles across the human intestine: Opportunities, limitations and studies in human tissues. Adv Drug Delivery Rev, 2016; 106:256–276; doi: 10.1016/j.addr.2016.07.007

Curlin

, Bursac

. Cervical mucus: From biochemical structure to clinical implications. Front Biosci (Schol Ed), 2013; 5(2):507–515; doi: 10.2741/S386

Hein

, Valore

, Helmig

, et al. Antimicrobial factors in the cervical mucus plug. Am J Obstet Gynecol, 2002; 187(1):137–144; doi: 10.1067/mob.2002.123034

Gipson

, Spurr-Michaud

, Tisdale

. Human conjunctival goblet cells express the membrane associated mucin MUC16: Localization to mucin granules. Exp Eye Res, 2016; 145:230–234; doi: 10.1016/j.exer.2015.12.009

Bansil

, Turner

. Mucin structure, aggregation, physiological functions and biomedical applications. Curr Opin Colloid Interface Sci, 2006; 11(2–3):164–170; doi: 10.1016/j.cocis.2005.11.001

Lai

, Wang

Y-Y

, Wirtz

, et al. Micro-and macrorheology of mucus. Adv Drug Deliv Rev, 2009; 61(2):86–100; doi: 10.1016/j.addr.2008.09.012

10.

Carlstedt

, Sheehan

. Structure and macromolecular properties of cervical mucus glycoproteins. Sym Soc Exp Biol, 1989; 43:289–316.

11.

Thornton

, Sheehan

. From mucins to mucus: Toward a more coherent understanding of this essential barrier. Proc Am Thorac Soc, 2004; 1(1):54–61; doi: 10.1513/pats.2306016

12.

Wickstrom

, Davies

, Eriksen

, et al. MUC5B is a major gel-forming, oligomeric mucin from human salivary gland, respiratory tract and endocervix: Identification of glycoforms and C-terminal cleavage. Biochem J, 1998; 334(3):685–693; doi: 10.1042/bj3340685

13.

Murgia

, Loretz

, Hartwig

, et al. The role of mucus on drug transport and its potential to affect therapeutic outcomes. Adv Drug Deliv Rev, 2018; 124:82–97; doi: 10.1016/j.addr.2017.10.009

14.

Schuster

, Suk

, Woodworth

, et al. Nanoparticle diffusion in respiratory mucus from humans without lung disease. Biomaterials, 2013; 34(13):3439–3446; doi: 10.1016/j.biomaterials.2013.01.064

15.

Adivitiya, Kaushik

, Chakraborty

, et al. Mucociliary respiratory epithelium integrity in molecular defense and susceptibility to pulmonary viral infections. Biology, 2021; 10(2):95; doi: 10.3390/biology10020095

16.

Leal

, Smyth

, Ghosh

. Physicochemical properties of mucus and their impact on transmucosal drug delivery. Int J Pharm, 2017; 532(1):555–572; doi: 10.1016/j.ijpharm.2017.09.018

17.

Lieleg

, Ribbeck

. Biological hydrogels as selective diffusion barriers. Trends Cell Biol, 2011; 21(9):543–551; doi: 10.1016/j.tcb.2011.06.002

18.

Sanders

, Rudolph

, Braeckmans

, et al. Extracellular barriers in respiratory gene therapy. Adv Drug Deliv Rev, 2009; 61(2):115–127; doi: 10.1016/j.addr.2008.09.011

19.

Cone

RA.

Barrier properties of mucus. Adv Drug Deliv Rev, 2009; 61(2):75–85; doi: 10.1016/j.addr.2008.09.008

20.

Suk

, Lai

, Boylan

, et al. Rapid transport of muco-inert nanoparticles in cystic fibrosis sputum treated with N-acetyl cysteine. Nanomedicine, 2011; 6(2):365–375; doi: 10.2217/nnm.10.123

21.

Edwards

, Ben-Jebria

, Langer

. Recent advances in pulmonary drug delivery using large, porous inhaled particles. J Appl Physiol, 1998; 85(2):379–385; doi: 10.1152/jappl.1998.85.2.379

22.

Yang

, Peters

, Williams III

. Inhaled nanoparticles—A current review. Int J Pharm, 2008; 356(1–2):239–247; doi: 10.1016/j.ijpharm.2008.02.011

23.

Azarmi

, Roa

, Löbenberg

. Targeted delivery of nanoparticles for the treatment of lung diseases. Adv Drug Deliv Rev, 2008; 60(8):863–875; doi: 10.1016/j.addr.2007.11.006

24.

Huckaby

, Lai

. PEGylation for enhancing nanoparticle diffusion in mucus. Adv Drug Deliv Rev, 2018; 124:125–139; doi: 10.1016/j.addr.2017.08.010

25.

Lai

, Wang

Y-Y

, Hanes

. Mucus-penetrating nanoparticles for drug and gene delivery to mucosal tissues. Adv Drug Deliv Rev, 2009; 61(2):158–171; doi: 10.1016/j.addr.2008.11.002

26.

Pearson

, Chater

, Wilcox

. The properties of the mucus barrier, a unique gel—How can nanoparticles cross it?. Ther Deliv, 2016; 7(4):229–244; doi: 10.4155/tde-2015-0002

27.

Schneider

, Xu

, Boylan

, et al. Nanoparticles that do not adhere to mucus provide uniform and long-lasting drug delivery to airways following inhalation. Sci Adv, 2017; 3(4):e1601556; doi: 10.1126/sciadv.1601556

28.

Sigurdsson

, Kirch

, Lehr

C-M

. Mucus as a barrier to lipophilic drugs. Int J Pharm, 2013; 453(1):56–64; doi: 10.1016/j.ijpharm.2013.05.040

29.

Suk

, Lai

, Wang

Y-Y

, et al. The penetration of fresh undiluted sputum expectorated by cystic fibrosis patients by non-adhesive polymer nanoparticles. Biomaterials, 2009; 30(13):2591–2597; doi: 10.1016/j.biomaterials.2008.12.076

30.

Sanders

, De Smedt

, Van Rompaey

, et al. Cystic fibrosis sputum: A barrier to the transport of nanospheres. Am J Respir Crit Care Med, 2000; 162(5):1905–1911; doi: 10.1164/ajrccm.162.5.9909009

31.

Duncan

, Kim

, Colon-Cortes

, et al. An adeno-associated viral vector capable of penetrating the mucus barrier to inhaled gene therapy. Mol Ther Methods Clin Dev, 2018; 9:296–304; doi: 10.1016/j.omtm.2018.03.006

32.

Lai

, Suk

, Pace

, et al. Drug carrier nanoparticles that penetrate human chronic rhinosinusitis mucus. Biomaterials, 2011; 32(26):6285–6290; doi: 10.1016/j.biomaterials.2011.05.008

33.

Chisholm

, Shenoy

, Shade

, et al. Nanoparticle diffusion in spontaneously expectorated sputum as a biophysical tool to probe disease severity in COPD. Eur Respir J, 2019; 54(2):1900088; doi: 10.1183/13993003.00088-2019

34.

Huang

, Chisholm

, Zhuang

, et al. Protein nanocages that penetrate airway mucus and tumor tissue. Proc Natl Acad Sci U S A, 2017; 114(32):E6595–E6602; doi: 10.1073/pnas.1705407114

35.

Newby

, Seim

, Lysy

, et al. Technological strategies to estimate and control diffusive passage times through the mucus barrier in mucosal drug delivery. Adv Drug Deliv Rev, 2018; 124:64–81; doi: 10.1016/j.addr.2017.12.002

36.

Crystal

, Randell

, Engelhardt

, et al. Airway epithelial cells: Current concepts and challenges. Proc Am Thorac Soc, 2008; 5(7):772–777; doi: 10.1513/pats.200805-041HR

37.

Rokicki

, Rokicki

, Wojtacha

, et al. The role and importance of club cells (Clara cells) in the pathogenesis of some respiratory diseases. Kardiochir Torakochirurgia Pol, 2016; 13(1):26–30; doi: 10.5114/kitp.2016.58961

38.

Mercer

, Russell

, Roggli

, et al. Cell number and distribution in human and rat airways. Am J Respir Cell Mol Biol, 1994; 10(6):613–624; doi: 10.1165/ajrcmb.10.6.8003339

39.

Pedersen

, Holstein-Rathlou

N-H

, Larsen

, et al. Fluid absorption related to ion transport in human airway epithelial spheroids. Am J Physiol Lung Cell Mol Physiol, 1999; 277(6):L1096–L1103; doi: 10.1152/ajplung.1999.277.6.L1096

40.

Button

, Cai

L-H

, Ehre

, et al. A periciliary brush promotes the lung health by separating the mucus layer from airway epithelia. Science, 2012; 337(6097):937–941; doi: 10.1126/science.1223012

41.

Rubin

BK.

Physiology of airway mucus clearance. Respir Care, 2002; 47(7):761–768.

42.

Beule

AG.

Physiology and pathophysiology of respiratory mucosa of the nose and the paranasal sinuses. GMS Curr Top Otorhinolaryngol Head Neck Surg, 2010; 9:Doc07; doi: 10.3205/cto000071

43.

Ugwoke

, Agu

, Verbeke

, et al. Nasal mucoadhesive drug delivery: background, applications, trends and future perspectives. Adv Drug Deliv Rev, 2005; 57(11):1640–1665; doi: 10.1016/j.addr.2005.07.009

44.

Patton

JS.

Mechanisms of macromolecule absorption by the lungs. Adv Drug Deliv Rev, 1996; 19(1):3–36; doi: 10.1016/0169-409X(95)00113-L

45.

Wine

JJ.

The genesis of cystic fibrosis lung disease. J Clin Invest, 1999; 103(3):309–312; doi: 10.1172/JCI6222

46.

Vasquez

, Jin

, Palmer

, et al. Modeling and simulation of mucus flow in human bronchial epithelial cell cultures–Part I: Idealized axisymmetric swirling flow. PLoS Comput Biol, 2016; 12(8):e1004872; doi: 10.1371/journal.pcbi.1004872

47.

Wanner

, Salathé

, O'Riordan

. Mucociliary clearance in the airways. Am J Respir Crit Care Med, 1996; 154(6):1868–1902; doi: 10.1164/ajrccm.154.6.8970383

48.

Bustamante-Marin

, Ostrowski

. Cilia and mucociliary clearance. Cold Spring Harb Perspect Biol, 2017; 9(4):a028241; doi: 10.1101/cshperspect.a028241

49.

King

, Gilboa

, Meyer

, et al. On the transport of mucus and its rheologic simulants in ciliated systems. Am Rev Respir Dis, 1974; 110(6P1):740–745.

50.

Fulford

, Blake

. Muco-ciliary transport in the lung. J Theor Biol, 1986; 121(4):381–402; doi: 10.1016/S0022-5193(86)80098-4

51.

Ruge

, Kirch

, Lehr

C-M

. Pulmonary drug delivery: From generating aerosols to overcoming biological barriers—Therapeutic possibilities and technological challenges. Lancet Respir Med, 2013; 1(5):402–413; doi: 10.1016/S2213-2600(13)70072-9

52.

Labiris

, Dolovich

. Pulmonary drug delivery. Part I: physiological factors affecting therapeutic effectiveness of aerosolized medications. Br J Clin Pharmacol, 2003; 56(6):588–599; doi: 10.1046/j.1365-2125.2003.01892.x

53.

Patton

, Fishburn

, Weers

. The lungs as a portal of entry for systemic drug delivery. Proc Am Thorac Soc, 2004; 1(4):338–344; doi: 10.1513/pats.200409-049TA

54.

Sedaghat

, Shahmardan

, Norouzi

, et al. Effect of cilia beat frequency on muco-ciliary clearance. J Biomed Phys Eng, 2016; 6(4):265.

55.

Sedaghat

, Sadrizadeh

, Abouali

. Three-dimensional simulation of mucociliary clearance under the ciliary abnormalities. J Nonnewton Fluid Mech, 2023; 316:15029; doi: 10.1016/j.jnnfm.2023.105029

56.

Sedaghat

, Shahmardan

, Norouzi

, et al. Numerical simulation of muco-ciliary clearance: immersed boundary-lattice Boltzmann method. Comput Fluids, 2016; 131:91–101; doi: j.compfluid.2016.03.015

57.

Sedaghat

, Farnoud

, Schmid

, et al. Nonlinear simulation of mucociliary clearance: A three-dimensional study. J Nonnewton Fluid Mech, 2022; 300:104727; doi: 10.1016/j.jnnfm.2021.104727

58.

Sedaghat

, George

, Abouali

. A nonlinear viscoelastic model of mucociliary clearance. Rheol Acta, 2021; 60(6–7):371–384; doi: 10.1007/s00397-021-01274-4

59.

Vanaki

, Holmes

, Jayathilake

, et al. Three-dimensional numerical analysis of periciliary liquid layer: Ciliary abnormalities in respiratory diseases. Appl Sci, 2019; 9(19):4033; doi: 10.3390/app9194033

60.

Lee

, Jayathilake

, Tan

, et al. Muco-ciliary transport: Effect of mucus viscosity, cilia beat frequency and cilia density. Comput Fluids, 2011; 49(1):214–221; doi: 10.1016/j.compfluid.2011.05.016

61.

Jayathilake

, Tan

, Le

, et al. Three-dimensional numerical simulations of human pulmonary cilia in the periciliary liquid layer by the immersed boundary method. Comput Fluids, 2012; 67:130–137; doi: 10.1016/j.compfluid.2012.07.016

62.

Pellicciotta

, Das

, Kotar

, et al. Cilia density and flow velocity affect alignment of motile cilia from brain cells. J Exp Biol, 2020; 223(24):jeb.229310; doi:10.1242/jeb.229310

63.

Pellicciotta

, Hamilton

, Kotar

, et al. Entrainment of mammalian motile cilia in the brain with hydrodynamic forces. Proc Natl Acad Sci U S A, 2020; 117(15):8315–8325; doi: 10.1073/pnas.1910065117

64.

Chatelin

, Poncet

. A parametric study of mucociliary transport by numerical simulations of 3D non-homogeneous mucus. J Biomech, 2016; 49(9):1772–1780; doi: 10.1016/j.jbiomech.2016.04.009

65.

Jayathilake

, Le

, Tan

, et al. A numerical study of muco-ciliary transport under the condition of diseased cilia. Comput Methods Biomech Biomed Eng, 2015; 18(9):944–951; doi: 10.1080/10255842.2013.864285

66.

Viallat

, Loiseau

, Gsell

, et al. Long-range self-organization of ciliary activity and flow patterns in reconstituted bronchial epithelium. In: 72nd Annual Meeting of the APS Division of Fluid Dynamics, American Physical Society, 2019, Seattle, WA, pp. Q31–004.

67.

Gsell

, Loiseau

, D'ortona

, et al. Hydrodynamic model of directional ciliary-beat organization in human airways. Sci Rep, 2020; 10(1):1–12; doi: 10.1038/s41598-020-64695-w

68.

Gsell

, Loiseau

, D'Ortona

, et al. Flow-driven orientational self-organization of the bronchial epithelium. In: APS March 2021 Meeting, American Physical Society, 2021, USA, pp. R09–002.

69.

Loiseau

, Gsell

, Nommick

, et al. Active mucus–cilia hydrodynamic coupling drives self-organization of human bronchial epithelium. Nat Phys, 2020; 16(11):1158–1164; doi: 10.1038/s41567-020-0980-z

70.

Mesdjian

, Wang

, Gsell

, et al. Longitudinal to transverse metachronal wave transitions in an in vitro model of ciliated bronchial epithelium. Phys Rev Lett, 2022; 129(3):038101; doi: 10.1103/PhysRevLett.129.038101

71.

Lafforgue

, Seyssiecq

, Poncet

, et al. Rheological properties of synthetic mucus for airway clearance. J Biomed Mater Res A, 2018; 106(2):386–396; doi: 10.1002/jbm.a.36251

72.

Lafforgue

, Seyssiecq

, Poncet

, et al. A time-dependent Herschel-Bulkley model for synthetic bronchial mucus. In: Conference: Proceedings of the Canadian Society for Mechanical Engineering International (CSME 2020). 2020.

73.

Lafforgue

, Poncet

, Seyssiecq

, et al. Rheological characterization of macromolecular colloidal gels as simulant of bronchial mucus. AIP Conference Proceedings, Dec 14, 2017, Vol. 14, No. 1, p. 110003. ALP Publishing LLC.

74.

Smith

, Gaffney

, Blake

. A viscoelastic traction layer model of muco-ciliary transport. Bull Math Biol, 2007; 69(1):289–327; doi: 10.1007/s11538-006-9177-6

75.

Sedaghat

, Shahmardan

, Norouzi

, et al. On the effect of mucus rheology on the muco-ciliary transport. Math Biosci, 2016; 272:44–53; doi: 10.1016/j.mbs.2015.11.010

76.

Galko

, Gsell

, D'Ortona

, et al. Numerical study of oscillating non-Newtonian flows, application to a pulmonary airway clearance device. In: 73rd Annual Meeting of the APS Division of Fluid Dynamics, American Physical Society, 2020; Nov. 22–24. USA, pp. H01–001.

77.

Wang

, Gsell

, D'Ortona

, et al. Dynamics of a cilium/cilia beating in 3D non-Newtonian flow. 72nd Annual Meeting of the APS Division of Fluid Dynamics, American Physical Society, 2019; Nov. 23–26. Seattle, WA, pp. Q31–002.

78.

Chatelin

, Anne-Archard

, Murris-Espin

, et al. Numerical and experimental investigation of mucociliary clearance breakdown in cystic fibrosis. J Biomech, 2017; 53:56–63; doi: 10.1016/j.jbiomech.2016.12.026

79.

Kim

, Greene

, Sankaran

, et al. Mucus transport in the airways by two-phase gas-liquid flow mechanism: Continuous flow model. J Appl Physiol, 1986; 60(3):908–917; doi: 10.1152/jappl.1986.60.3.908

80.

Kim

, Rodriguez

, Eldridge

, et al. Criteria for mucus transport in the airways by two-phase gas-liquid flow mechanism. J Appl Physiol, 1986; 60(3):901–907; doi: 10.1152/jappl.1986.60.3.901

81.

Rajendran

, Banerjee

. Mucus transport and distribution by steady expiration in an idealized airway geometry. Med Eng Phys, 2019; 66:26–39; doi: 10.1016/j.medengphy.2019.02.006

82.

Paz

, Suárez

, Parga

, et al. Glottis effects on the cough clearance process simulated with a CFD dynamic mesh and Eulerian wall film model. Comput Methods Biomech Biomed Engin, 2017; 20(12):1326–1338; doi: 10.1080/10255842.2017.1360872

83.

Paz

, Suárez

, Vence

, et al. Analysis of the volume of fluid (VOF) method for the simulation of the mucus clearance process with CFD. Comput Methods Biomech Biomed Eng, 2019; 22(5):547–566; doi: 10.1080/10255842.2019.1569637

84.

Kim

, Iglesias

, Sackner

. Mucus clearance by two-phase gas-liquid flow mechanism: asymmetric periodic flow model. J Appl Physiol, 1987; 62(3):959–971; doi: 10.1152/jappl.1987.62.3.959

85.

King

, Brock

, Lundell

. Clearance of mucus by simulated cough. J Appl Physiol, 1985; 58(6):1776–1782; doi: 10.1152/jappl.1985.58.6.1776

86.

, Wang

, Feng

. Computational analysis of obstructive disease and cough intensity effects on the mucus transport and clearance in an idealized upper airway model using the volume of fluid method. Phys Fluids, 2021; 33(2):021903; doi: 10.1063/5.0037764

87.

Modaresi

, Shirani

. Effects of continuous and discrete boundary conditions on the movement of upper-convected maxwell and Newtonian mucus layers in coughing and sneezing. Eur Phys J Plus, 2022; 137(7):846; doi: 10.1140/epjp/s13360-022-03067-x

88.

Boucher

RC.

Airway surface dehydration in cystic fibrosis: Pathogenesis and therapy. Annu Rev Med, 2007; 58:157–170; doi: 10.1146/annurev.med.58.071905.105316

89.

Kim

, Baumlin

, Yoshida

, et al. Losartan rescues inflammation-related mucociliary dysfunction in relevant models of cystic fibrosis. Am J Respir Crit Care Med, 2020; 201(3):313–324; doi: 10.1164/rccm.201905-0990OC

90.

Button

, Boucher

, Group

UoNCVL

. Role of mechanical stress in regulating airway surface hydration and mucus clearance rates. Respir Physiol Neurobiol, 2008; 163(1–3):189–201; doi: 10.1016/j.resp.2008.04.020

91.

George

, Pidaparti

. Airway pressure gradient may decrease the beating amplitude of cilia. Front Physiol, 2019; 7:157; doi: 10.3389/fphy.2019.00157

92.

Houtmeyers

, Gosselink

, Gayan-Ramirez

, et al. Regulation of mucociliary clearance in health and disease. Eur Respir J, 1999; 13(5):1177–1188; doi: 10.1034/j.1399-3003.1999.13e39.x

93.

Spagnolie

SE.

Complex fluids in biological systems. In: Biological and Medical Physics, Biomedical Engineering. Springer Science and Business Media: New York; 2015.

94.

Vanaki

, Holmes

, Saha

, et al. Muco-ciliary clearance: A review of modelling techniques. J Biomech, 2019; 109578; doi: 10.1016/j.jbiomech.2019.109578

95.

Levy

, Hill

, Forest

, et al. Pulmonary Fluid Flow Challenges for Experimental and Mathematical Modeling. The Society for Integrative and Comparative Biology; 2014.

96.

Sanderson

, Sleigh

. Ciliary activity of cultured rabbit tracheal epithelium: Beat pattern and metachrony. J Cell Sci, 1981; 47(1):331–347; doi: 10.1242/jcs.47.1.331

97.

Sleigh

. Cilia and Flagella. Academic Press: 1974.

98.

Gueron

, Levit-Gurevich

. Energetic considerations of ciliary beating and the advantage of metachronal coordination. Proc Natl Acad Sci U S A, 1999; 96(22):12240–12245; doi: 10.1073/pnas.96.22.12240

99.

Barton

, Raynor

. Analytical investigation of cilia induced mucous flow. Bull Math Biophys, 1967; 29(3):419–428; doi: 10.1007/BF02476581

100.

Blake

Mucus flows. Math Biosci, 1973; 17(3):301–313; doi: 10.1016/0025-5564(73)90073-4

101.

Blake

An active porous medium model for ciliary propulsion. J Theor Biol, 1977; 64(4):697–701; doi: 10.1016/0022-5193(77)90268-5

102.

Hancock

The self-propulsion of microscopic organisms through liquids. Proc Roy Soc London Ser A, 1953; 217(1128):96–121; doi: 10.1098/rspa.1953.0048

103.

Blake

, Winet

. On the mechanics of muco-ciliary transport. Biorheology, 1980; 17(1–2):125–134; doi: 10.1093/imamat/32.1-3.69

104.

Kurbatova

, Bessonov

, Volpert

, et al. Model of mucociliary clearance in cystic fibrosis lungs. J Theor Biol, 2015; 372:81–88; doi: 10.1016/j.jtbi.2015.02.023

105.

Bottier

, Peña Fernández

, Pelle

, et al. A new index for characterizing micro-bead motion in a flow induced by ciliary beating: Part II, modeling. PLoS Comput Biol, 2017; 13(7):e1005552; doi: 10.1371/journal.pcbi.1005552

106.

Chatelin

, Poncet

. A hybrid grid-particle method for moving bodies in 3D stokes flow with variable viscosity. SIAM J Sci Comput, 2013; 35(4):B925–B949; doi: 10.1137/120892921

107.

Peskin

CS.

Flow patterns around heart valves: A numerical method. J Comput Phys, 1972; 10(2):252–271; doi: 10.1016/0021-9991(72)90065-4

108.

Breugem

W-P.

A second-order accurate immersed boundary method for fully resolved simulations of particle-laden flows. J Comput Phys, 2012; 231(13):4469–4498; doi: 10.1016/j.jcp.2012.02.026

109.

Roma

, Peskin

, Berger

. An adaptive version of the immersed boundary method. J Comput Phys, 1999; 153(2):509–534; doi: 10.1006/jcph.1999.6293

110.

Guo

, Kanso

. A computational study of mucociliary transport in healthy and diseased environments. Eur J Comput Mech, 2017; 26(1–2):4–30; doi: 10.1080/17797179.2017.1321206

111.

Chateau

, d'Ortona

, Poncet

, et al. Transport and mixing induced by beating cilia in human airways. Front Physiol, 2018; 9:161; doi: 10.3389/fphys.2018.00161

112.

Chateau

, Favier

, D'ortona

, et al. Transport efficiency of metachronal waves in 3D cilium arrays immersed in a two-phase flow. J Fluid Mech, 2017; 824:931–961; doi: 10.1017/jfm.2017.352

113.

Chateau

, Favier

, Poncet

, et al. Why antiplectic metachronal cilia waves are optimal to transport bronchial mucus. Phys Rev E, 2019; 100(4):042405; doi: 10.1103/PhysRevE.100.042405

114.

Sedaghat

, Bagheri

, Shahmardan

, et al. A hybrid immersed boundary-Lattice Boltzmann method for simulation of viscoelastic fluid flows interaction with complex boundaries. Commun Comput Phys, 2021; 29(5):1411–1445; doi: 10.4208/cicp.OA-2019-0158

115.

Quek

YLR

, Lim

, Chiam

K-H

. Three-dimensional computational model of multiphase flow driven by a bed of active cilia. Comput Fluids, 2018; 170:222–235; doi: 10.1016/j.compfluid.2018.05.001

116.

Puchelle

, Zahm

, Quemada

. Rheological properties controlling mucociliary frequency and respiratory mucus transport. Biorheology, 1987; 24(6):557–563; doi: 10.3233/BIR-1987-24606

117.

Modaresi

Numerical investigation of mucociliary clearance using power law and thixotropic mucus layers under discrete and continuous cilia motion. Biomech Model Mechanobiol, 2023; 22:253–269; doi: 10.1007/s10237-022-01645-7

118.

Funk

, Hall

, Crandall

, et al. Linear and quasi-linear viscoelastic characterization of ankle ligaments. J Biomech Eng, 2000; 122(1):15–22; doi: 10.1115/1.429623

119.

Bird

, Armstrong

, Hassager

Dynamics of polymeric liquids. Vol. 1: Fluid Mechanics. Wiley: New York; 1987.

120.

Lauga

Propulsion in a viscoelastic fluid. Phys Fluids, 2007; 19(8):083104; doi: 10.1063/1.2751388

121.

Giesekus

A simple constitutive equation for polymer fluids based on the concept of deformation-dependent tensorial mobility. J Nonnewton Fluid Mech, 1982; 11(1–2):69–109; doi: 10.1016/0377-0257(82)85016-7

122.

Dillon

, Fauci

, Omoto

, et al. Fluid dynamic models of flagellar and ciliary beating. Ann N Y Acad Sci, 2007; 1101(1):494–505; doi: 10.1196/annals.1389.016

123.

Dauptain

, Favier

, Bottaro

. Hydrodynamics of ciliary propulsion. J Fluids Struct, 2008; 24(8):1156–1165; doi: 10.1016/j.jfluidstructs.2008.06.007

124.

Mitran

SM.

Metachronal wave formation in a model of pulmonary cilia. Comput Struct, 2007; 85(11):763–774; doi: 10.1016/j.compstruc.2007.01.015

125.

Sleigh

MA.

The nature and action of respiratory tract cilia. In: Respiratory Defense Mechanisms. ( Brain

, Proctor

, Reid

. eds.) Marcel Dekker: New York; 1977; pp. 247–288.

126.

Liron

, Rozenson

. Muco-ciliary transport. J Submicrosc Cytol Pathol, 1983; 15(1):317–321

127.

Puchelle

, Herard

A-L

, Zahm

J-M

. Airway mucociliary epithelium injury and repair. In: Cilia, Mucus, and Mucociliary Interactions. ( Baum

, Priel

, Roth

, et al. eds.) Dekker: New York; 1998; pp. 203–217.

128.

Blake

On the movement of mucus in the lung. J Biomech, 1975; 8(3–4):179–190.

129.

Matsui

, Randell

, Peretti

, et al. Coordinated clearance of periciliary liquid and mucus from airway surfaces. J Clin Invest, 1998; 102(6):1125; doi: 10.1172/JCI2687

130.

Salathe

, O'Riordan

, Wanner

Mucociliary clearance. In: The Lung: Scientific Foundations. ( Crystal

, West

, Barnes

, et al. eds.) Lippencott-Raven: Philadelphia, PA; 1997; pp. 2295–2308.

131.

Longest

, Bass

, Dutta

, et al. Use of computational fluid dynamics deposition modeling in respiratory drug delivery. Expert Opin Drug Deliv, 2019; 16(1):7–26; doi: 10.1080/17425247.2019.1551875

132.

Islam

, Paul

, Ong

, et al. A review of respiratory anatomical development, air flow characterization and particle deposition. Int J Environ Res Public Health, 2020; 17(2):380; doi: 10.3390/ijerph17020380

133.

Michaelides

, Crowe

, Schwarzkopf

(eds). Biological systems and biomimetics. In: Multiphase Flow Handbook (2nd Ed.). CRC Press: 2016; pp. 887–954; doi: 10.1201/9781315371924

134.

Rastegar

, Ahmadi

, Babu

. Filtration of aqueous colloidal ceria slurries using fibrous filters—An experimental and simulation study. Sep Purif Technol, 2017; 176:231–242; doi: 10.1016/j.seppur.2016.12.017

135.

, Ahmadi

. Computer simulation of particle deposition in the upper tracheobronchial tree. Aerosol Sci Technol, 1995; 23(2):201–223; doi: 10.1080/02786829508965304

136.

Ogston

, Preston

, Wells

. On the transport of compact particles through solutions of chain-polymers. Proc Roy Soc London Ser A, 1973; 333(1594):297–316.

137.

Perrins

, McKenzie

, McPhedran

. Transport properties of regular arrays of cylinders. Proc Roy Soc London Ser A, 1979; 369(1737):207–225.

138.

Amsden

An obstruction-scaling model for diffusion in homogeneous hydrogels. Macromolecules, 1999; 32(3):874–879; doi: 10.1021/ma980922a

139.

Renkin

EM.

Filtration, diffusion, and molecular sieving through porous cellulose membranes. J Gen Physiol, 1954; 38(2):225.

140.

Anderson

, Quinn

. Restricted transport in small pores: a model for steric exclusion and hindered particle motion. Biophys J, 1974; 14(2):130–150.

141.

Johansson

, Elvingson

, Loefroth

. Diffusion and interaction in gels and solutions. 3. Theoretical results on the obstruction effect. Macromolecules, 1991; 24(22):6024–6029; doi: 10.1021/ma00022a019

142.

Johansson

, Löfroth

. Diffusion and interaction in gels and solutions. 4. Hard sphere Brownian dynamics simulations. J Chem Phys, 1993; 98(9):7471–7479; doi: 10.1063/1.464686

143.

Olmsted

, Padgett

, Yudin

, et al. Diffusion of macromolecules and virus-like particles in human cervical mucus. Biophys J, 2001; 81(4):1930–1937; doi: 10.1016/S0006-3495(01)75844-4

144.

Yang

, Steukers

, Forier

, et al. A beneficiary role for neuraminidase in influenza virus penetration through the respiratory mucus. PLoS One, 2014; 9(10):e110026; doi: 10.1371/journal.pone.0110026

145.

, Saltzman

. Mathematical modeling of molecular diffusion through mucus. Adv Drug Deliv Rev, 2009; 61(2):101–114; doi: 10.1016/j.addr.2008.09.006

146.

Desai

, Mutlu

, Vadgama

. A study of macromolecular diffusion through native porcine mucus. Experientia, 1992; 48:22–26; doi: 10.1007/BF01923598

147.

Rutland

, Cole

, Fox

. Ciliary dyskinesia. Lancet, 1980; 316(8199):859; doi: 10.1016/S0140-6736(80)90203-2

148.

Ariane

, Kassinos

, Velaga

, et al. Discrete multi-physics simulations of diffusive and convective mass transfer in boundary layers containing motile cilia in lungs. Comput Biol Med, 2018; 95:34–42; doi: 10.1016/j.compbiomed.2018.01.010

149.

Yoo

S-J

, Ito

. Assessment of transient inhalation exposure using in silico human model integrated with PBPK-CFD hybrid analysis. Sustain Cities Soc, 2018; 40:317–325; doi: 10.1016/j.scs.2018.04.023

150.

Chari

, Sridhar

, Walenga

, et al. Computational analysis of a 3D mucociliary clearance model predicting nasal drug uptake. J Aerosol Sci, 2021; 155:105757; doi: 10.1016/j.jaerosci.2021.105757

151.

Corley

, Kabilan

, Kuprat

, et al. Comparative computational modeling of airflows and vapor dosimetry in the respiratory tracts of rat, monkey, and human. Toxicol Sci, 2012; 128(2):500–516; doi: 10.1093/toxsci/kfs168

152.

Corley

, Kabilan

, Kuprat

, et al. Comparative risks of aldehyde constituents in cigarette smoke using transient computational fluid dynamics/physiologically based pharmacokinetic models of the rat and human respiratory tracts. Toxicol Sci, 2015; 146(1):65–88; doi: 10.1093/toxsci/kfv071

153.

Rygg

, Longest

. Absorption and clearance of pharmaceutical aerosols in the human nose: development of a CFD model. J Aerosol Med Pulm Drug Deliv, 2016; 29(5):416–431; doi: 10.1089/jamp.2015.1252

154.

Rygg

, Hindle

, Longest

. Absorption and clearance of pharmaceutical aerosols in the human nose: effects of nasal spray suspension particle size and properties. Pharm Res, 2016; 33(4):909–921; doi: 10.1007/s11095-015-1837-5

155.

Rygg

, Hindle

, Longest

. Linking suspension nasal spray drug deposition patterns to pharmacokinetic profiles: A proof-of-concept study using computational fluid dynamics. J Pharm Sci, 2016; 105(6):1995–2004; doi: 10.1016/j.xphs.2016.03.033

156.

Shang

, Inthavong

, Tu

. Development of a computational fluid dynamics model for mucociliary clearance in the nasal cavity. J Biomech, 2019; 85:74–83; doi: 10.1016/j.jbiomech.2019.01.015

157.

Shang

, Inthavong

, Qiu

, et al. Prediction of nasal spray drug absorption influenced by mucociliary clearance. PLoS One, 2021; 16(1):e0246007; doi: 10.1371/journal.pone.0246007

158.

Hinds

WC.

Aerosol Technology: Properties, Behavior, and Measurement of Airborne Particles. John Wiley & Sons; 1999.

159.

Ounis

, Ahmadi

, McLaughlin

. Brownian diffusion of submicrometer particles in the viscous sublayer. J Colloid Interface Sci, 1991; 143(1):266–277; doi: 10.1016/0021-9797(91)90458-K

160.

Tian

, Ahmadi

. Particle deposition in turbulent duct flows—Comparisons of different model predictions. J Aerosol Sci, 2007; 38(4):377–397; doi: 10.1016/j.jaerosci.2006.12.003

161.

Saffman

PG.

The lift on a small sphere in a slow shear flow. J Fluid Mech, 1965; 22(2):385–400; doi: 10.1017/S0022112065000824

162.

, Ahmadi

. Dispersion and deposition of spherical particles from point sources in a turbulent channel flow. Aerosol Sci Technol, 1992; 16(4):209–226; doi: 10.1080/02786829208959550

163.

Kirch

, Guenther

, Schaefer

, et al. Computational fluid dynamics of nanoparticle disposition in the airways: Mucus interactions and mucociliary clearance. Comput Vis Sci, 2011; 14(7):301–308; doi: 10.1007/s00791-012-0184-x

164.

Kirch

, Guenther

, Doshi

, et al. Mucociliary clearance of micro-and nanoparticles is independent of size, shape and charge—An ex vivo and in silico approach. J Control Release, 2012; 159(1):128–134; doi: 10.1016/j.jconrel.2011.12.015

165.

Ernst

, John

, Guenther

, et al. A model for the transient subdiffusive behavior of particles in mucus. Biophys J, 2017; 112(1):172–179; doi: 10.1016/j.bpj.2016.11.900

166.

Singer

, Schuss

, Osipov

, et al. Partially reflected diffusion. SIAM J Appl Math, 2008; 68(3):844–868; doi: 10.1137/060663258

167.

Bajd

, Serša

. A Bond-fluctuation model of translational dynamics of chain-like particles through mucosal scaffolds. Biophys J, 2018; 114(11):2732–2742; doi: 10.1016/j.bpj.2018.04.031

168.

Vanaki

, Holmes

, Suara

, et al. Transport and fate of inhaled particles after deposition onto the airway surface liquid: A 3D numerical study. Comput Biol Med, 2020; 117:103595; doi: 10.1016/j.compbiomed.2019.103595

169.

Smith

, Gaffney

, Blake

. A model of tracer transport in airway surface liquid. Bull Math Biol, 2007; 69(3):817–836; doi: 10.1007/s11538-006-9163-z

170.

Lukens

, Yang

, Fauci

. Using Lagrangian coherent structures to analyze fluid mixing by cilia. Chaos, 2010; 20(1):017511; doi: 10.1063/1.3271340

171.

Roy

, Sivathanu

, Das

. Effect of bunching of cilia and their interplay on muco-ciliary transport. Comput Biol Med, 2013; 43(11):1758–1772; doi: 10.1016/j.compbiomed.2013.08.014

172.

Chatelin

, Poncet

, Didier

, et al. Mucus and ciliated cells of human lung: Splitting strategies for particle methods and 3D stokes flows. Procedia IUTAM, 2015; 18:114–122; doi: 10.1016/j.piutam.2015.11.012

173.

Knowles

, Boucher

. Mucus clearance as a primary innate defense mechanism for mammalian airways. J Clin Invest, 2002; 109(5):571–577; doi: 10.1172/JCI0215217

174.

Chang

MM-J

, Shih

, Wu

Pulmonary epithelium: Cell types and functions. In: The pulmonary epithelium in health and disease. John Wiley & Sons, Ltd: Chichester, UK; 2008; pp. 1–26.