Pulmonary drug delivery and retention: A computational study to identify plausible parameters based on a coupled airway-mucus flow model

Pulmonary drug delivery systems rely on inhalation of drug-laden aerosols produced from aerosol generators such as inhalers, nebulizers etc. On deposition, the drug molecules diffuse in the mucus layer and are also subjected to mucociliary advection which transports the drugs away from the initial deposition site. The availability of the drug at a particular region of the lung is, thus, determined by a balance between these two phenomena. A mathematical analysis of drug deposition and retention in the lungs is developed through a coupled mathematical model of aerosol transport in air as well as drug molecule transport in the mucus layer. The mathematical model is solved computationally to identify suitable conditions for the transport of drug-laden aerosols to the deep lungs. This study identifies the conditions conducive for delivering drugs to the deep lungs which is crucial for achieving systemic drug delivery. The effect of different parameters on drug retention is also characterized for various regions of the lungs, which is important in determining the availability of the inhaled drugs at a target location. Our analysis confirms that drug delivery efficacy remains highest for aerosols in the size range of 1-5 μm. Moreover, it is observed that amount of drugs deposited in the deep lung increases by a factor of 2 when the breathing time period is doubled, with respect to normal breathing, suggesting breath control as a means to increase the efficacy of drug delivery to the deep lung. A higher efficacy also reduces the drug load required to be inhaled to produce the same health effects and hence, can help in minimizing the side effects of a drug.


I. IDEALIZATION OF THE LUNG GEOMETRY
summarises the magnitudes of various parameters used while approximating the lung geometry (see Eqs 1-3 in the main manuscript). Table B lists the assumed fraction of airway area (γ) that is alveolated at each generation in the lung model. This is required in calculating aerosol deposition in the alveolar region (see Eqs S32 and S34).

A. Aerosol transport
The one-dimensional transport equation for aerosols at any location in the idealised lung geometry is expressed as where, c a represents the aerosol concentration, Q represents the volume flow rate of air in breathing, and D a represents the diffusivity of aerosols in air. Eq S1 is equivalent to Eq 4 in the main manuscript. The coefficient L D accounts for the aerosols deposited in the airway mucus. This equation is based on the trumpet model proposed by Taulbee & Yu 4 which has been later used by various authors to study different aspects of aerosol deposition in the lung [5][6][7] . The transport equation is formulated based on the assumption that the aerosols are monodispersed, do not undergo coagulation, and are decoupled from airflow in the lungs. It is also assumed that external forces (such as electrical and magnetic forces) do not have any influence on the aerosol dynamics. It is further assumed that there is no additional source of aerosols present within the lungs and the aerosols are either deposited in the airway mucus or washed out of the airways.
Eq S1 is presented in terms of airway length (x), while the lung model adopted (see Fig 1 in the main manuscript) is in terms of lung generation number (N ). As such, Eq S1 needs to be converted to a more appropriate form in terms of N . This requires an additional mathematical relation (Eq S2) connecting airway length x and the lung generation number (N ) given by which can be derived by considering Eq 2 in the main manuscript. Converting Eq S1 using Eqs S2 and A N = A 0 (2β) N , we get where, q(t) represents the temporal sinusoidal function accounting for airflow variation during breathing such that Q = Q max q(t). Eq S3 is reduced to its dimensionless form by multiplying and dividing Eq S3 with L 0 A 0 D a and − α ln(α) 1 − α , respectively, and using the following scaling parameters where, P e a and St a are the Peclet number for aerosols and Strouhal number for the airways, respectively. ϕ a and τ denotes the dimensionless aerosol concentration and time, respectively, while the quantities T a and T b represents the convective airflow time-scale and breathing time-scale, respectively. The expression of D a is based on the Stokes-Einstein relation 8 , where C s represents the Cunningham slip correction, T represents the ambient temperature, µ a denotes air viscosity, and d a denotes the aerosol diameter.
The dimensionless equation (equivalent to Eq 6 in the main manuscript), thus, obtained is used to analyse aerosol transport in the present study and is given by where, L ′ D represents the dimensionless form of aerosol deposition coefficient (L D ; see Section II B) and F a represents the total aerosol flux. These are expressed as follows -

B. Aerosol deposition models
The major mechanisms of aerosol deposition in the lungs have been identified in the literature as diffusion, sedimentation and impaction of the aerosols in the airways, as well as diffusion and sedimentation of the aerosols in the alveoli 5,9 . Different empirical models have been used in the past to estimate various depositions. However, these models need to be converted into a form consistent for use with the derived transport equation (Eq S5). This is discussed in the present section.
The probability of aerosol deposition in the airways by diffusion (P d ), sedimentation (P s ) and impaction (P i ) can be expressed following Yeh & Schaum 10 as Eq S8 can be re-written as where, the terms B d , B s and B i are the corresponding coefficients, and k d , k s and k i are the corresponding constants in the exponential functions for different deposition mechanisms as proposed by Yeh & Schaum 10 . Detailed expressions for the different deposition mechanisms can be found in the subsequent discussion. Taking the derivative of c a (Eq S9) with respect to x, we obtain - Equation S10 represents the droplet deposition flux in the airways. It is further converted to a dimensionally relevant form for use in the transport equation (Eq S3) as follows - The term L D represents the aerosol deposition coefficient which is determined using different empirical relations. The empirical relations are converted to a form relevant to Eq S11 and then reduced to their dimensionless forms for use in the final transport equation (Eq S5). These are discussed in the following sections for the various deposition mechanisms considered in this analysis.

Diffusional deposition in the airways
The probability of diffusional deposition (P d ) of the aerosols in the airways can be expressed following Yeh & Schaum 10 as The above equation can be simplified by expressing the coefficients in terms of effective magnitudes (B d , k d ) as follows where, B d and k d are determined as It is estimated that this simplification does not have any significant influence on the calculation for diffusional deposition (see Fig A). The simplified form is, as such, used in this analysis for calculation diffusional deposition in the airways. Using Eqs S12 and S13, we obtain - and where, R N and v N denotes the airway radius and airflow velocity of a particular lung generation, respectively. D a denotes aerosol diffusivity in air. Using this, droplet deposition in the airways due to aerosol diffusion is estimated as (refer to Eq S11) Conversion of Eq S17 to its dimensionless form gives us the following expression for dimensionless diffusional deposition of the aerosols in the airways -

Sedimentation deposition in the airways
The probability of deposition of the aerosols (P s ) due to sedimentation in the airways is expressed following Yeh & Schaum 10 as where, ρ a , g and ψ represents droplet density, gravitational acceleration and airway orientation angle considering horizontal as 90 • , respectively. Linearising the above equation using the approach followed in Eq.S9, we obtain - Aerosol deposition in the airways due to sedimentation can, then, be estimated as - Conversion of the dimensional deposition (L D,s ) to its dimensionless form gives us the following expression for dimensionless sedimentation deposition in the airways - where, S g is defined as the sedimentation parameter and expressed as

Impact deposition in the airways
The probability of deposition due to impaction (P i ) of the aerosols in the airways is given by Yeh & Schaum 10 as where, θ denotes the branching angle of the airways and St denotes the Stokes number The expression of P i is not in a form that can be directly linearized. As such, certain mathematical treatments are carried out in order to estimate the impact deposition.
Loss of aerosols i.e. the amount of aerosols deposited in one generation of the lungs can be determined based on the aerosol concentrations before and after the lung generation.
Mathematically, this can be expressed as In terms of lung generations, the above expression can be re-written as Differentiating with respect to generation number, we obtain - Converting the above derivative to a derivative in terms of x, we get - The impact deposition is estimated using the above expression as The dimensionless form of the impact deposition is obtained following a similar approach as in Sections II B 1 and II B 2. The dimensionless impact deposition, thus, obtained is expressed as

Diffusional deposition in the alveoli
Diffusional deposition of the aerosols in the alveoli is estimated using the following dimensionless expression - where, γ N denotes the fraction of alveolated area in the corresponding generation (see Table   B) and η d,alv denotes the diffusional deposition efficiency in the alveoli. η d,alv is expressed

Sedimentation deposition in the alveoli
Deposition of the inhaled aerosols due to their sedimentation in the alveoli are estimated using the following dimensionless expression - where, γ N and η s,alv denotes the fraction of alveolated area in the corresponding generation (see Table B) and sedimentation deposition efficiency in the alveoli, respectively. η s,alv is expressed as 5 η s,alv = 1 + min d s d eq , 1

Total deposition
The total amount of aerosols deposited at any lung generation is determined as the sum total of aerosol deposition through all deposition mechanisms. Mathematically, this is expressed as C. Drug molecule transport in mucus The corresponding 1D transport equation for the drugs deposited in the airway mucus (equivalent to Eq 7 in the main manuscript) is expressed as where, c d denotes the drug concentration in the airway mucus, Q m represents the volume flow rate of mucociliary clearance and D d denotes the diffusivity of drug molecules in the mucus layer.
The drug-laden aerosols deposited in the airway mucus serve as the only source of drugs in the lungs. The source term in Eq S37 is, therefore, equivalent in magnitude to the deposition term in Eq S1 (L D c a ) times the drug load in aerosols (ϕ l ). Mathematically, this is expressed as - where, ϕ l is defined as the amount of drug molecules contained by the aerosols per unit amount of the aerosols. Equation S37 is converted to a form in terms of N using Eqs S2 and A m = A m,0 (2 √ βζ) N in a similar manner as in Section II A as follows - The above equation is further reduced by multiplying and dividing by L 0 A m,0 D d and − α ln(α) 1 − α , respectively. The reduced equation is expressed as The following parameters are utilised to achieve the dimensionless form of the drug molecule transport equation in the airway mucus given by Eq S42.
where, ϕ d , P e d and St m represents the dimensionless drug concentration, Peclet number for the drug molecules and Strouhal number for the mucus layer, respectively. T m denotes the time-scale for mucociliary transport. Eq S42 is equivalent to Eq 9 in the main manuscript.
Drug diffusivity (D d ) is estimated using the Stokes-Einstein relation where µ m represents mucus viscosity and r d represents size of the drug molecules. The term F d in Eq S42 represents the total flux of the drug molecules and is expressed as

D. Implementation of the model and validation
The mathematical model discussed in Sections II A-II C is implemented for computational analysis using MATLAB ® . The governing transport equations(Eq S5 and Eq S42) are

Diffusion in Airways Diffusion in Alveoli Sedimentation in Airways Sedimentation in Alveoli Impaction
Whole lung Alveolar Region (c) (d)

III. PHYSIOLOGICAL BASIS FOR PARAMETER SELECTION
The magnitudes of different parameters used in the mathematical model are selected based on relevant physiological data. Physiological quantities pertinent to the lung model are tabulated in Table A. Other relevant physiological quantities are summarised in Table   C below.

IV. SUPPORTING RESULTS
A. Effect of aerosol size on drug deposition in the deep lungs  The increase in deposition with exposure time is linear (see Fig 4D in  to suggest ways to improve the drug delivery to the deep lung. Although the magnitude of drug deposition and drug concentration after a certain exposure duration can be determined by such extrapolations, it needs to be noted that this method is not a substitute for detailed simulations. Detailed simulations are still needed for drug retention calculation. The deposition characteristics can also change with variation in any one of the relevant parameters. Also, this knowledge does not provide information about the fraction of the inhaled aerosols that are deposited in the deep lungs. These informations can only be obtained from detailed simulations. Fig Fb shows the drug concentration within the lungs for various τ exp at the end of respective exposures. As expected, the drug concentration also increases due to larger aerosol deposition. In the upper airways (N < 18), mucuociliary transport clearance occurs simultaneously with aerosol deposition and as such, the effective drug concentration is the resultant of drug transport due to the deposition and clearance mechanisms. While drug concentration increases in these generations due to higher aerosol deposition, continuous mucus transport clears the drugs from these generations towards the 0 th generation and as a consequence, drug accumulates in the first few generations (leading to much higher ϕ d ) before being washed out of the lung.