A Partial Differential Equation Approach to Inhalation Physiologically Based Pharmacokinetic Modeling.

The heterogeneous nature of the lungs and the range of processes affecting pulmonary drug disposition make prediction of inhaled drugs challenging. These predictions are critical, as the local exposure cannot be measured and current inhalation physiologically based pharmacokinetic (PBPK) models do not capture all necessary features. Utilizing partial differential equations, we present an inhalation PBPK model to describe the heterogeneity in both lung physiology and particle size. The model mechanistically describes important processes, such as deposition, mucociliary clearance, and dissolution. In addition, simplifications are introduced to reduce computational cost without loss of accuracy. Three case studies exemplify how the model can enhance our understanding of pulmonary drug disposition. Specific findings include that most small airways can be targeted by inhalation, and overdosing may eradicate the advantage of inhalation. The presented model can guide the design of inhaled molecules, formulations, as well as clinical trials, providing opportunities to explore regional targeting.

WHAT IS THE CURRENT KNOWLEDGE ON THE TOPIC?

The heterogeneous nature of the lungs and the range of important processes (e.g., deposition, dissolution, and mucociliary clearance) make inhaled drug disposition challenging to predict. Current inhalation PBPK models do not capture all of these properties.

WHAT QUESTION DID THIS STUDY ADDRESS?

This study sought to develop an inhalation PBPK model accounting for the aforementioned processes as well as the heterogeneity in physiology and particle size.

WHAT DOES THIS STUDY ADD TO OUR KNOWLEDGE?

This study adds a sophisticated mathematical model for guiding the design of inhaled molecules and particle size distributions, as well as clinical trials, providing opportunities to explore regional lung‐targeting. Specific findings include: (i) inhalation leads to spatially varying free concentrations; (ii) most small airways can be targeted by inhalation; and (iii) overdosing may eradicate the advantage of inhalation.

HOW MIGHT THIS CHANGE DRUG DISCOVERY, DEVELOPMENT, AND/OR THERAPEUTICS?

This model can guide the design of inhaled molecules, and, furthermore, could aid the design and interpretation of preclinical/clinical studies.

The prospect of direct drug delivery to the target organ has made inhalation an attractive administration route for respiratory diseases. A well‐designed molecule and/or formulation could offer several advantages, including rapid onset of action and sustained local tissue concentrations.1 Equally important, high local concentrations may be obtained while minimizing the systemic exposure and, thus, potential systemic side‐effects.2 This is the aim of locally acting inhaled drugs: to obtain a lung‐targeted drug exposure.

Despite its widespread use, predicting inhaled pharmacokinetics (PK) and local pharmacodynamics (PD) has remained a long‐standing challenge for several reasons. For once, current methodologies do not allow for quantification of the target site concentration(s). For other administration routes, PK/PD evaluation is frequently facilitated by an assumption of equilibrium between free target site exposure and free plasma concentration, where the latter is easily measured. Obviously, this assumption cannot be made for a targeted drug after inhaled delivery. We note that indirect approaches, such as target occupancy, have been used.3, 4, 5 However, these methods are target specific, time demanding to develop, and do not currently provide a high regional resolution.

Pulmonary drug disposition is also recognized as complex because it involves a wide range of processes. The breathing pattern, airway morphometry, and particle size distribution (PSD) all combine to produce a heterogeneous deposition pattern throughout the lungs. Extrathoracic deposition also takes place; subjected particles may subsequently be swallowed and absorbed from the gastrointestinal tract. Prior to pulmonary absorption and/or target interaction, deposited particles have to dissolve in the epithelial lining fluid (ELF), which covers the epithelium. However, far from every particle in the tracheobronchial region is dissolved at its deposition site. Due to a self‐cleansing mechanism, mucociliary clearance (MCC), particles are transported toward the pharynx and eventually swallowed.6 The MCC transport rate varies across the lungs with a maximum velocity in the trachea.7 Further complexity arises from the highly heterogeneous physiology of the lungs with distinct differences between regions.8

All these factors render a legitimate question: is it correct to only consider a single free concentration in the lungs, or is there in fact a spectrum of free concentrations spanning the lungs after inhalation? In the latter case, how should molecules, formulations, and devices be designed to allow for a targeted delivery to the expected target site of the drug? Due to the vast complexity of inhaled drug disposition and the shortcomings of contemporary experimental methods, this poses a very challenging question to answer.

To help unravel the complexity of pulmonary drug disposition, all the various processes (deposition, dissolution, MCC, etc.) and the heterogeneous nature of the organ and particle distribution should be taken into account. Currently, there exist several different inhalation physiologically based pharmacokinetic (PBPK) models, each of which covers some but not all of these important phenomena.9, 10, 11, 12, 13

To model the regionally varying concentrations, expected due to, for example, the deposition pattern, MCC, and the heterogeneous physiology, a high spatial resolution is required. Although models discretized into as few as two (tracheobronchial and alveolar)9, 11 or three pulmonary compartments (thoracic, bronchiolar, and alveolar)13 have contributed to our understanding of inhaled drug disposition, their limited spatial resolution significantly affects the simulated result. Physiology and deposition vary richly across the lungs, and averaging these fail to capture important nonlinear phenomena, as conclusions from a lumped region are not necessarily valid across the entire space. Additionally, the simulation does not provide sufficient basis for evaluating targeting in specific pulmonary regions.

Recently, compartmental models with a higher spatial resolution have been published10, 12 addressing several shortcomings of previous models. However, these do not account for heterogeneity in particle size (i.e., only monodisperse particles are considered). As the dissolution rate for poorly soluble compounds is known to be largely affected by the PSD, it is also important to have a high spatial resolution in an additional dimension, particle size. Furthermore, a high particle size resolution enables in silico comparisons of different PSDs, which provides a tool to investigate different materials prior to conducting preclinical/clinical trials.

Systemic PK is another factor with important implications for lung‐targeted drug exposure, as theoretically demonstrated by, for example, refs. 9, 14. To describe the interaction between the lungs and the systemic circulation, it is important to allow for a back‐flow of drug from the blood to the lungs. Furthermore, regional differences in surface areas and blood flows are expected to affect the regional rate by which a molecule leaves/enters the lungs,15 accentuating the spatial heterogeneity in concentrations. Clearly, the sum of these features cannot be well described by unidirectional rate constants for pulmonary absorption rate, such as in refs. 12, 16, 17.

Given the aforementioned mechanisms, inhaled drug disposition can be summarized as: several concentration states, both position varying and invariant, with mass transport driven by concentration gradients. In addition, a two‐dimensional partial differential equation (PDE) describes a particle state, spanning size, and position. The coupling between particle and concentration states is nonlinear. Existing models could be regarded as highly discretized approximations of this very general description. We have identified two issues with approximating high resolution solutions to the model. The system of equations is stiff, and combined with its two‐dimensional nature makes it computationally demanding. In addition, the radius shrinkage rate of dissolving particles is singular as the particle radius tends to zero. Current models avoid numerical issues by compromising mass balance and adding computational burden.9, 13

In this article, we propose the first inhalation PBPK model, which treats the lung as a continuous heterogeneous organ, and also the particles as a continuously varying distribution over space and size. We address the issues related to constructing numerical approximations of solutions to the model. We reduce the system to a one‐dimensional PDE, and also remove the singular term in the shrinkage rate, significantly reducing the computational cost while maintaining mass balance. Three case studies are presented, wherein important features of the model are demonstrated, and different aspects of inhaled drug disposition are explored.

The presented model provides a sophisticated tool for understanding inhaled PK and can thereby guide the design of inhaled drug molecules, PSDs, and devices. Furthermore, it could potentially aid the design of clinical trials. In addition, the presented numerical improvements could also be used in gut absorption modeling, where high level of discretization is still used.13, 18

Method

A schematic overview of the model is provided in Figures 1 and 2. The former focuses on the lung and graphically explains different physiological parameters, whereas the latter shows the whole‐body PBPK model. Blood flows and tissue volumes of the included organs are presented in Table S1 in Appendix S1. The lung structure and airway geometry were obtained from Weibel19 with a functional residual capacity of 3,000 mL, Table S3 in Appendix S1. Interpolation is used to compute intermediate values.

Figure 1

Cross‐section of an airway showing the epithelial lining fluid (blue), epithelium (red), and subepithelium (light brown). Important physiological parameters are highlighted: (i) the cross‐sectional areas of the fluid, epithelial, and subepithelial layers (A f(x), A ep(x), and A sub(x), respectively), and (ii) the circumferences of the airway and epithelial layer as a function of their diameters (πD(x) and πD ep(x), respectively).

Figure 2

Structure of the whole‐body physiologically based pharmacokinetic model.

Particle deposition

Regional drug deposition was modeled using a typical path model, described in detail in ref. 20. Tidal breathing was used for the deposition calculations with a tidal volume of 1,000 mL, a breathing frequency of 15 breaths/minute, and the inspiratory:expiratory ratio was 1 without any pause between inhalation and exhalation. The airways were scaled as described in ref. 21 with the bolus volume set to 1,000 mL to reflect tidal breathing.

Extrathoracic deposition during oral inhalation was described according to Cheng.22 The model accounts for deposition in the airways by inertial impaction,23 gravitational sedimentation,24 and diffusion.25 The latter contains a modification under turbulent conditions (Re > 2,000), as described by ref. 23.

As for the computation of deposition in each airway generation, the well‐established closed form expression provided by Findeisen26 and Landahl27 is used. In short, volumes of inhaled air pass from generation to generation, depositing a fraction of available particles during each transition.

Inhalation PBPK model

Let p(x, r, t) describe the distribution of particles over lung depth x ∊ [x min, x max], particles radius and time t ∊ [0, T]. MCC and dissolution can both be modeled as nonlinear transport, moving particles along x and r, respectively. The first order PDE describing these two phenomena is where α(x) is the particle velocity along x at spatial position x due to MCC, and β(x, r, t) the rate at which a particle of radius r shrinks along r due to dissolution. We set the following sufficient, although not necessary, conditions for the system to be well‐posed: and β(x, r, t) ≤ 0 (i.e., there is no MCC at maximum depth, MCC only transports particles upward, and finally particles cannot precipitate).

The MCC term is described as proportional to the relative change in the cross‐sectional area of the airway A(x)7:

where α0 = 3.6 mm/minute is the MCC at the very top of the lung according to ref. 28. For well‐posedness, the equation is slightly shifted such that MCC is zero at the bottom, see Appendix S1. The change in particle radius can be derived from the Nernst‐Brunner equation29, 30 as where v diff is the diffusion coefficient, ρ particle density, v rr radius relative stagnant water layer, C s solubility, f u,f unbound fraction in the ELF, and finally C f(x, t) the concentration of drug in the ELF. We note that C f(x, t) varies over the depth of the lung, and that the equation is singular in r.

The mass of undissolved drug per unit x is

Consequently, the rate of drug exiting the lung to the gut is defined using the amount of drug and MCC velocity at the boundary:

Next, the drug mass transfer from dissolving particles is defined using the change in radius

The change in ELF concentration can be described as

where N(x) is the number of airways at depth x, A f(x) is the cross‐sectional ELF area, D(x) the airway diameter, P(x) the airway permeability (defined in Appendix S1), C ep(x, t) epithelial drug concentration, and K p,u,lung the unbound tissue‐plasma partition coefficient.

The rate of change in epithelial concentration is

where A ep(x) is the cross‐sectional epithelial area, D ep(x) is the epithelial layer diameter, and C sub(x, t) is the subepithelial drug concentration. Similarly, the rate of change in subepithelial concentration is

where A sub(x) is the cross‐sectional subepithelial area, R the blood‐plasma ratio, Q CO(x) is the cardiac output per unit lung depth, K p,lung is the tissue‐to‐plasma partition coefficient, Q br(x) is the bronchial blood flow per unit lung depth, C a(t) is the arterial blood concentration, and C ve(t) is the venous blood concentration.

Note that the latter two rows in Eq. (9) define the mass flow of drug to the traditional PBPK system described in Appendix S1. The model remains consistent with Figure 2 by setting Q br(x) and Q CO(x) to zero in the alveolar and tracheobronchial region, respectively. For multiple reasons, it can be beneficial to simulate these two regions separately. This, and the calculation of A f(x), A ep(x) and A sub(x) are also discussed in Appendix S1.

Numerical approximation

The main difficulty with constructing numerical approximation to the presented model is the singular particle dissolution rate (Eq. (3)). Additionally, approximation speed is slow due to stiffness, the two‐dimensional PDE (Eq. (1)), as well as the integral describing the rate of dissolved mass (Eq. (6)). In this section, we show how to mitigate these difficulties.

Rather than working with the radius, define the squared radius as a new spatial dimension z = r  2. The rate of change in squared radius for a particle is

and as such the dissolution is no longer singular. In addition, we find the dissolution is independent of z, this can be used to simplify the problem further.

Let us define the source (original) depth of particles currently at p(x, r, t) as

with initial condition x 0(x, 0) = x. In addition, because particles now move uniformly along z according to Eq. (10), we can define as the displacement along z of particles at position x at time t by

with initial condition δ(x, 0) = 0. Finally, let denote the relative thinning of particles due to transport along x:

with initial condition . Using Eqs. (11), (12), (13), we can describe the current number of particles as

where p 0(x, r) = p(x, r, 0) is the initial particle distribution.

The new PDE (Eq. (12)) describing the particle dissolution requires only a single spatial dimension, details on the approximation scheme are contained in Appendix S1. In addition, because both Eq. (11) and Eq. (13) have slow dynamics and are independent of the other system states, these can be precomputed to reduce the number of states. One could regard Eqs. (12) and (13) as applying the method of characteristics in only the z direction.

Finally, Eq. (14) is applied to the particle to fluid mass transfer (Eq. (6)) such thatwhere the integral term is a function of only x 0(x, t) and δ(x, t), and can, thus, be precomputed for speed. A detailed derivation is included in Appendix S1.

Results

To highlight important features of the model and explore different aspects of pulmonary drug disposition, we have conducted three case studies. The drug solubility (C s), the PSD, and the lung deposited dose (LDD) are varied, as specified in Table 1. All other drug‐specific parameters are held constant, as described in Table 2. The results are summarized in Figures 3 and 4. The former figure shows the total lung and plasma PK, whereas the latter focuses on spatial and temporal aspects of free epithelial concentrations or lung‐targeting.

Table 1

Settings for the case studies

Case study	C s (nM)	LDD (μg)	PSD
1: MCC	100	100	Large
2: PSD	250	50	Small, medium, large
3: Overdosing	100	100, 400, 1,600	Large

The PSDs are defined in case study 2.

C s, solubility; LDD, lung deposited dose; MCC, mucociliary clearance; PSD, particle size distribution.

Table 2

Drug‐specific parameters

Parameter	Value
Blood/plasma ratio	1
CLb a (L/hour)	70
CLp a (L/hour)	70
F	0.20
f u,p	0.75
f u,fluid	1
k a (hour−1)	0.6
K p,adipose	0.56
K p,gut	3.7
K p,hepatic	5.9
K p,lung	4.9
K p,poorly	2.3
K p,richly	3.1
K p,spleen	3.9
K p,u,lung	6.5
MW (g/mol)	250
ρ (g/cm3)	1
P app c (cm/s)	1.5 × 10−6
v diff b (cm2/s)	8.5 × 10−6
V ss a (L)	140

CLb, blood clearance; CLp, plasma clearance; F, oral bioavailability; f u,fluid, unbound fraction in the epithelial lining fluid; f u,p, unbound fraction in plasma; k a, oral absorption rate constant; K p, tissue‐plasma partition coefficient; K p,u,lung, unbound tissue‐plasma partition coefficient; MW, molecular weight; ρ, particle density; P app, apparent permeability; v diff, diffusion coefficient.

aBW = 70 kg. b v diff was calculated according to Stokes‐Einstein equation. cThe effective permeability was calculated from P app according to ref. 34.

Figure 3

(a) Lung (C lung) and plasma concentrations (C p) for simulations with the following settings: (i) mucociliary clearance (MCC) + wide particle size distribution (PSD; blue), (ii) MCC + narrow PSD (red), and (iii) no MCC + wide PSD (black). (b) C lung and C p for the following PSDs: (i) small PSD (blue), (ii) medium PSD (red), and (iii) large PSD (black). (c) Dose‐normalized C p for increasing lung deposited doses: (i) 100 (blue), (ii) 400 (red), and (iii) 1,600 μg (black). Solid and dashed lines represent lung and plasma, respectively. In addition, C lung is the mass of drug in undissolved particles, fluid, epithelium, and subepithelium, divided by the total lung volume.

Figure 4

Top row: Free concentrations in the epithelium for simulations with the following settings: (a) mucociliary clearance (MCC) + wide particle size distribution (PSD), (b) MCC + narrow PSD, and (c) no MCC + wide PSD. Middle row: Free epithelial concentrations for the following PSDs: (d) small PSD, (e) medium PSD, and (f) large PSD. Bottom row: Lung‐targeting, as defined by the ratio between free epithelial and free plasma concentrations, for the following lung deposited doses (LDDs): (g) 100, (h) 400, and (i) 1,600 μg.

The PDE was converted to ordinary differential equations using a forward difference approximation and simulated using the MATLAB ode23s routine, suited for stiff problems. The upper and lower lung regions were each discretized into n = 100 states resulting in a total of 10 + 8n states. On a modern desktop computer with a final time of 36 hours, a mean simulation time of 5 seconds was recorded, this figure grew to 26 seconds at n = 300. We note that the small particles in case study 1 were the most costly with 17 and 47 seconds at the two resolutions, respectively. The number of steps required for each simulation case did not vary significantly with n, although the computation time increased by approximately (n 2.3). The model code is available in Appendix S2.

Case study 1: Mucociliary clearance

This case study demonstrates the importance of including both MCC and a PSD. This is done by turning off the transport and using a narrow PSD, respectively. The narrow PSD is centered around the mass median diameter of the original distribution, and mimics the behavior in models, which only account for monodisperse particles. The wide and narrow PSDs are shown in Figure S3 in Appendix S1. The following three cases are considered: (i) MCC + wide PSD; (ii) MCC + narrow PSD; and (iii) no MCC + wide PSD. The drug is poorly soluble.

As can be seen in Figure 3 a, neglecting MCC has a pronounced influence on the PK: particularly on the lung but also on the plasma profile. The lung retention is erroneously predicted to be substantially longer in the absence of MCC. Distinct differences are also noted between a narrow and wide PSD. As the former does not account for the larger particles in the distribution, it will dissolve more rapidly leading to an earlier and more distinct drop in lung concentrations.

When looking at the spectrum of free epithelial concentrations, we note distinct differences with and without MCC (Figure 4 a–c). Without MCC, the model predicts more sustained levels of free concentrations in the epithelium. This phenomenon can be explained by Figure S2 in Appendix S1 showing the mass of undissolved drug m(x, t). In the absence of MCC, m(x, t) is solely influenced by drug dissolution. For slowly dissolving drugs, this creates an artificial drug depot, which gives rise to an unrealistically long drug coverage. Thus, neglecting MCC could ultimately lead to overestimation of the effect duration. In contrast, only accounting for a narrow PSD predicts a faster dissolution rate, thereby underestimating the free epithelial drug coverage and, thus, potentially also the effect duration.

Case study 2: Particle size distribution

The PSD is an important determinant of the deposition pattern and, for poorly soluble compounds, also the dissolution rate. In this case study, the influence of the PSD is evaluated by simulating the following PSDs: (i) small, (0.75, 0.15) μm; (ii) medium, (1.5, 0.3) μm; and (iii) large, (3,0.6) μm, where d = 2r is the particle diameter.

The PSD largely influenced both the lung and plasma profile (Figure 3 b). This can be attributed to the dissolution rate and the deposition pattern, both which are accounted for by the model. A faster dissolution and absorption process is observed for a smaller PSD. This is reflected by a more rapid decline in the lung profile as well as a faster absorption to plasma. The latter is characterized by an earlier and higher plasma peak. Clearly, this behavior is most pronounced for the small PSD.

As for the free concentrations in the epithelial (Figure 4 d–f) and subepithelial layers (Figure S1 in Appendix S1), we note a spectrum of free concentrations arises along the lung. Furthermore, distinct differences are noted between these two layers. This spatial variation is even more pronounced for the larger PSD, which dissolves more slowly. The smaller PSDs initially provide higher free levels as compared to the larger PSD. However, because of the slower dissolution process, the decline in the free levels is slower for the larger PSD, which thus provides a longer coverage at the expense of lower levels (Figure 4 d–f).

It is also noted that larger PSDs lead to higher extrathoracic deposition: 1.0%, 3.3%, and 10.1% for the small, medium, and large PSD, respectively. In case of oral bioavailability (i.e., F ≠ 0), this could lead to higher systemic exposure and, thus, reduced lung‐targeting.

Case study 3: Overdosing

In this case study, we evaluate how the dose affects lung‐targeting by varying the LDD: 100, 400, and 1,600 μg. Let us define lung‐targeting as the ratio between the free concentration in any pulmonary subregion (e.g., the epithelium) divided by the free plasma concentration (i.e., a ratio of 1 indicates no targeting).

Evaluation of dose‐normalized and non‐normalized plasma concentrations show a small nonlinearity, which is particularly pronounced at t > 24 hour (Figure 3 c, and Figures S4 and S5 in Appendix S1). Nevertheless, this would not be easily detected in a preclinical/clinical study as the plasma levels of the lower doses are at risk of falling below the limit of quantification. Coupled with individual variability, we stipulate that an extensive number of samples would be required to investigate such a nonlinearity, without which the PK would have been defined as dose linear.

A different situation arises when evaluating lung‐targeting. The extent of targeting markedly decreases with increasing doses, this is especially pronounced in the epithelium (Figure 4 g–i). For all investigated doses, the targeting vanishes (the ratio converges to 1) upon sufficiently long simulation times.

Discussion

Predicting free concentrations at relevant pulmonary target sites after inhalation is critical for human dose prediction and for design decisions in both drug discovery and development. In this article, we present a model that can guide such complex decisions by mechanistically describing the various processes and accounting for the heterogeneity in physiology and particle size. In contrast to earlier inhalation PBPK models, this model utilizes PDEs to describe both the particle size and the lung depth with high resolution. We have demonstrated that these features have a critical impact on the simulations and shown how the model can be used to help unravel complexities of inhaled drug disposition.

As expected, the importance of a high‐spatial resolution was highlighted by a spectrum of free concentrations that spanned the lung after inhalation with distinct differences between the epithelial (Figure 4 d–f) and subepithelial layer (Figure S1 in Appendix S1). This emphasizes our hypothesis that it would be erroneous to assume a single free concentration representing the entire lung and, equally important, that it is crucial to identify the pulmonary region(s) relevant for the effect. Failing the latter would hamper design decisions.

Another important observation resulting from the high‐spatial resolution was that lung‐targeting could be obtained in most small airways (Figure 4 g–i). In contrast, this could not be demonstrated earlier by a two‐compartment inhalation PBPK model.9 This shows that conclusions from a lumped region are not necessarily valid across the entire space. Nevertheless, the extent of targeting decreased with lung depth in the alveolar region, suggesting it is more challenging to target distal regions.

Simulations also indicate that overdosing risks to reduce the lung‐targeting (Figure 4 g–i), and possibly even eradicate the advantage of inhalation. Interestingly, a drastic decrease in targeting was observed at dose levels where this nonlinearity would have been challenging to detect from plasma observations in a preclinical/clinical study (Figure 3 c and Figures S4 and S5 in Appendix S1). Such a decrease might, thus, remain undetected. Hence, an optimal sampling scheme combined with a high analytical sensitivity is required to draw any conclusions of lung‐targeting based on plasma concentrations.

Choosing a suitable PSD is another important decision, whose significance increases for poorly soluble compounds. Simulations showed diverse spectra of free concentrations arose with varying PSDs. At the expense of a shorter duration, the smaller PSD initially provided higher free levels compared with larger particles (Figure 4 d–f). This highlights an important design decision: whereas larger particles might be advantageous for obtaining effect duration, this choice might be inappropriate if a high target site exposure is required to elicit a pharmacological response. Clearly, these decisions must be made on a case‐by‐case basis considering what target level and duration are needed to obtain the desired pharmacological effect. Furthermore, larger particles were more prone to deposit in the extrathoracic region, introducing additional complexity. This could lead to higher absorption from the gastrointestinal tract (if F ≠ 0) and reduced lung‐targeting.

Although contemporary methods do not allow for validation of all behaviors/states, some aspects of the model could potentially be validated by, for example, blood sampling in well‐designed studies. For instance, as discussed previously, detecting the nonlinear behavior of lung‐targeting from plasma measurements is challenging, potentially a study with a high time‐resolution and accurate pharmacodynamic readout could validate this behavior. Other phenomena may only be validated as the development of experimental methods progress, at some point we hope it will be possible to evaluate the spatial heterogeneity in concentrations predicted by the simulations.

Let us conclude the results analysis with an interesting observation. Usmani31 has, in fact, demonstrated that smaller particles achieve a higher and earlier plasma peak compared with larger particles. This behavior was reproduced by our simulations (Figure 3 b).

We recognize the challenges associated with the experimental methods providing some of the input parameters (e.g., solubility and airway permeability). In this article, the latter is estimated to be 10‐fold higher in the alveolar region based on its significantly thinner epithelial layer.8 There are currently no standardized methods for measuring the aforementioned values and it might not be feasible to represent the heterogeneity of the lung in a single assay.32 As such, it becomes interesting to revise the model as more information becomes available on, for example, regionally varying drug/physiology parameters, and/or pulmonary drug disposition mechanisms.

Also note that this model is developed for neutral small molecules relying on transcellular transport. Hence, it does not account for pH‐dependent solubility (acidic/basic compounds) or lung retention by lysosomal trapping (basic compounds). Neither does it consider active or paracellular transport, which can be incorporated if the investigated compound utilizes any of these transport mechanisms. Although only MCC is included as nonabsorptive clearance mechanism and we have neglected clearance by alveolar macrophages, the latter has been judged as less important for the investigated compound class.32 The aforementioned mechanisms might need consideration if the model is extended to other compound classes.

We emphasize that the presented model has applications beyond those exemplified by the case studies. For instance, it can guide the chemical design by evaluating how different drug‐specific properties affect the free concentrations in relevant regions. By subsequently relating these to expected potency values, it can provide directions for the design of a chemical series. Interesting extensions of the model could include pathophysiology and use of computational fluid dynamics for particle deposition.

At high resolutions, the simulation cost was approximately in proportion to that of matrix inversion ((n 2.376)33). The transformation to a one‐dimensional PDE reduced the complexity greatly, which would otherwise have been closer to ((n 2)2.376). In addition, without removing the singular term from the dissolution, simulation would have required significantly more iterations and restart of the algorithm. Importantly, the model also has applicability beyond the inhaled field. The mathematical description is general and may be extended to describe the important absorption from the gastrointestinal tract.

In conclusion, this article presents the first inhalation PBPK model utilizing PDEs to describe both the particle size and lung depth. As such, it provides opportunities for in silico evaluation of different materials and targeting of specific lung regions. Important applications include guiding the design of inhaled molecules and PSDs. Furthermore, it could aid the design and interpretation of preclinical/clinical studies. As the mathematical description is general, it may have utility in other fields, including absorption modeling from the gastrointestinal tract.

Funding

No funding was received for this work.

Conflict of Interest

E.B. is an employee of AstraZeneca and O.W. has no conflicts of interest to declare.

Author Contributions

E.B. and O.W. wrote the manuscript. E.B. and O.W. designed the research. E.B. and O.W. performed the research. E.B. and O.W. analyzed the data.

Appendix S1. Supporting text, tables and figures.

Click here for additional data file.

Appendix S2. Model code.

Click here for additional data file.