Toward Personalized Salbutamol Therapy: Validating Virtual Patient-Derived Population Pharmacokinetic Model with Real-World Data

Interindividual variability, influenced by patient-specific factors including age, weight, gender, race, and genetics, among others, contributes to variations in therapeutic response. Population pharmacokinetic (popPK) modeling is an essential tool for pinpointing measurable factors affecting dose-concentration relationships and tailoring dosage regimens to individual patients. Herein, we developed a popPK model for salbutamol, a short-acting β2-agonist (SABA) used in asthma treatment, to identify key patient characteristics that influence treatment response. To do so, synthetic data from physiologically-based pharmacokinetic (PBPK) models was employed, followed by an external validation using real patient data derived from an equivalent study. Thirty-two virtual patients were included in this study. A two-compartment model, with first-order absorption (no delay), and linear elimination best fitted our data, according to diagnostic plots and selection criteria. External validation demonstrated a strong agreement between individual predicted and observed values. The incorporation of covariates into the basic structural model identified a significant impact of age on clearance (Cl) and intercompartmental clearance (Q); gender on Cl and the constant rate of absorption (ka); race on Cl; and weight on Cl in the volume of distribution of the peripheral compartment (V2). This study addresses critical challenges in popPK modeling, particularly data scarcity, incompleteness, and homogeneity, in traditional clinical trials, by leveraging synthetic data from PBPK modeling. Significant associations between individual characteristics and salbutamol’s PK parameters, here uncovered, highlight the importance of personalized therapeutic regimens for optimal treatment outcomes.


Introduction
Therapeutic response varies considerably among patients, with drugs being pharmacologically effective in some while ineffective or even toxic in others [1,2].The precise determination of a safe and effective drug dose is directly dependent on understanding the PK and PD properties of that particular drug [3].Indeed, PK variability plays a crucial role in the success of drug treatments.Patient-specific factors such as age, weight, body mass index (BMI), gender, hormonal status, race, ethnicity, renal and hepatic function, genetic polymorphisms, disease status, concomitant therapies, smoking, and dietary habits can all contribute to variation in drug disposition [1,4].
As pharmacometrics advances as a cornerstone of precision medicine, there is a growing interest in understanding the impact of these factors on the PK profile of drugs.

Virtual Patient Data Collection
The physicochemical and PK properties of salbutamol were estimated using ADMET Predictor ® (Version 10.4; Simulation Plus Inc., Lancaster, CA, USA).Its chemical structure was drawn in MedChem Designer (Version 5.5; Simulation Plus Inc., Lancaster, CA, USA), according to its SMILES, and then imported into ADMET Predictor ® .
PBPK models developed in GastroPlus software (Version 9.8.3;Simulation Plus Inc., Lancaster, CA, USA) were used to generate virtual patient PK data.These models were built using the values predicted by ADMET Predictor ® .All patients were virtually treated with 600 µg of salbutamol DPI (3 successive inhalations of 200 µg).Thus, the dosage form selected in the software was PL:IT powder, which represents the delivery of inhaled drug in solid phase, as a single dose.The PK profile of salbutamol was modeled in different PBPK models with specific individual characteristics, namely, age, weight, race, and gender.In particular, subjects aged 5, 10, 20, 30, and 65 years were included.The body mass index (BMI) scale (BMI of 18.5-24.9 is normal, BMI of 25-29.9 is overweight, and BMI ≥ 30 is obese) was employed to determine weight.American, Japanese, and Chinese race groups were examined in this study.All patients were healthy, meaning that they had no adjacent pathological conditions.Detailed characteristics of these individuals are summarized in Supplemental Table S1.
Parameters such as bioavailability (Fa, fraction absorbed; FDp, fraction of the drug concentration in the portal vein; and F, fraction of the drug concentration in blood), maximum plasma concentration (C max ), time required to maximum plasma concentration (T max ), area under the curve (AUC), and maximum concentration in the liver (C max Liver ), were derived from ADMET Predictor ® .The drug disposition-based parameters were simulated in the virtual patients during 12 h.Quantitative and visual (plots) outputs of the drug concentration profile were extracted to establish a dataset for the subsequent phase of the study.
All PK parameter values were compared with data from the literature, and the plasma concentration profile plots were visually inspected to evaluate and validate all PBPK models.

Noncompartmental Analysis of PK Data
To establish initial PK metrics, a noncompartmental analysis (NCA) of the virtual dataset was performed in PKAnalix 2023R1 (Lixoft, Antony, France).The integral method employed was linear trapezoidal linear, with equal weighting assigned to each data point.As an acceptance criterion, we utilized adjusted R 2 to select the number of points for the terminal phase, enabling the determination of the slope of the linear regression for calculating λ z (elimination rate constant).Subsequently, all PK parameters were determined.

Population Pharmacokinetic Modeling
The virtual salbutamol data were inputted to build a nonlinear mixed-effects model (NLME) through estimation by maximum likelihood using the stochastic approximation expectation-maximization (SAEM) algorithm in Monolix Suite 2023R1 (Lixoft, Antony, France).The determination of conditional means and standard deviations for individual population PK parameters involved the utilization of Markov chain Monte Carlo (MCMC) convergence assessment.The objective function value (OFV), expressed as −2 × log likelihood (−2LL), and the Akaike information criterion (AIC) and Bayesian information criterion (BIC), were determined to select the best popPK model.Significant reductions in these parameters indicated better model fitting.

Structural and Statistical Models
Several PK models for extravascular administered salbutamol were investigated (one-, two-, or three-compartment models with first-or zero-order absorption with/without lag time or transit and with linear or Michaelis-Menten elimination).PK models correspond to a system of ordinary differential equations (ODEs) which describes transfers between compartments and elimination from the central compartment.For the observation model, the default combined 1 error model (additive and proportional term) were initially applied.For the individual model, we also keep the default: log-normal distributions for PK parameters.In a preliminary analysis, we ignored the inclusion of covariates in the popPK models.The optimal combination of structural and statistical models was systematically selected by evaluating the OFV, model diagnostics, and ensuring biologically plausible and accurately population mean estimates (as measured by relative standard error (RSE)).

Model Evaluation and Covariate Selection
After running several popPK models, the assessment of the model fitting was based on various evaluative metrics.BIC, the precision of estimates (RSE), and the goodnessof-fit (GOF) allowed the final decision for selecting the best popPK model.The model with the lowest corrected BIC (BICc) value was selected.In addition, complementary criteria including AIC, to estimate the model quality, and OFV, for identifying smaller values indicating superior fit, were also employed.The estimated population parameters' standard errors and the random effects error models' standard errors were also computed.Further, several diagnostic plots were used to visually test the model's fit, including observation vs. individual predictions (IPREDs), population-weighted residuals (PWRESs) vs. time/predictions, individual weighted residuals (IWRESs) vs. time/predictions, and the distribution of the empirical and theoretical standard Gaussian probability density functions (PDFs) against IWRESs and normalized prediction distribution errors (NPDEs).A visual predictive check (VPC) using a 90% prediction interval was performed to graphically assess misspecifications in structural, variability, and covariate models.
Age, height, and weight were the continuous independent covariates tested to be included in the popPK model.Gender and race were designated as categorical covariates.The selection of covariates was based on the model proposed by the software program.Using the ANOVA statistical test for the categorical covariate and the Pearson's correlation test for the continuous covariate, a p-value can be calculated.Regardless of their incorporation into the model, the random effect-covariate associations are sorted using the p-values.The forward and backward method was employed to select covariates.Until there are no correlation p-values over a threshold, the covariate with the smallest correlation p-value is included in the model, or the next smallest if the smallest has already been attempted.Until there are no correlation p-values below a threshold, the covariate with the highest correlation p-value is disregarded, or the next highest if the highest value has already been attempted.Therefore, covariates with a p-value less than 0.05 that improved the fit while lowering BIC remained in the model.
For the evaluation of the proposed final model, we employed statistical tests using individual parameters drawn from the conditional distribution.Correlation and Wald test were used to assess whether covariates should be removed from the model.The normal distribution of random effects and the symmetry around 0 of residuals' distribution were examined through Shapiro-Wilk test and symmetry test, respectively.

External Validation
External validation was conducted using data from an open-label, randomized, crossover, two-cohort, single-dose study in healthy volunteers to evaluate the unit dose dry powder inhaler (UD-DPI) for the delivery of salbutamol and to compare the PK profile with the MDI and Diskus presentations (NCT01984086) [27], with access granted by Glaxo-SmithKline (GSK).The focus of data extraction was on Part A of the study, encompassing 30 individuals subjected to Treatment A. This treatment involved a single administration of salbutamol (as sulfate) (200 µg per blister of 1.6% blend), delivered via UD-DPI through the inhalation of 3 blisters, resulting in a cumulative dose of 600 µg.The selection of this treatment aimed to align with the dose administered in the development of the popPK model using virtual patients.

of 18
The study cohort comprised healthy nonsmoking male/female individuals, aged 18 to 65 years, with a body weight ≥ 50 kg, and a body mass index (BMI) within the range of 19.0 to 34.0 kg/m 2 .Participants had to provide written informed consent.Main exclusion criteria included current or history of chronic liver disease, a history of salbutamol sensitivity, positive tests for Hepatitis B, Hepatitis C, and HIV, pregnancy, breastfeeding, or active attempts to conceive.Further details can be found in the study protocol.
Upon acquisition of the dataset for input into Monolix software, the popPK model derived from virtual data was applied, following the previously outlined steps.The validation of the virtual data-based model was performed through visual inspection of outcomes after applying the same popPK model to the clinical dataset.A general scheme of the workflow is shown in Figure 1.
30 individuals subjected to Treatment A. This treatment involved a single administration of salbutamol (as sulfate) (200 µg per blister of 1.6% blend), delivered via UD-DPI through the inhalation of 3 blisters, resulting in a cumulative dose of 600 µg.The selection of this treatment aimed to align with the dose administered in the development of the popPK model using virtual patients.
The study cohort comprised healthy nonsmoking male/female individuals, aged 18 to 65 years, with a body weight ≥ 50 kg, and a body mass index (BMI) within the range of 19.0 to 34.0 kg/m 2 .Participants had to provide written informed consent.Main exclusion criteria included current or history of chronic liver disease, a history of salbutamol sensitivity, positive tests for Hepatitis B, Hepatitis C, and HIV, pregnancy, breastfeeding, or active attempts to conceive.Further details can be found in the study protocol.
Upon acquisition of the dataset for input into Monolix software, the popPK model derived from virtual data was applied, following the previously outlined steps.The validation of the virtual data-based model was performed through visual inspection of outcomes after applying the same popPK model to the clinical dataset.A general scheme of the workflow is shown in Figure 1.

Results and Discussion
During the process of developing a precise popPK model, pharmacometricians encounter significant challenges, particularly regarding the substantial quantity and access to PK data and personal information required, due to several limitations imposed by company policies and ethical considerations.In addition, data derived from clinical trials are often incomplete, especially datasets collected in large, late-phase trials or during routine

Results and Discussion
During the process of developing a precise popPK model, pharmacometricians encounter significant challenges, particularly regarding the substantial quantity and access to PK data and personal information required, due to several limitations imposed by company policies and ethical considerations.In addition, data derived from clinical trials are often incomplete, especially datasets collected in large, late-phase trials or during routine healthcare or follow-up visits.Most of these data may be missing or inaccurate due to a variety of factors such as study site noncompliance, patient noncompliance, inappropriate sample handling, data entry errors, and analytical issues.Consequently, how these missing or erroneous data are handled can impact their interpretation.Although popPK modeling requires a smaller volume of data points, that is, it is possible to develop an accurate model with sparse data (few observations per subject) [30], sometimes even this is difficult to obtain in clinical studies (especially in Phase I trials), because volunteers have to undergo a number of blood draws, raising serious ethical questions, especially in special populations such as children.Indeed, there are some guidelines that help pharmacometricians in dealing with these challenging data, but the truth is that this will always have a significant impact on the final model and outcomes [32][33][34].
In this context, PK data for salbutamol are scarce, as we previously stated [35], and to circumvent all the aforementioned issues, our study proposes a new way to study therapeutic regimen optimization through the generation of a synthetic dataset, which encompasses a diverse set of virtual patients.This approach has already been employed in recent studies and involves the artificial generation of data that mimics real-world data [36][37][38][39].For example, artificial intelligence (AI) algorithms have been proposed in this scope.In the context of our study, a virtual dataset was generated using PBPK models to simulate salbutamol PK profiles across various demographic characteristics.

Demographic Characteristics
Thirty-two virtual subjects were included in this study.Among these subjects, there are 18 males and 14 females, with a median age of 20.0 years and a BMI median of 21.3 kg/m 2 .Furthermore, the virtual population is composed of two groups, according to race: American Indian or Alaskan Native represent 38% of total patients, with East Asian (Japanese and Chinese) patients representing 62% of the sample.These demographic characteristics are summarized in Table 1.ND-not determined.
Notable differences across various parameters are observed between virtual patients and patients enrolled in the clinical study.The mean age is 26.8 years.Gender distribution varies, with virtual patients having a more balanced ratio, while clinical study patients skew towards a higher percentage of males (73%).In terms of body composition, virtual patients, on average, have a lower BMI, height, and weight.This is attributed to the inclusion of Japanese and Chinese individuals in the virtual dataset, who inherently present lower values.Consequently, racial composition significantly differs between the two datasets.Virtual patients are predominantly East Asian or American, whereas patients included in the clinical study are mostly Caucasian.In fact, studies of this nature do not accurately represent the "real-world" population [40].Clinical trials often impose stringent criteria to ensure sample homogeneity and minimize confounding variables, employing limited analyses tailored to address specific research questions.Further, voluntary participation introduces a potential bias, as participants are self-selected and willing to engage.Additionally, due to resource constraints, many studies exhibit limited representation of racial and ethnic groups, resulting in underrepresentation or inadequate portrayal of population diversity, creating a gap between the trial world and the real world [41][42][43][44].
Efforts within the scientific community have been undertaken to address this gap by enhancing diversity in clinical trial recruitment [40].Therefore, this virtual dataset aims to present a more diverse ethnic group.Ideally, it would include Caucasians as well; however, the PBPK modeling software mimicking the patient cohort has limitations concerning ethnic groups.
Subsequently, an NCA of both datasets was conducted, allowing the collection of PK parameters to compare their distribution between virtual and real populations.Table 2 displays the C max , AUC, T max , and the volume of plasma from which drug is removed from the human body (Cl).The observed differences in values may be attributable to variations in the characteristics of the datasets.The exploration of virtual patient data is depicted in Figure 2. The data were stratified based on some crucial individual characteristics: age, weight, race, and gender.A noticeable difference is observed between the younger and the adult and elderly populations-younger subjects exhibit higher concentrations of salbutamol over time.The pediatric population is recognized as a special group for drug therapy, as many physiological changes occurring at these ages have a significant impact on the PK of several compounds.Factors such as pH, differences in bile luminal concentrations, intestinal permeability, and variations in enzymatic activity may explain higher plasma concentrations of salbutamol in younger individuals [45].Consequently, toxicity may occur in this age range, highlighting the importance of dose adjustment in children.In the context of this study, the dose administered to virtual patients under 18 years is equivalent to the dose for children (up to 100-200 µg as a single dose on demand) [46].
Weight also influences the PK profile of this drug, with individuals with more than 75 kg (overweight and obese) displaying reduced values.Regarding race and gender, with this preliminary analysis, it is not visually discernible to determine the impact that different groups may have on the kinetics of this β-agonist.
Data exploration also uncovered several correlations between PK parameters and covariates.There is a negative correlation between the AUC and both height and weight, implying that individuals with greater height or weight tend to have lower overall exposure to the drug over time.Similarly, a negative correlation was observed between C max and height/weight.On the other hand, there is a positive correlation between drug clearance and height or weight.Individuals with higher height or weight tend to eliminate the drug from their bodies at a faster rate.
An interesting observation is that Asians exhibit higher AUC and C max values and lower clearance values compared to the American population, suggesting that, on average, Asians may experience higher overall drug exposure and slower elimination of the drug from their bodies.
An interesting observation is that Asians exhibit higher AUC and Cmax values and lower clearance values compared to the American population, suggesting that, on average, Asians may experience higher overall drug exposure and slower elimination of the drug from their bodies.

Final popPK Model
The plasma concentration-time profiles of inhaled salbutamol in DPI form were best described by a two-compartment model, with first-order absorption (no lag time), and linear elimination.During the development of the structural model, several models were created and evaluated based on criteria detailed in Section 2.3.The model with the lowest BICc value, following visual inspection of the VPC, was selected.BICc values for each tested model are provided in the Supplementary Material (Table S2).In particular, for this model, a −2LL of −6480.18 and a BICc of −6420.66 were obtained.
Therefore, the most suitable model for this dataset included first-order absorption, as salbutamol is absorbed through passive diffusion [47], a common mechanism in drugs following first-order kinetics, without lag time, considering salbutamol is inhaled and has an immediate onset of action at its target.Two compartments were also incorporated, a central compartment and a peripheral compartment, although the literature reports a onecompartment model to describe salbutamol PK.In fact, after inhalation of salbutamol, a

Final popPK Model
The plasma concentration-time profiles of inhaled salbutamol in DPI form were best described by a two-compartment model, with first-order absorption (no lag time), and linear elimination.During the development of the structural model, several models were created and evaluated based on criteria detailed in Section 2.3.The model with the lowest BICc value, following visual inspection of the VPC, was selected.BICc values for each tested model are provided in the Supplementary Material (Table S2).In particular, for this model, a −2LL of −6480.18 and a BICc of −6420.66 were obtained.
Therefore, the most suitable model for this dataset included first-order absorption, as salbutamol is absorbed through passive diffusion [47], a common mechanism in drugs following first-order kinetics, without lag time, considering salbutamol is inhaled and has an immediate onset of action at its target.Two compartments were also incorporated, a central compartment and a peripheral compartment, although the literature reports a one-compartment model to describe salbutamol PK.In fact, after inhalation of salbutamol, a significant portion of the dose is swallowed, and only a small fraction (around 20%) reaches the specific target (lungs) [48].Thus, the inclusion of two compartments in the PK model may be logical, with the central compartment being the lungs, where drug distribution occurs, binding to the β 2 -adrenergic receptors present in lung tissue, and the peripheral compartment being the gut, where the drug is delivered due to swallowing that occurs in inhaled drugs.Additionally, the most appropriate elimination type is linear elimination, which means that the elimination rate is directly proportional to the drug concentration.
Covariate analysis, employing the model proposed by Monolix ® , encompassed gender, age, height, weight, and racial population.The incorporation of all covariates significantly improved the base structural model.Therefore, we highlighted the influence of individual patient characteristics on salbutamol's kinetics: age exhibited a significant impact on Cl and intercompartmental clearance (Q); gender influences absorption constant rate (k a ) and Cl; race has an impact on Cl, and weight showed a considerable influence on Cl, Q, and on the volume of distribution of peripheral compartment (V2).The interindividual variability (IIV) in these PK parameters was considered by applying a normal distribution and a combined residual model, consisting of proportional and additive terms.Parameter estimates for the final structural model are shown in Table 3.All structure parameters were estimated with good precision, with relative standard error (RSE) < 30% for fixed effects, except for the volume of distribution of the central compartment (V1), which exhibited an infinitely large standard error.Regarding the RSE of random effects, there is insufficient accuracy in some estimates.A correlation between k a and Q was identified from scatter plots for each pair of random effects, along with the Pearson correlation coefficient.This correlation was introduced into the model to achieve a more accurate estimation of these parameters.

Model Assessment Using the Diagnostic Plots
The low p-values (p < 0.05) in Pearson's correlation test confirmed that the correlation between parameters and covariates was significant.These tests were in agreement with Wald test results, where low p-values (p < 0.05) indicated the relevance of covariate effects.
The goodness-of-fit plots of the final PK model with covariates are shown in Figure 3.The individual predicted values, obtained through the conditional model (empirical Bayes estimates, EBEs) demonstrated strong agreement with the observed salbutamol plasma concentrations.The proportion of outliers, determined by assessing the balance of points on each side of the y = x line and the points lying outside the 90% prediction interval, was 4.62%, suggesting an insignificant proportion.Both NPDE and IWRES plots revealed that the residuals are randomly scattered around the horizontal zero-line, indicating that residuals behave as independent standardized normal random variables.Figure 4 illustrates the VPC plot of the final model.It reveals a close alignment between observed (empirical) and predicted (theoretical), within their respective prediction intervals, along with a reduced proportion of outliers.This alignment suggests an adequate fit of the model to the observed data.These diagnostic plots allow for a clearer detection of potential misspecifications in the model.Overall, the developed model is well suited to the salbutamol PK data generated using virtual patients.
terval, was 4.62%, suggesting an insignificant proportion.Both NPDE and IWRES plots revealed that the residuals are randomly scattered around the horizontal zero-line, indicating that residuals behave as independent standardized normal random variables.Figure 4 illustrates the VPC plot of the final model.It reveals a close alignment between observed (empirical) and predicted (theoretical), within their respective prediction intervals, along with a reduced proportion of outliers.This alignment suggests an adequate fit of the model to the observed data.These diagnostic plots allow for a clearer detection of potential misspecifications in the model.Overall, the developed model is well suited to the salbutamol PK data generated using virtual patients.

External Validation of the Final Structural popPK Model
External validation was conducted by applying the structural model obtained previously to the GSK-derived dataset.The variation of plasma concentration profile over time in the patients enrolled in the clinical study can be visualized in Figure 4, according to stratifications by age, weight, race, and gender.No significant differences are observed, as was noted in the virtual dataset, likely due to the homogeneity of the population.Indeed, this remains a challenge in clinical trials, which has been discussed over the years.Clinical trials are the cornerstone of evidence-based medicine; however, the quality of evidence provided is often poor, particularly regarding participant recruitment.Many clinical researchers tend to perceive heterogeneity as a burden that must be eliminated or controlled through refined statistical approaches to reduce the effects of interindividual variability and achieve consistent results [41][42][43][44].Nevertheless, the population does not accurately represent the diversity of the overall population.In a 2020 analysis of global participation in clinical trials, the FDA revealed significant disparities among various racial groups: 76% of participants were white, 11% were Asian, and only 7% were black [49].With the shift in healthcare towards precision medicine, variability in responses across different subpopulations has been recognized.For instance, data from the clinical trial sponsored by GSK predominantly include Caucasian patients, as previously discussed.Several reasons may explain this occurrence.Inclusion and exclusion criteria are usually limited due to ethical and scientific considerations.Including children, for example, in general clinical trials is not very common unless they are pediatric studies designed for that purpose.Participant recruitment alone significantly limits the study, making recruitment of a diverse population in terms of characteristics even more challenging.Attempts to minimize potential confounding variables are also one of the reasons behind the homogeneity of these studies.Ethical and policy issues involving stakeholders may also explain this.Thus, enrolling a homogenous population is not beneficial in any instance for studies seeking to understand the impact of interindividual variability.

External Validation of the Final Structural popPK Model
External validation was conducted by applying the structural model obtained previously to the GSK-derived dataset.The variation of plasma concentration profile over time in the patients enrolled in the clinical study can be visualized in Figure 4, according to stratifications by age, weight, race, and gender.No significant differences are observed, as was noted in the virtual dataset, likely due to the homogeneity of the population.Indeed, this remains a challenge in clinical trials, which has been discussed over the years.Clinical trials are the cornerstone of evidence-based medicine; however, the quality of evidence provided is often poor, particularly regarding participant recruitment.Many clinical researchers tend to perceive heterogeneity as a burden that must be eliminated or controlled through refined statistical approaches to reduce the effects of interindividual variability and achieve consistent results [41][42][43][44].Nevertheless, the population does not accurately represent the diversity of the overall population.In a 2020 analysis of global participation in clinical trials, the FDA revealed significant disparities among various racial groups: 76% of participants were white, 11% were Asian, and only 7% were black [49].With the shift in healthcare towards precision medicine, variability in responses across different subpopulations has been recognized.For instance, data from the clinical trial sponsored by GSK predominantly include Caucasian patients, as previously discussed.Several reasons may explain this occurrence.Inclusion and exclusion criteria are usually limited due to ethical and scientific considerations.Including children, for example, in general clinical trials is not very common unless they are pediatric studies designed for that purpose.Participant recruitment alone significantly limits the study, making recruitment of a diverse population in terms of characteristics even more challenging.Attempts to minimize potential confounding variables are also one of the reasons behind the homogeneity of these studies.Ethical and policy issues involving stakeholders may also explain this.Thus, enrolling a homogenous population is not beneficial in any instance for studies seeking to understand the impact of interindividual variability.
The two-compartment model with first-order absorption, no delay, and linear elimination was found to fit the data when analyzing diagnostic plots (Figures 5-7).Despite an outlier proportion of 7.62% (also visually confirmed in the VPC plot), the model is considered suitable for the data.A −2LL of −7138.63 and a BICc of −7058.17 were obtained.The estimates of the popPK parameters are outlined in Table 4.The two-compartment model with first-order absorption, no delay, and linear elimination was found to fit the data when analyzing diagnostic plots (Figures 5-7).Despite an outlier proportion of 7.62% (also visually confirmed in the VPC plot), the model is considered suitable for the data.A −2LL of −7138.63 and a BICc of −7058.17 were obtained.The estimates of the popPK parameters are outlined in Table 4.Despite the identified demographic differences between the virtual patients and the clinical trial participants, the range of values in the two datasets is comparable, allowing the validation of the popPK model developed using virtual patients with real-world data.In addition, the need for more heterogeneous clinical trials implies more extensive and costly studies [41].Therefore, by validating the model developed from synthetic data with real-world data, we can employ this methodology in therapeutic regimen optimization studies.Its advantages are evident, including flexibility, cost-effectiveness, and the ability to simulate diverse patient populations.

Impact of Covariates on Salbutamol PK Parameters
The developed model revealed the influence of covariates on the PK parameters of salbutamol.Thus, we proceeded to analyze these parameters considering various stratifications for age, weight, gender, and race (Table 5).Firstly, the model identified that age influenced Cl and Q, which, as observed, presented significantly different values when comparing the two subgroups of asthmatic patients (5-22 years and 23-65 years).Q is a parameter that specifies the drug transfer rate between compartments and mainly impacts the drug volume in compartments, which consequently differs significantly between the two subgroups.In general, the V2 is much higher than the volume in V1, suggesting a lower drug availability to target its site of action, increased risk of ADRs due to accumulation in peripheral tissues, or prolonged drug action.At this point, it is not possible to distinguish whether it is a positive or negative effect.In contrast, Cl is higher in the younger population (5-22 years), potentially indicating a need for closer observation in younger individuals.Maximum efficacy may not be reached, and a dose adjustment may be required.Regarding the absorption of this β 2 -agonist, younger people display increased k a , following the PK profile observed in Figure 2A.Indeed, several underlying processes related to aging influence drug disposition, such as reduced first-pass metabolism, decreased organ/tissue mass, increased body fat, and decreased body water content [50].Concerning the latter one, salbutamol being a hydrophilic drug, with a greater affinity for water, tends to have a lower distribution volume, as confirmed by the presented results.
Weight, in turn, proved to be an important characteristic in the disposition of salbutamol.For individuals weighing over 75 kg (overweight and obese), markedly reduced values were observed for all PK parameters.Regarding the popPK model, relationships between this covariate and Cl, Q, and V2 were identified.Variations in Cl may be explained by several factors.For instance, fat accumulation in the liver (the main organ responsible for elimination) may alter hepatic blood flow [51,52].Changes in enzyme expression of cytochrome P450 (CYP) enzymes, especially CYP2D6 and CYP2C19, may also underlie these decreased values [53].Further, the discrepancy in drug distribution values is reasonable, considering the hydrophilicity of salbutamol-individuals with higher weight have more adipose tissue.Tissue perfusion may also be reduced in obese patients, impacting Vd (volume of distribution) [52].The drug absorption rate was also much lower in individuals > 75 kg compared to individuals with normal weight (<75 kg, according to software guidelines).This difference can be attributed once more to changes in blood flow and body composition.
According to the initial population analysis, gender influences k a and Cl, which is consistent with the results presented in Table 5.K a and Cl are both reduced in females.However, Cl is, indeed, decreased in women (approximately half of the Cl value displayed in men).These variations may be explained by anatomical differences, namely, weight and total body water volume [4,54].Consequently, alterations in drug distribution are also expected.
Finally, significant differences were found in Vd and k a between the two racial groups, in contrast to the popPK model results.In particular, American people exhibit higher values for k a and Vd.Among many other factors, well-recognized polymorphisms between racial groups may explain the variation in these PK parameters [55].Several single-nucleotide polymorphisms (SNPs) have been identified in the coding region of the gene, the gene encoding for the β 2 -adrenoceptor.The most common SNPs result from three missense mutations, one of which involves the substitution of isoleucine (Ile) for threonine (Thr) at codon 164 [56].According to in vitro studies conducted by Chung et al. [57], this alteration has been associated with decreased receptor binding to the ligand and coupling to Gs proteins in response to different SABAs, namely, salbutamol.This results in a reduction in receptor activity and agonist-stimulated activation.In this regard, the ethnic disparities reflected by the polymorphisms significantly influence the PK parameters of salbutamol, and it is crucial to take them into account when prescribing the drug.

Conclusions
This study addresses critical challenges in popPK modeling, particularly regarding data scarcity, incompleteness, and homogeneity in traditional clinical trials.Although the landscape of clinical trials has been changing in recent years, pharmacometricians still face a burden when it comes to PK data.Therefore, by leveraging synthetic data generated through PBPK models, we have overcome these limitations and provided insights into salbutamol's PK profile across diverse patient populations.External validation using real data from clinical trials is crucial for ensuring the accuracy and reliability of the final models.
Through this comprehensive analysis, significant associations between individual characteristics and salbutamol's PK parameters have been identified.Age, weight, gender, and race have indeed a great impact on salbutamol's ADME processes.This underscores the relevance of personalized dosing strategies, minimizing adverse effects and maximizing therapeutic efficacy.Although the popPK model fails to identify some relationships, these findings revealed great consistency.Further studies using this methodological approach should be conducted, to prove its reliability and capability of predicting interindividual variability in drug responses.

Figure 2 .
Figure 2. Observed data of plasma salbutamol concentration (µg/mL) through time (h), stratified by (A) age, (B) weight, (C) race, and (D) gender.The absorption phase cannot be visualized.

Figure 2 .
Figure 2. Observed data of plasma salbutamol concentration (µg/mL) through time (h), stratified by (A) age, (B) weight, (C) race, and (D) gender.The absorption phase cannot be visualized.

Figure 3 .
Figure 3. Goodness-of-fit plots of the final popPK model: (A) Observed salbutamol concentration vs. individual predictions using conditional mode (EBEs); and (B) IWRES vs. time and individual predicted, and NPDE vs. time and population predicted.

Figure 3 .
Figure 3. Goodness-of-fit plots of the final popPK model: (A) Observed salbutamol concentration vs. individual predictions using conditional mode (EBEs); and (B) IWRES vs. time and individual predicted, and NPDE vs. time and population predicted.

Figure 4 .
Figure 4. VPC for the two-compartment model with first-order absorption, no delay, and linear elimination.The VPC of the final covariate PK model displays the prediction intervals for the 10th, 50th, and 90th percentiles (shaded areas, from bottom to top, respectively).Outliers are visualized as red dots and areas.

Figure 4 .
Figure 4. VPC for the two-compartment model with first-order absorption, no delay, and linear elimination.The VPC of the final covariate PK model displays the prediction intervals for the 10th, 50th, and 90th percentiles (shaded areas, from bottom to top, respectively).Outliers are visualized as red dots and areas.

Figure 6 .
Figure 6.Observed salbutamol concentration vs. individual predictions using conditional mode (EBEs), for the two-compartment model with first-order absorption, no delay, and linear elimination applied to the clinical study dataset for external validation.The dotted line represents 90% prediction interval, the black line refers to y = x line, and the yellow line represents the spline.

Figure 6 .
Figure 6.Observed salbutamol concentration vs. individual predictions using conditional mode (EBEs), for the tw compartment model with first-order absorption, no delay, and linear elimination applied to the clinical study dataset for external validation.The dotted line represents 90% prediction interval, the black line refers to y = line, and the yellow line represents the spline.

Figure 7 .
Figure 7. VPC for the two-compartment model with first-order absorption, no delay, and linear elimination applied the clinical study dataset for external validation.

Figure 7 .
Figure 7. VPC for the two-compartment model with first-order absorption, no delay, and linear elimination applied to the clinical study dataset for external validation.

Table 1 .
Subject characteristics (mean or median ± standard deviation, SD or interquartile range, IR) from virtual dataset and clinical trial.

Table 2 .
PK parameters of the virtual and clinical datasets using NCA.
C max -maximum observed concentration; AUC-area under the curve from the time of dosing to the last measurable positive concentration; T max -time of maximum observed concentration; Cl-clearance; NA-not applicable.

Table 3 .
Estimates of the population pharmacokinetic parameters of the final model for salbutamol DPI formulation.
RSE: relative standard error; k a : absorption constant rate; Cl: clearance; V1 and V2: the volume of distribution of the compartments one (central), and two (peripheral); Q: intercompartmental clearance; NaN: infinitely large standard error.

Table 4 .
Estimates of the popPK parameters of the final model for salbutamol DPI formulation (using GSK-deriv data).

Table 4 .
Estimates of the popPK parameters of the final model for salbutamol DPI formulation (using GSK-derived data).

Table 5 .
Differences in PK parameters according to the stratification of the virtual population in terms of age, weight, gender, and race.Values presented in geometric mean ± geometric SD.
a : absorption constant rate; Cl: clearance; V1 and V2: the volume of distribution of the compartments one (central), and two (peripheral); Q: intercompartmental clearance.