A Pediatric Upper Airway Library to Evaluate Interpatient Variability of In Silico Aerosol Deposition

The airway of pediatric patients’ changes through development, presenting a challenge in developing pediatric-specific aerosol therapeutics. Our work aims to quantify geometric variations and aerosol deposition patterns during upper airway development in subjects between 3.5 months–6.9 years old using a library of 24 pediatric models and 4 adult models. Computational fluid–particle dynamics was performed with varying particle size (0.1–10 μm) and flow rate (10–120 Lpm), which was rigorously analyzed to compare anatomical metrics (epiglottis angle (θE), glottis to cricoid ring ratio (GC-ratio), and pediatric to adult trachea ratio (H-ratio)), inhaler metrics (particle diameter, dp\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${d}_{p}$$\end{document}, and flow rate, Q), and clinical metrics (age, sex, height, and weight) against aerosol deposition. Multivariate non-linear regression indicated that all metrics were all significantly influential on resultant deposition, with varying influence of individual parameters. Additionally, principal component analysis was employed, indicating that dp\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${d}_{p}$$\end{document}, Q, GC-ratio, θE, and sex accounted for 90% of variability between subject-specific deposition. Notably, age was not statistically significant among pediatric subjects but was influential in comparing adult subjects. Inhaler design metrics were hugely influential, thus supporting the critical need for pediatric-specific inhalable approaches. This work not only improves accuracy in prescribing inhalable therapeutics and informing pediatric aerosol optimization, but also provides a framework for future aerosol studies to continue to strive toward optimized and personalized pediatric medicine.


Introduction
Current strategies to deliver airways to the pediatric lung fall short of efficient airway delivery, necessitating new approaches to inhaler personalization that take into account the age of the child [1]. For children, as much as 80% of therapeutic inhaled aerosols end up in the oropharyngeal region of the upper airways and throat, where aerosols that do not deposit continue to the lungs [2]. Thus, pediatric lung delivery is often underestimated compared to adults, while off-target upper airway delivery is enhanced [1,3]. Anatomical variations between patients are suspected to influence fluid dynamics, aerosol transport, and ultimately therapeutic delivery; thus, the variations between adult and pediatric subjects have been long hypothesized to dictate inhalable therapeutic delivery efficiency [1]. The overall anatomy of the airway increases in size as children age, but additional developmental features further influence shapes, angles, and constrictions throughout the upper airways [4][5][6][7]. Development is a transient process in which a subject may exhibit a range of features, where patient characteristics (such as age, sex, height, and weight) are intertwined with anatomical parameters (such as trachea size) in ways that have yet to be decoupled to predict aerosol deposition [8]. Determining the ideal set of anatomical characteristics needed a priori to predict distribution of inhalable therapeutics would be an important step toward developing approaches for pediatricspecific inhalable therapeutic delivery.
Anatomic features have been shown to influence air and particle flow patterns that affect aerosol therapeutic delivery [9,10]. For example, after the narrowing of the vocal cords, the air flow is focused into a high velocity jet (i.e., "glottal jet"), which is influenced by the geometry variations between adult and pediatric subjects and leads to varied aerosol deposition [8,11]. However, the relationship between anatomical feature variation during development and resultant aerosol deposition lacks thorough investigation. Such knowledge would avoid unwanted off-target upper airway delivery and improve lung delivery for current pediatric therapeutics, as well as provide a framework for targeted upper airway delivery for future therapeutics. Additionally, pharmacokinetic models rely on such deposition information for accurate modeling and prediction of therapeutic response [12]. A meticulous investigation of anatomical geometry, fluid dynamics, and aerosol transport relationships is a critical piece in the holistic puzzle of how and why inhalable therapeutic delivery differs in pediatric subjects. The major challenge in experimental pediatric studies is the barrier of ethical concerns [3]; to circumvent this, computational fluid particle dynamics (CFPD) studies allow researchers to study the complex flow of air and aerosols within patient-derived models of the respiratory system, avoiding ethical concerns of direct in vivo testing of pediatric patients while utilizing a high-throughput approach. CFPD involves computer simulation to model the behavior of aerosolized particles in an anatomically-relevant three-dimensional geometry, providing information on flow profiles and deposition under specified boundary conditions. This enables broader interpatient analyses and can capture wider set of patient features, as opposed to use of a single or an idealized model. CFPD also allows for high-throughput studies and application of rigorous techniques such as principal component analysis (PCA) to be utilized [13,14]. Using these techniques, researchers can create detailed models of the respiratory system and simulate the deposition of aerosols at different ages and stages of development [8,[14][15][16].
In this work, we created a library of designated healthy pediatric upper airway models between 3.5 months and 6.9 years of age and generated CFPD simulated aerosol deposition data. We evaluated aerosol deposition as a function of three categories: 1) relevant clinical metrics (age, sex, height, and weight), 2) aerosol therapeutic device metrics (particle diameter d p and mouth flow rate Q), and 3) patient anatomical parameters typically not immediately available to a clinician (trachea size ratio (H-ratio), glottis and cricoid ring area size ratio (GC-ratio), and epiglottis angle ( E )). The impaction parameter, d 2 Q, was used in analysis to compare to literature mathematical models predicting deposition. Additionally, deposition profiles in four adult upper airway models were generated as a comparison to the pediatric subjects. Rigorous statistical analyses (i.e., multivariate non-linear regression (MVNLR) and PCA methods) were performed to investigate dominating influential parameters in aerosol deposition. Parameters d p and Q were repeatedly identified as dominating parameters in predicting aerosol deposition. Additionally, the anatomical metrics of GC-ratio and E and the clinical metric of sex were identified as key parameters in aerosol deposition through PCA analyses; age-and development-related variations are potentially incorporated into these decoupled parameters. To our knowledge, this work represents the largest pediatric upper airway CFPD study of this understudied age range to date. We hope it provides improved insight into deposition patterns in the developing airways of pediatric patients needed to advance therapeutic aerosols for this age range.

IRB Approval and Ethical Considerations
A dataset of healthy pediatric patient CT scans was originally obtained through a prior study exempt from consent under Internal Review Board (IRB) approval and is used retrospectively in this work. Subjects were exempt from consent under Category 4 of secondary studies on deidentified data. Access to these CT records was obtained by approval through the University of Delaware and Nemours Children's Hospital-Delaware IRB Review. No new data or information was collected from the human pediatric subjects involved in this work. Subject-specific CT scans have been deidentified such that the information known to the primary researchers are the CT scan data, patient age, sex, height, and weight. Subjects have been assigned an identification number for internal reference.

CT Reconstruction, Mesh Generation, and CFPD Simulations
CT reconstruction was performed as previously described [8,17,18]. Briefly, twenty-four CT scans of pediatric subjects between 3.5 months to 6.9 years old (14 male (M), 10 female (F), reported in Table I) were converted to 3D renderings and compared to 4 adult models previously rendered (see Fig. 1) [16]. Meshes of various element counts were constructed to confirm results were mesh independent. CFPD simulations were performed using Ansys Fluent 2021 R1 (Ansys, Inc., Canonsburg, PA, USA) in congruence with previously published work [8,17,18].  Briefly, a Brownian motion in-house C program enhanced Euler-Lagrange method and transition shear stress transport (SST) model were employed to simulate unsteady and incompressible airflow and the transport and deposition of spherical particles with negligible thermophoretic forces and a high particle-to-air density ratio [17]. Gravity was defined by generating an averaged vector based on the model wall and the greatest axial vector was defined by the scalar value for the gravitational constant, 9.81 m/s 2 . The particle velocity, diameter, and physical characteristics, as well as the fluid assumptions of values of Q, inlet pressure boundary conditions, the trapped wall particle-wall interaction condition, and triplicate file injection particle tracking method were assumed based on previous research [8,[17][18][19]. Of note, flow rates during tidal breathing for pediatric subjects typically range between 10 and 15 Lpm that will depend on the subject age [20,21]; however, peak inhalation flows can reach much higher. Often these peak flows could occur during distress or may be necessary for actuating dry powder devices [22,23]. Thus, we evaluated a wide range of flow rates spanning 10-120 Lpm for all models based on prior works [8,[24][25][26], while acknowledging that certain subsets of flow rates may be more statistically relevant for patients of different ages. For further information, see Supplemental Methods S1 and S2; mesh independence analyses are shown in Supplemental Figures S1-S7.

Anatomical Analysis
CT scans of the 24 subjects were reviewed by a radiologist and otolaryngologist to confirm the models were representative of healthy subjects. A previously established inhouse program was utilized to characterize the model and to evaluate the anatomical features of interest based on a published protocol [8,27,28]. The narrowest point of the glottis was compared to the cricoid ring area to determine the relative narrowest feature and to calculate the diameter ratio, GC-ratio. The angle of the epiglottis relative to the trachea ( E ) was measured to evaluate if it was more acute (< 30°) or more obtuse (> 35°) for the expected range of 20-45° for pediatric subjects [6,8,29]. Finally, the trachea was compared to an average adult trachea by determining the homothety ratio of the airways based on literature methods of comparing trachea diameters [30]. The pediatric trachea diameters were normalized to the average adult trachea diameter of the four adult models utilized in this study, resulting in a ratio describing the relative scale of the pediatric airway to the adult airway referred to as the H-ratio. Values of GC-ratio, H-ratio, and E against age for all subjects are reported in Table I.

Statistical Analysis
Two methods of statistical analysis were utilized during this study: seven MVNLRs and three PCAs were performed using R programming software (R Foundation for Statistical Computing, Vienna, Austria). The built-in multivariate function 'lm' and the CRAN package emmeans was utilized in conjunction with variables assembled in non-linear arrangements to create the MVNLR model and to visualize the results [31]. The R 2 value, regression model p-value, and individual variable p-values were assessed to evaluate the relative significance of each variable. The CRAN packages ggfortify, gridExtra, readxl, and dplyr were utilized to perform the PCA analysis. The optimal number of components to include was deduced by interrogating the scree plot of each PCA performed and explaining variance greater than 80%. The resultant significance of the overall analysis and key variables of interest are reported in Supplemental Tables S1-S3. One model was excluded from these analyses due to missing information in their clinical file.

Results
CT scans from 24 pediatric subjects were successfully rendered, and their geometric features were characterized by their H-ratio, GC-ratio, and E . These values, as well as subject age, sex, height, and weight, are reported in Table I. Additionally, graphs of age versus H-ratio, GC-ratio, and E are shown in Fig. 2, demonstrating no clear trends as a function of age. It was previously identified that aerosol deposition is mostly in the upper-most regions of the airways of interest, namely, the naso/oropharynx and supraglottis [8]. It can be seen that this is the case in these models, as shown in Fig. 3. Figure 3 shows pediatric subjects with a range of features, from a 1-year-old with typical "pediatric"-like features to a 6-year-old with a combination of "pediatric"like features and "adult"-like features, which exemplifies that development is a progressive process. Deposition appears restricted to a higher position in the youngest model ( Fig. 3a) while the 3-and 4-year-old subjects ( Fig. 3b and c) exhibit particle deposition at lower positions, with a growing number of concentrated deposition "hot spots" compared to 1-year-old subjects. This visualization indicates that a relationship between development and aerosol deposition could potentially be extracted.
To decrease the dimensionality of our deposition data, we next plotted the impaction parameter, d 2 Q, for each of the models, thus representing oropharyngeal deposition relationships. The relationship between d 2 Q and deposition fraction is often used in aerosol research to indicate variations between subjects and inform future practices [32][33][34]. We visualized the deposition of each model across a range of d 2 Q values, either individually or averaged into a group by parameters of age, sex, and cricoid ring area shape. Results were compared to literature values, as well as adult model values as a control. Based on evidence that upper airway development occurs as a subject ages, d 2 Q was represented as an average smoothed function for all pediatric models tested (blue curve, Fig. 4a) and broken out to each age group (Fig. 4b). For nearly all pediatric models, deposition started at 0% for low values of d 2 Q, a sharp transition to increased deposition occurred between 1E-2 and 1E-3, and deposition plateaued by 1E-4 at nearly 90%. By visual inspection, no obvious correlations emerge with age when subjects were grouped, and thus separated smoothed representations were also visualized (see Fig. 5).
These data aligned well with three pediatric mathematical models integrated from literature (black curves, Fig. 4) and served as key validation benchmarking for these studies. The mathematical models followed the expected trend with age, such that the model representing the youngest age group exhibited the highest deposition at lowest d 2 Q and the model representing the oldest age group exhibited the least deposition (see Fig. 4a). The mathematical model presented by Tavernini et al. represents young infants and predicts higher deposition than the models shown, as would be expected based on the theoretical prediction that younger subjects would exhibit higher deposition [35]. The Storey-Bischoff et al. model included a wider range of ages, 3-18, but a lower average age of 8.5 years old than that of the Golshahi et al. model, which represents subjects between 6 and 14 years old with an average age of 10 years old [36,37]. There appeared to be closer values to the data with the mathematical model presented by Storey-Bischoff et al., which represented a wide age range; however, 45% of subjects were in the age range of this work, the highest overlap of all the models considered [36]. Additionally, the Golshahi et al. model, which had the highest number of subjects outside of the age range considered for this work, exhibited lower predicted deposition compared to the data [37]. In stark contrast to the pediatric models, the adult models exhibited very low deposition in the range observed (green curve, Fig. 4), with a sharp increase at d 2 Q values around 1E-4, no plateau, and a maximum of 76% deposition. The literature-based mathematical model from Stahlhofen et al. for adult deposition exhibited similarities to the simulated data collected [38]. In comparison to the four adult model geometries utilized in this study, the pediatric models exhibited consistently greater deposition compared to adults for both the mathematical model and collected data.
We next evaluated different variables within our dataset, beginning by analyzing particle deposition across d 2 Q as an average smoothed LOESS (locally estimated scatterplot   Table I for 1.4, 3.3, 4.1, and 6.1, respectively. Airway reconstructions scaled to fit the height of the figure. The 1-yearold subject exclusively has deposition in the upper regions of the airway, while the 3-and 4-year-old subjects had deposition throughout. The number of deposition "hot spots" also increases with age smoothing) curve for each sex (see Fig. 6a). Both male (red) and female (blue) subjects exhibited an average increase in deposition starting around a d 2 Q of 1E2, with an inflection or midpoint around 1E3 and plateauing by 1E4, with males having a higher deposition fraction. Additionally, it can be seen from Fig. 6b and c that when the individual data points are visualized, there is evidently much more scatter amongst males compared to females, indicating that there might be more variability in the data of the male subjects that contributes to the qualitative analysis. While deposition varies between sexes in adults, geometry does not vary between sexes in the pediatric age range observed here [6,39]. Thus, the correlation to sex was unexpected; conversely it is expected that differences in anatomic geometries would be the main driving force of deposition variations across the developmental range evaluated here. In Fig. 7, we plotted deposition fraction and d 2 Q as a function of epiglottis angle E (Fig. 7a), trachea H-ratio (Fig. 7b), and GC-ratio (Fig. 7c). Each curve is color-coded, such that increasing values of each metric would correspond to a transition from violet to yellow. By visual inspection, no obvious correlations emerge among physical geometric metrics of anatomical features from this evaluation. These observations further necessitate rigorous statistical analysis.
MVNLR analysis was utilized to evaluate deposition and determine if there were statistically significant trends with the considered parameters, as described in "Methods." Fig. 8 shows the results of the MVNLR analysis, where the top plots in each panel indicate the MVNLR results, while the bottom plots show all of the individual data points along with a LOESS curve used to locally average the data. The three sets of parameters considered were those relating to Fig. 4 Evaluation of aerosol deposition as it varies by subject age. a Literature mathematical models developed across a range of ages compared to our pediatric data (blue) and adult data (green) with 95% confidence intervals in gray and the average age of each cohort in parenthesis in the legend in units of years. The data largely aligns with the literature models, with the qualitative trend expectedly showing an apparent decrease in average deposition with an increase in age. b The smoothed representation of deposition fraction across impaction parameter, d 2 Q, averaged together for each age group such that each curve represents 4 subjects. Representative mathematical models of deposition in pediatric subjects from literature are shown in black, with no clear correlation apparent with age. The mathematical models included are from Tavernini et al. [35] Storey-Bischoff et al. [36], Golshahi et al. [37], and Stahlhofen et al. [38] Fig. 5 Individual evaluation of aerosol deposition as it varies by subject age. The smoothed LOESS representation of deposition fraction across impaction parameter, d 2 Q, of each subject separated into panels by age group, including groups: a 1-year-old, b 2 year old, c 3-year-old, d 4-year-old, e 5-year-old, f 6-year-old, and g adult. Age is described by color map from dark blue to yellow for the pediatric subjects, while the four adult models are represented as four shades of red assigned at random clinical metrics (i.e., age, sex, height, and weight), those relating to anatomical development metrics (i.e., H-ratio, GC-ratio, and E ), and those relating to inhaler metrics (i.e., d p and Q). One MVNLR included all possible variables and an individual MVNLR was performed on each set of variables. Each combination pairing (i.e., inhaler and clinical, inhaler and anatomical, and anatomical and clinical) metrics were combined in three additional MVNLRs for a total of seven regressions. The resultant R 2 values of the model, p-values of the model, and individual p-values of each metric are shown in Supplemental Table S2 and S3. Notably, the highest R 2 value was from the regression including all variables, with an R 2 value of 0.7986 and a p-value of less than 0.0001. All variables were statistically significant in predicting deposition with p-values less than 0.02. Visualization of the regression with sex and GC-ratio for a range of d p values is a confirmation that these factors are significant in influencing deposition. Of the MVNLRs performed on each set individually, the p-values were all significant but the R 2 values were too low (less than 0.55), indicating that the models were not useful in predicting deposition. This is evident from the high variability within the dataset, as present in Fig. 8 bottom panels. Notably, the highest R 2 value was from the MVNLR including inhaler metrics, indicating an expected trend of deposition with d p and Q that is visualized in Fig. 8a and b, respectively. Additionally, when combining sets of variables, the R 2 were all below 0.7 and less useful in predicting deposition compared to the model including all variables. While these MVNLR results are indicative of which parameters are significant in predicting aerosol deposition trends, the relative significance and optimal combination of parameters was more accurately identified through PCA.
PCA is an unsupervised machine learning technique that allows for reduction in data dimensionality, allowing for additional insights to be generated of significant components in complex datasets, such as the one generated here [40]. Three PCAs were performed, one PCA including all possible variables (Fig. 9a), one including clinical metrics and inhaler metrics (Fig. 9b), and one PCA including metrics available Fig. 6 Evaluation of aerosol deposition as it varies by subject sex. a The smoothed representation of deposition fraction across impaction parameter, d 2 Q, averaged together for both sexes (red represents males and blue represents females). Confidence intervals of 95% are shown in gray. Each curve represents the average of the samples of males and females; thus, the male curve represents the average of 14 samples and the female curve represents the average of 10 samples.
The smoothed representation of deposition fraction across impaction parameter, d 2 Q, averaged together for b males and c females with the individual data points for all subjects. While both curves exhibit similar behavior, with a deposition increase at a d 2 Q of 1E2, an inflection at 1E3, and a plateau at 1E4, there is much more variability in males than females Fig. 7 Evaluation of aerosol deposition as it varies by subject anatomical feature geometric metric. a The smoothed LOESS representation of deposition fraction across impaction parameter, d 2 Q, for each subject such that colors ranging from dark blue to yellow indi-cate increasing value of the geometric metric: A epiglottis angle E , b trachea H-ratio, and c glottis and cricoid ring diameter (GC-ratio). The high degree of variability within the data is apparent from these curves to a clinician (Fig. 9c). In these plots, each point is a projection of the multi-dimensional data along the directions of the two principal components, while the labeled lines show how strongly each characteristic influences a principal component (i.e., loading). For the first PCA analysis of all possible variables (Fig. 9a), it was found that the first five components accounted for 86% of variability in the data (see Supplemental Figure S8 for the scree plots of the PCAs variances and Supplemental Table S3 for the principal components eigenvectors for each PCA). Loading of PC1 (35.1% of the overall variance) was dominated by clinical metrics such as age, height, and weight, while loading of PC2 (17.8% of the overall variance) was dominated by deposition and d p , as evidenced by the almost vertical vectors for deposition, d p , and Q (black, yellow, and blue, respectively). Indeed, the colorimetric gradient showing relative deposition across each data point similarly separates along PC2. Components PC4 and PC5 (10% and 9% of the overall variance, respectively) were most heavily loaded by age, H-ratio, E , and weight, indicating that these components are expressing the relationship between developing anatomical features.
For the second PCA, it was shown that the optimal number of components was three and that these components accounted for 81.4% of variability in deposition (see Fig. 9b). Here again, PC1 (42.2% of the overall variance) was loaded with clinical metrics such as age, height, and weight, while PC2 (24.9% of the overall variance) was loaded with d p and a lesser extent, Q and sex. For the third PCA, it was shown that the optimal number of components was two and that these components accounted for 81% of variability in the data (see Fig. 9c). Similar trends were noted in the variable loading for PC1 (59.0% of the overall variance), with sex dominating the loading of PC2 (22.1% of the overall variance).

Discussion
In this work, we evaluated 24 different pediatric airway scans from healthy children ranging from 3.5 months to 6.9 years old. This is a unique dataset comprising, to our knowledge, the largest CFPD pediatric upper airway data set to date and represents a transitional age where anatomical development is highly variable. As a subject ages from infancy to adulthood, deposition in the upper airways is generally expected to decrease with age based on in vivo data [36]; however, the limited information about this subset of pediatric subjects ages 3.5 months to 6 year prompted our investigation. From the mathematical models in literature, it is evident that younger subjects have higher deposition at lower d 2 Q (see Fig. 4). The mathematical model derived from infants (Tavernini et al.), adolescents (Golshahi et al., ages 6-14 and on average 10 years old) and a broad age group (Storey-Bischoff et al. with ages 3-18 average 8.5 with 45% of patients between 3 and 6 years old) exemplifies the overall broad alignment of our data with prior pediatric mathematical models [35][36][37]. Comparing our deposition as a function of diameter (Fig. 8a), 50% deposition efficiency will occur at ~ 2 µm for the mouth-throat regions modeled . In all panels, regression model is shown in dashed lines while the LOESS averaging is shown in solid lines. The inhaler metrics of d p and flow rate had a higher R 2 than the clinical metrics of age and GC-ratio here; this is in contrast to adults, where 50% deposition efficiency in the same region occurs at ~ 10 µm [41]. This shift in deposition efficiency supports those observations and suggests that the off-label prescription of therapeutics will not solve the issue of upper airway deposition. Instead, smaller particle manufacturing for pediatric-specific therapeutics is likely required [1].
The greater deposition in adults in both the collected data and mathematical model demonstrates that the overall trend of deposition decreases with age in the upper airways. Due to the increase in airway generation, volume, and geometry between the ages 3.5 months and 6.9 years, the resultant aerosol deposition in the upper airways is expected to be significantly greater in younger subjects compared to older subjects [42]. However, our results imply that the relationship between deposition and age within this dataset is much more nuanced. Indeed, all of the pediatric subjects had greater deposition compared to the adults modeled. The adult deposition was only apparent at the highest d 2 Q considered here and had significant overlap with literature, such that the LOESS smoothed model was nearly identical to the adult mathematical model (see Fig. 4). However, the youngest subjects of 3.5 months old presented data on deposition as a function of d 2 Q that was similar in trend to that of a 4-year-old (see Fig. 5a and d). Age was not identified as a significant variable with a PCA analysis, indicating that it did not account for significant variability within the data. It is possible that previous identifications of age as a significant factor in predicting aerosol deposition do not sufficiently decouple age from related metrics, such as GCratio, nor adequately account for inhaler metrics, and thus a falsely identified significance is placed on age. It could also indicate the high degree of anatomic variability that occurs for patients within this age range, further contributing to the lack of correlation strictly based on age. Regardless, these data indicate that the variable of age does not entirely account for the variation in deposition and employing a metric of age alone is not completely useful in distinguishing the need for pediatric-specific inhalation devices. Thus, it can be established that these geometric features are related to development, yet there is a complex relationship such that age is not an all-encompassing predictive parameter.
Interestingly, we observed more deposition in male subjects than female subjects, as indicated by the LOESS smoothed functions (see Fig. 6). Thus, particle deposition as a function of d 2 Q was separated by sex and indicated that, visually, males have more deposition than females in the upper airways of the pediatric subjects represented. This was unexpected, as it has been shown in literature that sex does not play a role in the overall deposition of pediatric subject upper airways in vivo [43]. Additionally, sex is not associated with variation in upper airway geometry (e.g., the diameter of the glottis, cricoid ring area, and trachea), [44] so it has been hypothesized that this would result in similar airway function and aerosol deposition between sexes of this age group [10]. Thus, sex was not expected to be a variable which statistically significantly influenced deposition in the range of pediatric subjects represented here. This with clinical and inhaler metrics, and C with only clinically available metrics. The relative significance of the parameters indicated that several variables were influential on resultant deposition, listed in Supplemental Table 3. For the PCA including all variables (a), the five important variables were d p , Q, GC-ratio, E , and sex. For the PCA shown in (b), three components dominated, being d p , Q, and sex. The PCA using only clinical metrics (c) showed a weak relationship with aerosol deposition, highlighting the importance of inhaler metrics over clinical may just be an artifact of our particular dataset; expanding to a greater number of patients would further corroborate if this is a representative trend. Notably, it has been seen that sex can be predictive of asthma severity [45], yet has no influence on intubation metrics [46], so there have been established instances of in vivo airway discrepancies in the influence of sex.
The d p and Q range utilized in this study encompass many commercially available aerosol therapeutics. For example, the Turbohaler® is a budesonide aerosol device for asthma treatment with particle size distribution roughly centered around an MMAD of 3 μm, and the Q suggested for use is 60 Lpm [47][48][49][50]. Deposition in these subjects at a particle size of 3 μm and the Q of 60 Lpm ranges from 0 to 100% (see Figure S9). This distribution of results indicates that pediatric subjects prescribed an inhaler based on age would have high variability in upper airway off-target delivery, even without the confounding effects of patient compliance. Bioavailability and pharmacokinetic modeling to identify variability in biological response to deposition and the resultant dosimetry requirements would be a promising step toward improved understanding of current inhalable device functionality [51] and prove insightful into improved future inhalable therapeutic designs.
The MVNLR model results reflect all parameters, suggesting that all variables were statistically significant when creating a model that included all parameters. The low R 2 values for other MVNLRs indicates that these models were not accurate at predicting aerosol deposition, and these combinations of metrics are insufficient to determine inhalable therapeutics for pediatric subjects. This is visualized for a clinical metric of age (see Fig. 8c) and an anatomical metric of GC-ratio (Fig. 8d). The insignificance of clinical metrics, such as age, indicates that broad generalizations cannot be made, such as all 6-year-olds having less upper airway deposition than infants. It is evident that inhaler metrics of d p and Q show the greatest increases in R 2 (0.2476 to 0.6960 for clinical metrics without and with inhaler metrics, respectively) suggesting that they are essential parameters in predicting deposition, which has been shown in previous work as well [8]. This work supports the previously identified increase in aerosol deposition between d p of 1-3 μm, suggesting that inhaler choice is more influential on upper airway deposition than simply the stage of development of the subject. Finally, while the inclusion of anatomical variables resulted in the MVNLR with the highest R 2 value and thus the most accuracy in modeling the data, the individual MVNLR had the lowest R 2 value; it is expected that, if development is intricately tied to deposition, these metrics were not successful in elucidating that fact.
Three PCAs were calculated-one which included all available parameters (Fig. 9a), one with clinical and inhaler metrics only (Fig. 9b), and one which only included clinically relevant metrics (Fig. 9c)-in an analysis intended to be similar to the series of MVNLR analyses. Across all of these analyses, we find that PC1s are loaded with physical attributes such as age, height, and weight, while loading of PC2s were dominated by deposition and d p , and to a lesser extent, Q and sex. This highlights the significance of interpatient variability that drives the largest source of variance in our data. For the first PCA (Fig. 9a), components PC1, PC2, and PC3 exhibit loading indicative of relationships between the parameters analyzed and the resultant aerosol deposition, such that the key parameters correlating to deposition emerge as Q, d p , GC-ratio, E , and sex. It is notable that d p is repeatedly identified as a key parameter, suggesting that it is the most significant factor in predicting aerosol deposition. Indeed, the corresponding eigenvectors in Fig. 9a for deposition and d p are almost overlapping. The relative loading of deposition compared to GC-ratio, E , and the other key parameters of principal component PC1 (i.e., d p and Q) indicate that the other parameters are less significant than GC-ratio and E . Thus, the first PCA results demonstrate that the most significant factor in predicting aerosol deposition is d p , followed by Q, GC-ratio, E , and lastly sex.
The second PCA analysis included only-clinically available metrics as well as inhaler metrics (Fig. 9b), both of which would be information currently available to clinicians prescribing inhalable therapeutics. The first three principal components accounted for 81.4% of variability in the data and were dominated by the relationship between deposition and d p , Q, and, to a lesser extent, sex. The third PCA including only clinically available metrics of age, sex, height, and weight was additionally performed. Here, the first three principal components accounted for 97% of the variability in the data; however, there is a weak relationship between these parameters and aerosol deposition, as indicated by the lack of similar deposition grouping in the plot (Fig. 9c). Within all three components, the parameters of age, height, and weight are more closely related to each other than they are to aerosol deposition or sex, a relationship innate to development. Overall, similar conclusions can be drawn from the PCAs as were drawn from the MVNLRs: clinical metrics alone are insufficient in predicting aerosol deposition and inhaler metrics markedly improve deposition prediction.
In this work, we have developed a library of upper airway models which addresses developmental anatomical geometry changes and predicts aerosol deposition in pediatric subjects; however, there are a few notable limitations with the use of CFPD. The particles are assumed to be in a dilute phase, and particle bounce and particle-particle interactions are not considered. Hygroscopic growth and humid environments of the lungs are not accounted for, and the particles are assumed to be water droplets but behave as rigid spheres. These necessary assumptions for CFPD simulation all limit the accuracy of these predictions to true inhalable therapeutics. Additionally, only inhalation was simulated and no particles were exhaled; a cyclic breathing profile with variable velocities and exhaled particle fractions would represent a more realistic respiration model. Our simulations utilized in this study are static; however, the glottis is known to be highly dynamic, and airway motion can be influential of aerosol deposition [52][53][54]. Exploring the implementation of a moving glottis could help to generate more realistic deposition data. Furthermore, the CT scan available area was inconsistent amongst subjects and was simulated with available anatomy. In literature, a subject-specific artificial cavity at the naso-oropharynx area has been shown to allow for fully developed flow and altered predicted deposition at the inlet, but such a section was not included in these models; future studies including this cavity could help to refine our current data at the inlet. Finally, aerosol delivery to this age range of children faces many practical challenges that are not well represented in our findings, such as pediatric subjects often rejecting the inhalation device, the need for use of passive nebulizers, lacking adherence to protocols, and their inability to follow instructions; in combination, this often leads to deposition in esophagus and stomach [2]. Improvements on physiological accuracy would serve to better understand this relationship, though the work stands on its own as a data library of otherwise underrepresented subjects.

Conclusions
Pediatric subjects exhibit drastically different responses to inhalable therapeutic delivery compared to adults, but ethical concerns have limited the extent to which the relationship between development and aerosol transport could be studied. In this work, we have presented a library of simulations based on CT scans of the upper airways of healthy subjects up to 6 years old through CFPD models of inhaled aerosols. It is evident from this work that these parameters are indicative of airway development, which is shown here to heavily influence aerosol deposition. The relationship between d p , Q, GC-ratio, E , and sex are the key parameters in predicting aerosol deposition. The anatomical features of the ratio between the glottis and cricoid ring diameter and the angle of the epiglottis vary as a pediatric subject develops and influence fluid dynamics and aerosol transport phenomena. The d p and Q are often constrained by the inhalation device and system, thus illustrating the need for inhalation devices to be developed intentionally for pediatric subjects. In this work, we have presented the largest pediatric upper airway CFPD library to date, which has enabled unparalleled statistical analysis revealing the intricately intertwined relationship between deposition and development and exemplifying the critical need for pediatric-specific inhalation devices.