A mechanistic framework for a priori pharmacokinetic predictions of orally inhaled drugs

The fate of orally inhaled drugs is determined by pulmonary pharmacokinetic processes such as particle deposition, pulmonary drug dissolution, and mucociliary clearance. Even though each single process has been systematically investigated, a quantitative understanding on the interaction of processes remains limited and therefore identifying optimal drug and formulation characteristics for orally inhaled drugs is still challenging. To investigate this complex interplay, the pulmonary processes can be integrated into mathematical models. However, existing modeling attempts considerably simplify these processes or are not systematically evaluated against (clinical) data. In this work, we developed a mathematical framework based on physiologically-structured population equations to integrate all relevant pulmonary processes mechanistically. A tailored numerical resolution strategy was chosen and the mechanistic model was evaluated systematically against data from different clinical studies. Without adapting the mechanistic model or estimating kinetic parameters based on individual study data, the developed model was able to predict simultaneously (i) lung retention profiles of inhaled insoluble particles, (ii) particle size-dependent pharmacokinetics of inhaled monodisperse particles, (iii) pharmacokinetic differences between inhaled fluticasone propionate and budesonide, as well as (iv) pharmacokinetic differences between healthy volunteers and asthmatic patients. Finally, to identify the most impactful optimization criteria for orally inhaled drugs, the developed mechanistic model was applied to investigate the impact of input parameters on both the pulmonary and systemic exposure. Interestingly, the solubility of the inhaled drug did not have any relevant impact on the local and systemic pharmacokinetics. Instead, the pulmonary dissolution rate, the particle size, the tissue affinity, and the systemic clearance were the most impactful potential optimization parameters. In the future, the developed prediction framework should be considered a powerful tool for identifying optimal drug and formulation characteristics.

1 Derivation of the PDE model 1

.1 Derivation of the dissolution model
The Noyes-Whitney equation [1] describes the dissolution flux dW dt in terms of properties of the dissolving particles and the dissolution medium, where D is particle diffusivity, SA particle surface area, h height of the diffusion layer, C s particle solubility and C flu concentration of dissolved substance in the medium. Through geometric assumptions on particles, this equation can be turned into a differential equation describing the change of volume of a dissolving particle. We assume particles to be spherical in shape, with radius r, surface area SA = 4πr 2 , volume s = · 4 3 πr 3 and mass W = ρs. Furthermore, as suggested previously [2], we assume the height of the diffusion layer to equate particle radius, h ≈ r.
Since parametrizing the model in terms of radius r leads to a singularity of the dissolution model when r 0, in contrast to [2] we choose particle volume s = · 4 3 πr 3 as a size descriptor instead of particle radius. Differentiating the particle mass equation, and equating Eqs. (1) and (2) yields We opt to parametrize the dissolution model in terms of maximum dissolution rate k diss = D · C s rather than diffusivity D, since dissolution rate can be identified more directly from in vitro experiments (see Section Evaluation of dissolution model against in vitro data). The resulting dissolution model reads d(s, C flu ) = 4π k diss ( 4 3 π) 1/3 ρ · 1 − C flu C s · s 1/3 , ds dt (t) = −d s(t), C flu .
The concentration of dissolved substance, C flu , also changes during dissolution. These processes are coupled in the PDE model described below.

Derivation of the mucociliary clearance model
As explained in the main text, a continuous representation of airway radius r(x) depending on location x within the conducting airways is derived by interpolation. Using the Hofmann/Sturm model v = 0.12553 cm min · d 1 cm

2.808
, we obtain a location-dependent mucociliary clearance model for a particle at location x(t) at time t:

Individual and population states
Physiologically-structured models describe the time evolution of a set of individuals/particles, each exhaustively described by a vector of characteristics called state, denoted z, and which changes over time. The time evolution of the state of any individual is assumed to be governed by a law G, i.e. dz dt (t) = G(t, z(t)), z(0) = z 0 .
Assuming that a population consists of a large number of individuals, it is natural not to describe each single individual but rather the time evolution of a density ρ(t, z) of individuals over the state space. In this representation, the total number of particles is given by and the number of particles within a particular subregion ω of the state space is given by For such a domain ω, we set ω(t) = {z(t) : z 0 ∈ ω}. Assuming that the number of individuals is conserved in the state space, we obtain d dt N ω(t) (t) ≡ 0 as long as ω(t) does not touch the state space boundary. From this expression, a so-called continuity equation can be derived (see [3]): ∂ t ρ(t, z) + div z G(t, z)ρ(t, z) = 0. (3)

Derivation of physiologically-structured population models (PSPMs)
In our application context, the population consists of inhaled undissolved drug particles of different sizes, deposited at different locations within the conducting airways or within the alveolar space. The number of particles can only change if particles are (i) cleared to the GI tract by the mucociliary elevator (mucociliary clearance beyond the trachea, x(t) = 0) or (ii) completely dissolved (s(t) = 0).

Conducting airways
The particle state z = (x, s) ∈ [0, x TB ] × [0, s max ] can change by mucociliary clearance or dissolution (illustrated in Fig 1): , and Eq. (3) yields the location-and size-structured bronchial PSPM at time ! # > ! " Figure 1: Phase plane representation of a drug particle in the conducting airways. Each particle is characterized by its location and size. Over time, particles move within this two-coordinate system until they are either cleared to the GI tract or completely dissolved.

Mass balances
When coupling the PSPM models to equations for dissolved drug in lining fluids, the number of molecules (not particles) have to be conserved during dissolution and mucociliary clearance. This model feature is ensured by deriving dissolution and mucociliary clearance rates directly from the PSPMs, which is shown in the following. The number of undissolved molecules in the conducting airways / the alveolar space are given by  We illustrate the derivation for the conducting airways, using integration by parts at step ( * ):

Notation
We consider a uniform time discretization step ∆t > 0, a location discretization and a size discretisation 0 = s 1/2 < ... < s L+1/2 = s max These discretization points are understood as vertices of mesh elements (k, l) = [x k−1/2 , x k+1/2 ] × [s l−1/2 , s l+1/2 ] within which unknowns (approximations of ρ br , C br flu , etc.) are defined; they appear in the discretization of the location-and size-structured model in the conducting airways. The same size grid is also used when discretizing the size-structured model in the alveolar space. Furthermore, we define the center (x k , s l ) of mesh element (k, l) from the above discretization points, We use the following notation: • ∆s l := s l+1/2 − s l−1/2 (size length of mesh element (·, l)); we also define ∆s l+1/2 := s l+1 − s l (this expression will appear later during computations) • Abbreviations for location-structured physiology in conducting airways: λ k := λ mc (x k ), r br k := r br (x k ), q br k := q br (x k ), a br flu,k := a br flu (x k ), a br tis,k := a br tis (x k ) • ρ br,n k,l as the numerical approximation of ρ br (t n , x k , s l ) • ρ alv,n l as the numerical approximation of ρ alv (t n , s l ) • C br,n flu,k as the numerical approximation of C br flu (t n , x k ) • C br,n tis,k as the numerical approximation of C br tis (t n , x k ) • C alv,n flu as the numerical approximation of C alv flu (t n ) • C alv,n tis as the numerical approximation of C alv tis (t n )

Upwind discretization of physiologically-structured population equations
Upwind discretizations, i.e. non-centered finite difference approximations depending on the flow direction, are well tailored to PSPMs, resulting in stable discretizations as long as the timestep ∆t is small enough (called a CFL condition). The upwind discretization of the conducting airway PSPM is given by .., L} (with ρ br,n K+1,l = ρ br,n k,L+1 = 0, i.e. no inflow condition). Similarly the upwind discretization of the alveolar PSPM Within this framework, the number of undissolved drug molecules is approximated by

Implicit discretization of linear processes
Recognizing that all processes except for dissolution and mucociliary clearance are linear, we propose an implicit discretization to ensure unconditional stability of these other processes, too. The numerical scheme is formulated in terms of local amounts (in bronchial/alveolar fluid/tissue) rather than concentrations. To this end, we define Furthermore, it will be useful to define PS br k := ∆x k 2πr br k P app (permeability-surface area product at k-th location grid cell) Q br k := ∆x k q br k (perfusion of k-th location grid cell).
To arrive at a numerical scheme formulated on the computational grid, integrals are discretized as follows:  and noting that these two terms are matched in the equations for dissolved drug in the lining fluid and of cleared drug, we can conclude that mass is conserved during dissolution and mucociliary clearance. An analogous computation shows mass conservation during dissolution in the alveolar space. Mass balance was checked systematically during all simulations shown.

Projections onto the computational grid
Deposition patterns, as well as several parameters used in the PSPMs, are not resolved at the same scale as the computational grid. Therefore, a projection step is necessary prior to being able to integrate these quantities into the model.

Deposition patterns
Deposition data are given for each airway generation g 1 , ..., g K and for a fixed set of reference particle sizes S 1 , ..., S L , resulting in a discrete deposition pattern (D k,l ). The dose should be conserved, equivalent to conservation of number of molecules, but not number of particles.
We proceeded as follows (see Fig. 2 for an illustration): • We define a region S ε k around S k , given byS ε k = [S k − ε, S k + ε], with small ε such that all such regions are disjoint.
• From the discrete values D k,l , we define a continuous function • We define the initial condition on the computational grid by

Per-generation parameters
For a per-generation parameter (e.g., airway radius, blood flow, ...), generically denoted P , we construct a location-resolved representation using the previous construction only in the location coordinate, i.e.: • From the discrete values P k , we define a continuous function such that g k P (x)dx = P k .
• We define the location-resolved representation on the computational grid by

Evaluation of dissolution model against in vitro data
We evaluated the dissolution model against in vitro data from a dissolution study [4], where the authors evaluated the dissolution kinetics of fluticasone propionate and budesonide particles with defined particle sizes (see Table 1). Based on the in vitro data, we compared different dissolution models: • a first-order dissolution model (estimated empirically; size-independent) • an unsaturable dissolution model (formally corresponding to C s = +∞ in the dissolution model) • saturable dissolution models with different solubilities The results are shown in Fig. 3. A particle size-dependency is clearly visible, as well as a saturation effect. Among the different saturable dissolution models, the parametrization using in house data resulted in a qualitatively better description than the values reported in [4].

Location within airways
Particle size

Representation of lining fluid height in conducting airways based on literature data
Different values for the thickness of the lining fluid layer in the conducting airways have been reported. After reviewing the literature, we concluded that the linear relationship shown in Fig. 4 Table 1: Aerodynamic and geometric particle sizes corresponding to the experimental protocols of [4]. Particles within defined ranges of aerodynamic particle sizes were obtained from different stages of Anderson cascade impactors (ACI). For simulation of dissolution kinetics, we took the geometric diameter corresponding to the mean aerodynamic particle diameter within each impactor stage. assumptions (≈ 1.2 mL) was smaller than the ones given in the literature (10-70 mL).

Evaluation of Usmani data
As stated in the main text, we could not reproduce the fluticasone propionate exposure indices reported by Usmani et al. [9] based on the provided study information. Here we provide full details for this statement.  Table 2. In conclusion, the reported AUC value is approximately 2-4 times larger than what could be reasonably expected. Accordingly, C max values are also much higher than predicted by the PDE model. 4 Simulation of pulmonary deposition patterns 4

.1 Settings in MPPD software
In order to predict the pulmonary deposition patterns, the MPPD software v2.1 was used [11]. This software allows to predict the generation-dependent pulmonary deposition of inhaled particles, where generations 1-17 represent the conducting airways (generation 1 = trachea) and generations 18-25 the alveolar space. Three types of input data are required in the MPPD software: (1) airway morphometry, (2) particle properties, and (3) exposure condition, as outlined below. The MPPD software was only applied to simulate the deposition patterns but not used to investigate the clearance of particles from the lung.
Airway morphometry. For all predictions performed with the MPPD software, the airway morphometry was represented by the human "Yeh/Schum 5-Lobe" model [12]. The inhalation flow characteristics were assumed to be represented by uniform expansion of the lung so that consequently also the inhalation and exhalation flow were constant over time.  residual capacity (3300 mL) and upper respiratory tract volume (50 mL) were used [13].
Particle properties. The inhaled particle properties were defined based on the information in the respective publications, or alternatively for the respective inhalation device (references are provided in Table 3). For all predictions, the density of the particles was set to 1 g/cm 3 ; and the particle diameter was defined as the mass median aerodynamic diameter, which is typically provided in literature. As described in the main manuscript, the difference between aerodynamic and geometric diameters was accounted for, such that the real surface area could be used as an input parameter to the dissolution model. The MPPD software was only used to predict pulmonary deposition patterns of monodisperse particles. To predict the deposition patterns for the monodisperse gold/polystyrene particles (Study I) and the inhaled monodisperse fluticasone propionate particles (Study II), this information was sufficient. Whenever pulmonary deposition patterns of a particle size distribution were required (Studies III/IV), these were generated in a two-step approach. First, all relevant monodisperse particles size bins of the particle size distribution were simulated as monodisperse particles with the MPPD software. In a second step, the complete deposition pattern was calculated by normalizing the deposited amount per particle size bin by the dose in this respective bin. The two additional options of the MPPD software, namely the "Nanoparticle Model" and "Inhalability Adjustment" were not applied to predict the deposition patterns.
Exposure conditions. The exposure scenario was set to constant exposure and the body orientation during the inhalation process was assumed "upright". Furthermore, for all predictions, it was assumed that the breathing scenario was represented by oral breathing, which is the typical inhalation route for drugs delivered to the lungs. Breathing frequency, tidal volume, inspiratory fraction as well as pause fraction were all defined based on the inhalation flow properties provided in the respective publications (see Table 3).

Adaptation of deposition patterns for asthmatic patients
Since the MPPD software predicts deposition patterns in healthy volunteers, it cannot directly be used to predict deposition patterns in asthmatic or COPD patients. In these patients, due to narrowed airways, deposition is more central in comparison to healthy volunteers. Whenever patients were considered in a study rather than healthy volunteers, deposition patterns had to be adapted adequately. To this end, the fraction of the inhaled dose deposited in any specific airway generation was increased by an adjustment factor such that the deposited fraction of the lung dose in the alveolar space was 2-fold lower than in healthy volunteers. This number was derived from published data on conducting airway to alveolar deposition ratios [17].  5 Generation of in-house data 5.1 In vitro solubility determination in surfactant containing medium For in vivo relevant characterization of the drug solubility, the surfactant-containing medium Alveofact R , a commercially available product, was taken. Alveofact R contains phospholipids obtained from bovine lung (i.e., surfactants) and is available as dry powder ampoules ready for reconstitution. As reconstitution medium, a 0.1 mol/l sodium dihydrogenecarboante buffer with pH 7.4 was used. A suspension with 50 mg/ml Alveofact R was produced according to the information and instruction for use of the commercial product. At these concentrations, Alveofact R forms a micellar system. 1 mg of drug (either budesonide or fluticasone propionate) is suspended in 1 ml of this medium and shaken for 24 h at 37 • C. Afterwards, the suspension is filtered with a commercially available Whatman Mini-UniPrep syringeless filter containing a 0.45 µm filter membrane out of glass microfibers. As the micelles pass this membrane and as the concentration of phospholipids is too high to be directly injected in the HPLC system for analysis of the solubilized amount of drug, the micelles are destroyed by adding DMSO in a 1:1 ratio to the filtered micellar solution.
The phospholipids can be separated by an additional 5 -10 minutes centrifugation step. A small aliquot of the remaining solution is taken and injected into a HPLC system for quantitative analysis of the solubilized amount of drug.

Blood to plasma ratio determination
To determine the Blood:Plasma (BP) ratio, the respective amount of the drug (i.e., fluticasone propionate) was added to 490 µL human blood and to 490 µL plasma samples to obtain a drug concentration of 10 µM. Both the plasma (plasma sample #2) and the blood samples were incubated with the drug for 15 minutes at 37 • C (n=3). Afterwards the blood sample was centrifuged at 3000 rpm to separate the blood cells from the plasma sample (plasma sample #1). The respective plasma concentrations were determined by MS-based analysis. In a last step, the BP ratio was calculated by dividing the drug concentration in plasma sample #2 by the drug concentration in plasma sample #1. To ensure quality of the measurement, the degree of hemolysis was determined and considered negligible for all BP experiments. In addition, a control experiment without any drug was performed in parallel to determine the hematocrit of all three samples.