Li Shen1,2†
Mengdi Hu1†Xiaoyong Xu1
Yuxuan Zhou3
Wei Wu1,4,5
Xilin Ge1Guangfei Wang1
Yi Wang6*
Zhiping Li1,7*- 1Department of Clinical Pharmacy, Children’s Hospital of Fudan University, National Children’s Medical Center, Shanghai, China
- 2Department of Pharmacy, Suzhou Hospital, Affiliated Hospital of Medical School, Nanjing University, Suzhou, Jiangsu, China
- 3Global Health Research Center, Duke Kunshan University, Kunshan, Jiangsu, China
- 4Department of Pharmacy, The First Affiliated Hospital of Fujian Medical University, Fuzhou, China
- 5Department of Pharmacy, National Regional Medical Center, Binhai Campus of the First Affiliated Hospital, Fujian Medical University, Fuzhou, China
- 6Department of Neurology, Children’s Hospital of Fudan University, National Children’s Medical Center, Shanghai, China
- 7Department of Pharmacy, Kunshan Women and Children’s Healthcare Hospital, Children’s Hospital of Fudan University, Kunshan Branch, Kunshan, Jiangsu, China
Objective: This study aimed to develop an individualized dosing strategy for voriconazole (VRZ) in children under 2 years of age by integrating machine learning (ML) and population pharmacokinetic (PopPK) modeling.
Methods: This retrospective observational study included 76 eligible pediatric patients for model development, analyzing their baseline characteristics and laboratory parameters. A population pharmacokinetic (PopPK) model using NONMEM® software was performed to assess the clearance (CL) and volume of distribution (V) of VRZ. The individual CL and V were included as input variables. The Boruta algorithm was employed for feature selection, after which six machine learning algorithms were applied. The models were evaluated using Mean Squared Error (MSE), Root Mean Squared Error (RMSE), Mean Absolute Error (MAE), and coefficient of determination (R2) to identify the optimal algorithm, which then underwent independent external validation. The selected final model was analyzed for interpretability using Shapley Additive Explanations (SHAP).
Results: A total of 76 pediatric patients were enrolled for model development, consisting of 58 males (76.3%) and 18 females (23.7%), with a median age of 11 months and a median weight of 8.05 kg. We analyzed 110 therapeutic drug monitoring (TDM) samples of VRZ from these participants. A one-compartment model with first-order absorption and elimination described the population pharmacokinetics of VRZ. Population estimates for apparent clearance (CL/F) and volume of distribution (V/F) were 17.9 L/h/70kg (RSE, 10.8%) and 788 L/70kg (RSE, 15.4%), respectively. An XGBoost model accurately predicted voriconazole concentrations (R2 = 0.81, RMSE = 0.53) with a relative error of ±20% for most observations. In the external validation, the XGBoost model demonstrated an R2 of 0.75, RMSE of 0.14. SHAP analysis identified clearance, weight, and laboratory values as significant predictors.
Conclusion: This study emphasized the importance of personalized treatment in utilizing VRZ for children under 24 months. The XGBoost model demonstrated potential in identifying an initial dose recommendation for VRZ.
1 Introduction
Voriconazole (VRZ), a triazole antifungal, is employed extensively in the management of invasive fungal infections among pediatric patients (Kadam and Van Den Anker, 2016), particularly those with weakened immune systems due to conditions such as hematologic malignancies or after undergoing hematopoietic stem cell transplantation (HSCT). Although approved by the US FDA and European Medicines Agency (EMA) for children aged two and older, a scarcity of comprehensive pharmacokinetic data complicates its use in infants under two (Chen et al., 2018). This data gap requires that clinicians carefully weigh the potential therapeutic benefits for an individual patient against the known pharmacokinetic variability, leading them to frequently resort to its “off-label” use after a careful risk-benefit assessment. In such cases, therapeutic drug monitoring (TDM) is essential for addressing the high pharmacokinetic variability in this young population and for optimizing both safety and efficacy (Touw and van den Anker, 2022). Pediatric patients, especially those below 24 months, exhibit different pharmacokinetic profiles compared to older children and adults, necessitating careful consideration of dosing strategies (Yan et al., 2018; Gastine et al., 2023). Although effective, VRZ’s pharmacokinetics (PK) in children showed considerable variability, frequently influenced by age, body weight, and genetic polymorphisms (Soler-Palacin et al., 2012). This challenge is magnified in children under 2 years old, where profound PK variability - driven by developmental changes and genetic polymorphisms - makes optimal dosing difficult.
The pharmacokinetics of VRZ may vary significantly among patients, especially in critically ill individuals, pediatrics, the elderly, and the obese (Wang et al., 2025; Jiang et al., 2022). Understanding VRZ population pharmacokinetics is essential for optimal dosing, given its narrow therapeutic index and patient variability from factors like C-reactive protein levels and genetic polymorphisms (van den Born et al., 2023; Li et al., 2023). Population pharmacokinetic analysis is a powerful approach that effectively addresses dosing challenges by understanding voriconazole behavior and providing insights for optimizing clinical dosing strategies, even with limited sampling data. However, classical PopPK methods rely on fixed compartmental models that may limit the ability to fully capture all factors influencing drug behavior, particularly when working with sparse sampling data (Liu et al., 2025).
ML, a core branch of artificial intelligence, is adept at capturing complex patterns and nonlinear relationships within big data (Asnicar et al., 2024; Greener et al., 2022). ML techniques are emerging as valuable complements to PopPK approaches in optimizing voriconazole dosing. For instance, a recent study (Cheng et al., 2023) demonstrated that machine learning algorithms, particularly XGBoost, could predict toxic plasma trough concentrations of VRZ (>5 mg/L) with an accuracy of 78.8%, identifying key factors such as albumin and total bile acid for dosage optimization. Utilizing pharmacokinetic (PK) parameters in machine learning models is an innovative strategy that has been proven effective with various drugs (Li et al., 2024; Ma et al., 2024; Huang et al., 2022), leading to improved model predictive accuracy and enhanced model interpretability.
As VRZ is a crucial option for managing invasive fungal diseases, and considering the limited data on dosing and exposure in children under 24 months, this study aims to share experiences with individualized voriconazole therapy. We developed a model that integrate PopPK with ML approaches to predict the steady-state plasma concentration of oral VRZ in this age group. The study workflow was presented in Figure 1.
Figure 1. Study process workflow.
2 Methods
2.1 Study design
We conducted a retrospective observational study analyzing data from pediatric patients aged under 24 months who had been hospitalized at the Children’s Hospital of Fudan University from January 2020 to June 2025. Inclusion criteria included patients aged below 2 years who received oral VRZ for a minimum of 3 days during their hospitalization and underwent therapeutic drug monitoring (TDM) (Chen et al., 2018; Resztak et al., 2021). Exclusion criteria included newborns, preterm infants, patients requiring advanced life support therapies, and patients whose medical information was unavailable or missing more than 30%. For independent external validation, we also collected data from 10 patients under 24 months who received oral VRZ and underwent TDM at our hospital between July and August 2025.The ethics committee of Children’s Hospital of Fudan University approved the study protocol [No. (2024) 321]. Informed consent was waived for the collection and analysis of the anonymous data, as there was no intervention involved, and the study posed minimal risk (Filion et al., 2016; Strutz et al., 2025).
For patients who satisfied the enrollment criteria, a comprehensive set of data was collected from their electronic medical records. Clinical characteristics included age, sex, height (HT), and weight (WT), from which the body surface area (BSA) was calculated using the formula: BSA (m2) =
Treating physicians determined VRZ doses individually for each patient, considering clinical status and therapeutic drug monitoring. While referencing the standard 9 mg/kg q12h regimen for older children (2 to <12 years), oral doses for children under 2 years ranged from 4 to 10 mg/kg every 12 h. This dosing strategy is consistent with recommendations from UpToDate, a clinical decision support system (Isaac et al., 2012). Blood samples were collected to determine steady-state trough concentrations, with sampling performed within the 30-min window immediately preceding the next scheduled dose. In cases of missing covariate data, which accounted for less than 3% of the dataset, median imputation was employed (Shen et al., 2024).
2.2 VRZ therapeutic drug monitoring
The plasma concentration of VRZ was determined by high-performance liquid chromatography (HPLC) using a Shimadzu LC-20 series HPLC system (Shimadzu, Japan). A Waters XBridge C18 chromatographic column (4.6 × 150mm, 5 μm) was used. The mobile phase A consisted of water, while mobile phase B was acetonitrile. The injection volume was 20 μL, with the detection wavelength fixed at 256 nm. The standard curve for VRZ covered a linear range of 0.25–10 mg/L. The target concentration for VRZ was set between one and 5.5 mg/L in plasma trough (Pappas et al., 2016).
2.3 Population pharmacokinetic modeling
The PopPK analysis of VRZ was carried out in NONMEM® software (version 7.5.0, ICON plc., USA) executing the first-order conditional estimation method with interaction (FOCE-I). The pharmacokinetic profile of VRZ was characterized using the ADVAN2 TRANS2 subroutine to develop the basic structural model, representing a one-compartment disposition with first-order absorption and elimination. The pharmacokinetic parameters assessed were the apparent clearance (CL/F), the apparent volume of distribution (V/F), and the absorption rate constant (Ka). Due to the sampling design primarily focusing on steady-state trough concentrations, individualized absorption parameters could not be reliably assessed. Ka was fixed to 1.19 h-1, reflecting published literature values (Gastine et al., 2018). An exponential model was used to quantify interindividual variability in CL/F and V/F. Various residual error models (proportional, additive or combined) were evaluated to assess intra-individual variability. Using allometric scaling, we examined the effects of body weight on PK parameters.
Covariates were evaluated utilizing a stepwise forward-addition/backward-elimination approach based on objective function value (OFV) changes. During forward selection, covariates were included if they reduced the OFV by > 3.84 (p < 0.05). In the backward elimination phase, covariates were retained in the final model only if their removal increased the OFV by > 6.63 (p < 0.01). Final model evaluation was conducted through diagnostic goodness-of-fit (GOF) plots, bootstrap analysis, and visual predictive check (VPC). Parameter precision was evaluated, and 95% confidence intervals (95%CI) were generated using bootstrap resampling (n = 1000). VPC analysis was performed using 1000 simulated datasets to characterize the distribution characteristics of observed versus simulated data. The final PopPK model enabled the calculation of pediatric individual PK parameters CL and V via the empirical Bayesian estimates method. These PK parameters were subsequently introduced as features in the ML model.
2.4 Machine learning modeling
For model development and validation, the dataset was randomly partitioned into a training set (70%) and a validation set (30%). The Boruta algorithm was employed for feature selection to identify all relevant features (Shen et al., 2024). The Boruta algorithm automates feature selection by evaluating the significance of original features against shadow features (Li et al., 2025). The Boruta was executed for 200 iterations with a significance threshold of p < 0.01, applying Bonferroni correction for multiple comparisons to identify statistically significant features. Six ML algorithms were chosen to construct models: extreme gradient boosting (XGBoost), light gradient boosting machine (LightGBM), gradient boosting decision tree (GBDT), categorical boosting (CatBoost), adaptive boosting (AdaBoost), and random forest (RF). To evaluate the stability and generalizability of our final model, we performed a 10-iteration bootstrap cross-validation and an independent external validation.
Model performance was assessed using four standard regression metrics: Mean Squared Error (MSE), Root Mean Squared Error (RMSE), Mean Absolute Error (MAE), and coefficient of determination (R2), providing a comprehensive evaluation of predictive accuracy. Lower MSE, RMSE, and MAE values signify enhanced model performance, whereas R2 values approaching one indicate greater explanatory power. After establishing the optimal model, we used relative accuracy (RA) as an indicator to evaluate the model’s predictive performance (Huang et al., 2020). In this study, RA was measured by whether the relative error between VRZ predicted values and observed values falls within a specified percentage threshold of 20%. The optimal model was chosen for interpretability analysis using the Shapley Additive Explanations (SHAP) method. The approach assessed each feature’s impact on specific predictions, offering insights into both local and global model behavior (Qi et al., 2025).
2.5 Statistical analysis
Categorical variables were presented using frequencies and percentages [N (%)]. Continuous variables were described using both medians with interquartile ranges [Median (IQR)] and means with standard deviations [Mean ± SD]. The Kruskal–Wallis test was used to analyze continuous variables, while Fisher’s exact test was employed for categorical variables. A p-value of below 0.05 was established as the threshold for statistical significance. VRZ plasma concentrations below the lower limit of quantitation (LLOQ) of 0.25 mg/L were labeled as below quantitation limit (BQL). These BQL values were substituted with half the LLOQ value (0.125 mg/L) (Chabala et al., 2022). In this study, BQL data accounted for less than 4% of the total observations. R software (v 4.5.0) and Python (v 3.10.6) were utilized for all analyses.
3 Results
3.1 Clinical characteristics of pediatric patients
The clinical and demographic characteristics of the 76 pediatric patients, from whom 110 samples were collected, were summarized in Table 1. The study population consisted of 58 males (76.3%) and 18 females (23.7%). The median age, height, and weight were 11.0 months (IQR, 7.38–17.00), 70.00 cm (IQR, 66.00–74.00), and 8.05 kg (IQR, 6.95–9.00), respectively. The median voriconazole plasma concentration was 1.25 mg/L (IQR, 0.66–2.58). At the time of sampling, patients were receiving a median total daily VRZ dose of 100.00 mg (IQR, 100.00–133.25) for a median duration of 9.50 days (IQR, 5.00–17.75). Crucial laboratory values included a median WBC of 8.04 × 109/L (IQR, 4.71–12.38) with a median neutrophil percentage of 52.58% (IQR, 34.80–66.62), and a median HGB level of 105.50 g/L (IQR, 94.92–116.25). The median ALT and SCR were 39.50 U/L (IQR, 20.31–81.80) and 18.92 μmol/L (IQR, 15.40–22.54), respectively. The most common co-administered drugs were glucocorticoids (47 patients, 61.8%), proton pump inhibitors (34, 44.7%) and tacrolimus (25, 32.9%). For model development, the dataset was randomly partitioned into a training set and a validation set at a 7:3 ratio, with comparable characteristics detailed between the two sets in Supplementary Table S1. We assessed hepatic function parameters (ALT, AST, ALP) and found no significant changes between baseline and maximum values during voriconazole therapy (Supplementary Figure S1). Additionally, breakthrough fungal infections were not observed.
Table 1. Clinical characteristics of the study population.
3.2 PopPK analysis
In the final PopPK model assessed the inter-individual variability (IIV) for CL/F through the variance component ω2CL/F, estimated at 0.674, indicating a standard deviation of approximately 0.821, with an RSE% of 10.7%. However, IIV was not estimated for V/F due to high shrinkage (>90%) and an RSE exceeding 30%. The final model yielded a population typical value of CL/F at 17.9 (L/h/70kg) and a V/F of 788 (L/70kg), with RSE% values of 10.8% and 15.4%, respectively, as detailed in Supplementary Table S3 and Equations 1, 2. The GOF plots and VPC for our final model were illustrated in Supplementary Figures S2, S3. Bootstrapping validation (100% minimization success) of the population typical value estimates for CL/F indicated a median of 17.65 (L/h/70kg) with a 95% CI of 9.11–22.00 (L/h/70kg). Meanwhile, the median for V/F was 787.75 (L/70kg), with a 95% CI of 94.82–1050.85. The final PopPK model was as follows:
3.3 Feature selection and model construction
The Boruta algorithm, employed for covariate screening across 200 iterations, identified 10 clinically significant variables that outperformed the shadow attributes (Figure 2). These variables included PK parameters such as CL and V, hematologic indices (INR, WBC, RBC, and HGB), as well as age, BSA, weight, and total daily dose. These ten variables were subsequently selected for the further development of the ML model.
Figure 2. Boruta feature filtering each variable importance box plot. Abbreviations: CL, clearance; INR, international normalized ratio; WBC, white blood cell count; BSA, body surface area; TDOSE, Total daily dose; V, volume of distribution; WT, weight; RBC, red blood cell count; HGB, hemoglobin; L, lymphocyte percentage; TAC, tacrolimus; ALB, albumin; ALP, alkaline phosphatase; TBIL, total bilirubin; PLT, platelet count; SCR, serum creatinine; CSA, Cyclosporine A; AST, aspartate aminotransferase; SIRO, sirolimus; PPI, proton pump inhibitor; DDR, D-dimer; GCS, glucocorticoids; N, neutrophil percentage; CRP, C-reactive protein; DBIL, direct bilirubin; ALT, alanine aminotransferase; EGFR, estimated glomerular filtration rate.
The evaluation of model performance was shown in Figure 3 and Supplementary Table S4. Among the six algorithms assessed, XGBoost emerged as the optimal model, with an R2 value of 0.81, MSE of 0.28, RMSE of 0.53, and MAE of 0.40. Random Forest (RF) demonstrated comparable accuracy with an MSE of 0.29 and an R2 of 0.80, indicating a strong predictive capability. The robustness of the final model was confirmed via a 10-iteration bootstrap cross-validation, which yielded a stable mean R2 of 0.805 (range: 0.776–0.867), as detailed in Supplementary Table S5. The characteristics of the external validation population (n = 10) were summarized in Supplementary Table S2, including a median age of 14 months (IQR, 5.50–21.50) and a median weight of 7.75 kg (IQR, 6.00–10.23). Furthermore, on the external validation dataset, the XGBoost model achieved an R2 of 0.75, RMSE of 0.14, MSE of 0.02, and MAE of 0.12. The VRZ plasma concentration observations were contrasted with predictions based on the XGBoost algorithm, as shown in Figure 4. The majority of data points fall within the indicated ±20% relative error thresholds. Considering the model’s prediction accuracy and goodness of fit, the XGBoost was selected for further interpretability analysis.
Figure 3. Comparison of Model Performance Metrics. The ranges for MSE, RMSE, MAE, and R2 are all between 0 and 1, where a higher R2 indicates better model performance, and lower values of MSE, RMSE, and MAE indicate better accuracy. Abbreviations: GBDT, gradient boosting decision tree; RF, random forest. MSE, mean squared error; RMSE, root mean squared error; MAE, mean absolute error.
Figure 4. Voriconazole (VRZ) plasma concentration observations versus predicted concentration based on the XGBoost algorithm. In the scatter plot, the two dashed lines represented the ±20% relative error thresholds.
3.4 Model interpretation
The XGBoost model was employed to predict plasma concentrations of oral VRZ in children under 2 years old, alongside an analysis of model interpretability. Global feature contributions, as determined by SHAP analysis, were presented in Figure 5. Clearance (CL) was identified as the most important feature, followed by weight (WT), hemoglobin (HGB), international normalized ratio (INR), and age, based on mean absolute SHAP values (Figure 5A). The beeswarm plot (Figure 5B) visualized the magnitude and direction of these contributions for each individual, confirming the substantial and heterogeneous impact of CL and WT, as evidenced by their wide distribution of SHAP values. Figure 6 displayed the SHAP dependence plots, illustrating the relationship between feature values and their influence on predictions made by the XGBoost model. The findings demonstrated a negative correlation, suggesting that higher CL values were associated with lower VRZ concentrations.
Figure 5. The SHAP analysis of the XGBoost. (A) Bar plot of feature importance in XGBoost model predictions based on SHAP Values. (B) The beeswarm plot showed the SHAP values for ten features in the XGBoost model. The darker the color, the more important the variable was. (C) The SHAP heat force plot. (D) The SHAP waterfall plot.
Figure 6. SHAP dependence plot of the XGBoost model. The x-axis showed the SHAP values of each variable, while the y-axis displayed the raw values. SHAP dependence plots for (A) CL, clearance, (B) weight, (C) HGB, hemoglobin, (D) INR, (E) age, (F) RBC, red blood cell count (RBC), (G) TDOSE, total daily dose, (H) WBC, white blood cell count, (I) BSA, body surface area.
Figure 5C displayed a SHAP waterfall plot, which deconstructed the XGBoost model’s prediction for a pediatric patient (age 19 months, weight 10 kg) receiving oral VRZ. This visualization illustrated the stepwise derivation of the patient’s individualized concentration prediction (f(x) = 1.13 mg/L) from the population mean baseline (E [f(x)] = 1.44 mg/L). The analysis identified VRZ clearance (CL = 3.5 L/h/70kg) as the predominant driver of the prediction, contributing a SHAP value of −0.313. Conversely, the patient weight = 10 kg exerted a positive influence (SHAP = +0.16). Other covariates with negative contributions included RBC = 6.19 × 109/L (SHAP = −0.128), INR = 0.95(SHAP = −0.0635), and HGB = 102 g/L (SHAP = −0.054). Patient age (19 months) had a modest positive impact (SHAP = +0.0954). As a complementary visualization, the SHAP force plot (Figure 5D) provides a consolidated view, illustrating how these positive (red) and negative (blue) forces balance to produce the final output.
3.5 Clinical application
Figure 7 illustrated the process of predicting VRZ plasma concentration using our model. For instance, the model utilized patient-related parameters, including age (23 months), weight (11 kg), BSA (0.5 m2), WBC (6.9 × 109/L), RBC (3.95 × 1012/L), HGB (89 g/L), INR (1.07), CL/F (4.46 L/h/70kg), V/F (123.83 L/70kg), and total daily dose (180 mg). Given the initial dose of 180 mg daily, administered as 90 mg twice daily, this corresponds to a minimum dosage of 8 mg/kg (8.18 mg/kg). By inputting these parameters and running the model, the VRZ concentration was predicted to be 1.06 mg/L. Three days later, TDM results indicated an actual VRZ plasma concentration of 1.21 mg/L. This model is expected to assist clinicians in making initial VRZ dosage recommendations to rapidly achieve the therapeutic concentration window before TDM results are available.
Figure 7. Clinical application of the model.
4 Discussion
VRZ is a critical therapeutic agent for invasive fungal infections; however, evidence-based dosing guidelines and drug exposure data in children younger than 24 months are scarce. To address this knowledge gap, we developed a model integrating PopPK with ML to estimate steady-state plasma concentrations of oral VRZ in this vulnerable pediatric population. This study enrolled 76 pediatric patients, including 58 males (76.32%) and 18 females (23.68%), with a median age of 11 months and a weight of 8.1 kg. Using a Bayesian PopPK approach, we estimated individual clearance and volume of distribution, incorporating these PK parameters as features in our ML model. Prior to model development, the Boruta algorithm was used for feature selection, resulting in 10 variables being included in the final ML model. The XGBoost algorithm demonstrated optimal predictive ability for VRZ plasma concentrations in children under 2 years, with an R2 of 0.81 and an MSE of 0.28. SHAP analysis identified clearance, weight, and HGB as the top three significant predictors of VRZ concentrations.
Previous studies have investigated VRZ dosing and plasma concentrations in diverse patient populations. For instance, Gastine et al. examined VRZ plasma concentrations in 17 children under 24 months across 50 distinct treatment episodes, reporting a median trough concentration of 0.63 mg/L, with only 34.2% of samples reaching the recommended therapeutic range (Gastine et al., 2023). In contrast, our study involved a larger population of 76 pediatric patients, potentially enhancing the reliability and general applicability of our findings. Huang et al. assessed the efficacy of ML models versus traditional PopPK models for predicting VRZ trough concentrations in critically ill patients (Huang et al., 2025). Utilizing a dataset of 244 concentrations from 62 patients, they developed 6 ML models, with the XGBoost model demonstrating superior predictive performance (R2 = 0.73). Similary, Liu et al. developed a real-time ML ensemble model to forecast VCZ plasma concentrations in elderly patients, achieving an R2 of 0.828 by reducing features from 31 to nine while maintaining predictive accuracy (Liu et al., 2025). In our analysis, the XGBoost model achieved an R2 of 0.81 and an MSE of 0.28, indicating comparable predictive capability. Furthermore, our SHAP analysis evaluated the influence of several variables, including clearance, body weight, and hemoglobin levels, on VCZ concentrations, providing clinical interpretations of feature importance.
In our PopPK model, body weight was incorporated as a covariate influencing VRZ clearance. Numerous factors have been reported to affect VRZ metabolism, contributing to variability in blood concentrations. Patient demographics significantly influence voriconazole pharmacokinetics. Age-related differences are notable, with pediatric patients requiring higher weight-based dosing due to faster clearance compared to adults (Tilen et al., 2022). The PopPK study by Wang et al. identified clearance as a pivotal parameter, significantly correlating with age, ALP levels, and CYP2C19 genotype (Wang et al., 2014). Johnson et al. found that CL was significantly negatively correlated with INR, total bilirubin, and AST in liver transplant patients (P < 0.05), suggesting that liver function indicators may indirectly regulate drug concentration by affecting clearance (Johnson et al., 2010). Consistent with these findings, our final PopPK model estimated a population typical value of clearance (CL/F) at 17.9 L/h/70kg (RSE, 10.8%). Bootstrapping validation indicated a median CL/F of 17.65 L/h/70kg (95% CI, 9.11–22.00 L/h/70kg). These results further underscore the importance of clearance as a key factor of VRZ concentrations, particularly in the pediatric population.
SHAP analysis identified CL as the most significant factor influencing VRZ concentration, consistent with the substantial interindividual variability observed in clearance among individuals. This variability likely reflects underlying physiological differences such as liver function, which is primarily responsible for metabolizing voriconazole (Hashemiza et al., 2017), as well as interindividual differences in CYP450 enzyme activity. Studies have confirmed that liver dysfunction can complicate clearance metrics, necessitating adjustments to dosing to avoid subtherapeutic levels and toxicity (Gorski et al., 2011; Wang et al., 2018). While CL was the predominant factor, other features also contributed to VRZ concentration variability. Body weight (WT) plays a significant role in individualizing VRZ dosing regimens, particularly in pediatric populations. A clinical study observed that obese individuals demonstrated reduced VRZ metabolism compared to their non-obese counterparts, which was associated with higher VRZ trough concentrations (Takahashi et al., 2020). Conversely, Patients with lower body weight may experience insufficient drug exposure, leading to lower plasma drug concentrations (Tilen et al., 2022), particularly in those with a body surface area under 1 m2. Children metabolize VRZ at a faster rate due to higher metabolic activity and different body compositions compared to adults, necessitating more intensive dosing strategies to reach target plasma concentrations effectively (Neely et al., 2010; Shima et al., 2010). Hemoglobin (HGB) levels can influence pharmacokinetics indirectly through their association with blood volume and systemic inflammatory responses. Although the effect size of HGB was smaller than that of CL in our study, studies have suggested associations between HGB levels and altered pharmacokinetics, particularly in populations at risk for anemia (Hsu et al., 2015). Anemia or elevated WBC and RBC counts might reflect systemic inflammation or altered drug distribution, affecting VRZ concentrations (Neely et al., 2010). Age also significantly affects VRZ pharmacokinetics, with pediatric patients requiring different dosing strategies due to metabolic differences compared to adults. As patients age, VRZ clearance can change, necessitating tailored dosing to maintain efficacy without increasing toxicity (Chen et al., 2022; Bartelink et al., 2013). While CL remains the dominant factor, VRZ concentration variability, as revealed by our SHAP analysis, is meaningfully influenced by a combination of body weight, hemoglobin levels, age, and inflammatory markers; thus, these factors warrant consideration in individualized dosing, especially for vulnerable populations.
The study had several limitations. First, as a single-center, retrospective analysis, the findings may have limited generalizability. Furthermore, the study population comprised exclusively of Asian children; PK differences across ethnic groups suggest that the model’s applicability to non-Asian pediatric populations requires further investigation (Chen et al., 2018; Weiss et al., 2009). Additionally, Voriconazole plasma concentrations were measured using an assay established within our laboratory. Variations in reagents, instrumentation, and potential matrix effects between laboratories may lead to systematic differences in model performance at other centers.
The trough-only sampling design resulted in a wide 95% confidence interval for the apparent volume of distribution (V/F 94.82–1050.85L/70kg). While this approach optimized the estimation of voriconazole clearance (CL/F), it provided limited information regarding drug distribution. Nevertheless, the stability of the overall model structure and the reliability of key parameter estimates were supported by a 100% success rate in bootstrap validation. Moreover, our findings were consistent with previously reported high inter-individual variability and allometric scaling in pediatric populations (Friberg et al., 2012).
Additionally, the model did not include genotype features related to VRZ metabolism, such as CYP2C19, which may limit a comprehensive understanding of drug metabolism (Moriyama et al., 2017; Caudle et al., 2025). However, it is important to emphasize that CYP2C19 genotyping has not yet been established as a standard clinical practice in most healthcare settings, presenting a challenge for the widespread adoption of genotype-guided VRZ dosing strategies. Therefore, the primary objective of this study was to develop a clinically applicable model using readily available clinical data to predict VRZ concentrations in children under 2 years of age. This model is intended to help clinicians rapidly achieve therapeutic VRZ concentrations by informing initial dosage recommendations before therapeutic drug monitoring (TDM) results are available. As a National Children’s Medical Center, we are leading a multicenter, prospective study to systematically evaluate the efficacy and safety of VRZ in young children. Ultimately, we aim to develop a user-friendly web application to facilitate the model’s implementation in clinical practice.
5 Conclusion
This study represented the first effort to merge classic PopPK with ML to develop a predictive model for oral VRZ plasma concentrations in children under 24 months. A one-compartment model of voriconazole was established for this age group, revealing that weight significantly impacts CL/F and V/F. Individual CL and V were estimated using the empirical Bayesian method, which were then incorporated into the ML modeling process.
We used the Boruta algorithm for feature selection and constructed prediction models for VRZ plasma concentrations employing various ML algorithms. Among these, XGBoost demonstrated the best performance, achieving an R2 value of 0.81 and a mean squared error of 0.28. SHAP explainability analysis indicated that clearance, weight, hemoglobin, and other factors were most influential in the model.
Our findings provide a significant reference for individualized voriconazole dosing in clinical practice for children under 2 years old. A potential benefit of the model lies in its capacity to suggest an initial dose before TDM results are available, which may contribute to a more timely achievement of effective therapeutic concentrations.
Data availability statement
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding authors.
Ethics statement
The studies involving humans were approved by the ethics committee of Children’s Hospital of Fudan University [No. (2024) 321]. The studies were conducted in accordance with the local legislation and institutional requirements. The ethics committee/institutional review board waived the requirement of written informed consent for participation from the participants or the participants’ legal guardians/next of kin because this study was a retrospective analysis, and all data were sourced from the hospital’s electronic medical record system. As this study involved no interventions and posed minimal risk, informed consent was waived for the collection and analysis of the anonymous data.
Author contributions
LS: Methodology, Data curation, Writing – original draft, Formal Analysis. MH: Visualization, Data curation, Writing – original draft. XX: Validation, Supervision, Writing – review and editing. YZ: Visualization, Writing – review and editing. WW: Validation, Visualization, Writing – review and editing. XG: Supervision, Software, Writing – review and editing. GW: Supervision, Writing – review and editing. YW: Writing – review and editing, Supervision, Conceptualization. ZL: Writing – review and editing, Conceptualization, Funding acquisition, Supervision.
Funding
The author(s) declare that financial support was received for the research and/or publication of this article. This research was funded by the Shanghai Medicine and Health Development Foundation (No.20221128), Shanghai Municipal Human Resources and Social Security Bureau (No. EK00000861).
Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Generative AI statement
The author(s) declare that no Generative AI was used in the creation of this manuscript.
Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.
Publisher’s note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Supplementary material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphar.2025.1671652/full#supplementary-material
References
Asnicar, F., Thomas, A. M., Passerini, A., Waldron, L., and Segata, N. (2024). Machine learning for microbiologists. Nat. Rev. Microbiol. 22 (4), 191–205. doi:10.1038/s41579-023-00984-1
Bartelink, I. H., Wolfs, T., Jonker, M., de Waal, M., Egberts, T. C. G., Ververs, T. T., et al. (2013). Highly variable plasma concentrations of voriconazole in pediatric hematopoietic stem cell transplantation patients. Antimicrob. Agents Chemother. 57 (1), 235–240. doi:10.1128/AAC.01540-12
Caudle, K. E., Whirl-Carrillo, M., Relling, M. V., Hoffman, J. M., Donnelly, R. S., Haidar, C. E., et al. (2025). Advancing clinical pharmacogenomics worldwide through the clinical pharmacogenetics implementation consortium (CPIC). Clin. Pharmacol. Ther., cpt.70005. doi:10.1002/cpt.70005
Chabala, C., Turkova, A., Hesseling, A. C., Zimba, K. M., van der Zalm, M., Kapasa, M., et al. (2022). Pharmacokinetics of first-line drugs in children with tuberculosis, using world Health Organization-recommended weight Band doses and Formulations. Clin. Infect. Dis. 74 (10), 1767–1775. doi:10.1093/cid/ciab725
Chen, K., Zhang, X., Ke, X., Du, G., Yang, K., and Zhai, S. (2018). Individualized Medication of voriconazole: a practice guideline of the Division of therapeutic drug monitoring, Chinese Pharmacological Society. Ther. Drug Monit. 40 (6), 663–674. doi:10.1097/FTD.0000000000000561
Chen, J., Wu, Y., He, Y., Feng, X., Ren, Y., and Liu, S. (2022). Combined effect of CYP2C19 genetic polymorphisms and C-reactive protein on voriconazole exposure and dosing in immunocompromised children. Front. Pediatr. 10, 846411. doi:10.3389/fped.2022.846411
Cheng, L., Zhao, Y., Liang, Z., You, X., Jia, C., Liu, X., et al. (2023). Prediction of plasma trough concentration of voriconazole in adult patients using machine learning. Eur. J. Pharm. Sci. 188, 106506. doi:10.1016/j.ejps.2023.106506
Filion, K. B., Azoulay, L., Platt, R. W., Dahl, M., Dormuth, C. R., Clemens, K. K., et al. (2016). A multicenter observational study of Incretin-based drugs and heart failure. N. Engl. J. Med. 374 (12), 1145–1154. doi:10.1056/NEJMoa1506115
Friberg, L. E., Ravva, P., Karlsson, M. O., and Liu, P. (2012). Integrated population pharmacokinetic analysis of voriconazole in children, adolescents, and adults. Antimicrob. Agents Chemother. 56 (6), 3032–3042. doi:10.1128/AAC.05761-11
Gastine, S., Lehrnbecher, T., Muller, C., Farowski, F., Bader, P., Ullmann-Moskovits, J., et al. (2018). Pharmacokinetic modeling of voriconazole to develop an Alternative dosing regimen in children. Antimicrob. Agents Chemother. 62 (1), e01194-17. doi:10.1128/AAC.01194-17
Gastine, S. E., Rauwolf, K. K., Pieper, S., Hempel, G., Lehrnbecher, T., Tragiannidis, A., et al. (2023). Voriconazole plasma concentrations and dosing in paediatric patients below 24 months of age. Mycoses 66 (11), 969–976. doi:10.1111/myc.13643
Gorski, E., Esterly, J. S., Postelnick, M., Trifilio, S., Fotis, M., and Scheetz, M. H. (2011). Evaluation of hepatotoxicity with off-label oral-treatment doses of voriconazole for invasive fungal infections. Antimicrob. Agents Chemother. 55 (1), 184–189. doi:10.1128/AAC.01078-10
Greener, J. G., Kandathil, S. M., Moffat, L., and Jones, D. T. (2022). A guide to machine learning for biologists. Nat. Rev. Mol. Cell Biol. 23 (1), 40–55. doi:10.1038/s41580-021-00407-0
Hashemizadeh, Z., Badiee, P., Malekhoseini, S. A., Raeisi Shahraki, H., Geramizadeh, B., and Montaseri, H. (2017). Observational study of associations between voriconazole therapeutic drug monitoring, toxicity, and Outcome in liver transplant patients. Antimicrob. Agents Chemother. 61 (12), e01211-17. doi:10.1128/AAC.01211-17
Hsu, A. J., Dabb, A., and Arav-Boger, R. (2015). Autoinduction of voriconazole metabolism in a child with invasive pulmonary aspergillosis. Pharmacotherapy 35 (4), e20–e26. doi:10.1002/phar.1566
Huang, J., Mao, B., Bai, Y., Zhang, T., and Miao, C. (2020). An integrated Fuzzy C-means method for missing data imputation using Taxi GPS data. Sensors (Basel) 20 (7), 1992. doi:10.3390/s20071992
Huang, Q., Lin, X., Wang, Y., Chen, X., Zheng, W., Zhong, X., et al. (2022). Tacrolimus pharmacokinetics in pediatric nephrotic syndrome: a combination of population pharmacokinetic modelling and machine learning approaches to improve individual prediction. Front. Pharmacol. 13, 942129. doi:10.3389/fphar.2022.942129
Huang, Y., Zhou, Y., Liu, D., Chen, Z., Meng, D., Tan, J., et al. (2025). Comparison of population pharmacokinetic modeling and machine learning approaches for predicting voriconazole trough concentrations in critically ill patients. Int. J. Antimicrob. Agents 65 (2), 107424. doi:10.1016/j.ijantimicag.2024.107424
Isaac, T., Zheng, J., and Jha, A. (2012). Use of UpToDate and outcomes in US hospitals. J. Hosp. Med. 7 (2), 85–90. doi:10.1002/jhm.944
Jiang, Z., Wei, Y., Huang, W., Li, B., Zhou, S., Liao, L., et al. (2022). Population pharmacokinetics of voriconazole and initial dosage optimization in patients with talaromycosis. Front. Pharmacol. 13, 982981. doi:10.3389/fphar.2022.982981
Johnson, H. J., Han, K., Capitano, B., Blisard, D., Husain, S., Linden, P. K., et al. (2010). Voriconazole pharmacokinetics in liver transplant recipients. Antimicrob. Agents Chemother. 54 (2), 852–859. doi:10.1128/AAC.00429-09
Kadam, R. S., and Van Den Anker, J. N. (2016). Pediatric clinical Pharmacology of voriconazole: role of pharmacokinetic/Pharmacodynamic modeling in Pharmacotherapy. Clin. Pharmacokinet. 55 (9), 1031–1043. doi:10.1007/s40262-016-0379-2
Li, G., Li, Q., Zhang, C., Yu, Q., Zhou, X., et al. (2023). The impact of gene polymorphism and hepatic insufficiency on voriconazole dose adjustment in invasive fungal infection individuals. Front. Genet. 14, 1242711. doi:10.3389/fgene.2023.1242711
Li, G., Sun, Y., and Zhu, L. (2024). Application of machine learning combined with population pharmacokinetics to improve individual prediction of vancomycin clearance in simulated adult patients. Front. Pharmacol. 15, 1352113. doi:10.3389/fphar.2024.1352113
Li, X. H., Yang, X. L., Dong, B. B., and Liu, Q. (2025). Predicting 28-day all-cause mortality in patients admitted to intensive care units with pre-existing chronic heart failure using the stress hyperglycemia ratio: a machine learning-driven retrospective cohort analysis. Cardiovasc Diabetol. 24 (1), 10. doi:10.1186/s12933-025-02577-z
Liu, R., Ma, P., Chen, D., Yu, M., Xie, L., Zhao, L., et al. (2025). A real-time plasma concentration prediction model for voriconazole in elderly patients via machine learning combined with population pharmacokinetics. Drug Des. Devel Ther. 19, 4021–4037. doi:10.2147/DDDT.S495050
Ma, P., Shang, S., Huang, Y., Liu, R., Yu, H., Zhou, F., et al. (2024). Joint use of population pharmacokinetics and machine learning for prediction of valproic acid plasma concentration in elderly epileptic patients. Eur. J. Pharm. Sci. 201, 106876. doi:10.1016/j.ejps.2024.106876
Moriyama, B., Obeng, A. O., Barbarino, J., Penzak, S. R., Henning, S. A., Scott, S. A., et al. (2017). Clinical Pharmacogenetics implementation Consortium (CPIC) guidelines for CYP2C19 and voriconazole therapy. Clin. Pharmacol. Ther. 102 (1), 45–51. doi:10.1002/cpt.583
Neely, M., Rushing, T., Kovacs, A., Jelliffe, R., and Hoffman, J. (2010). Voriconazole pharmacokinetics and pharmacodynamics in children. Clin. Infect. Dis. 50 (1), 27–36. doi:10.1086/648679
Pappas, P. G., Kauffman, C. A., Andes, D. R., Clancy, C. J., Marr, K. A., Ostrosky-Zeichner, L., et al. (2016). Clinical practice guideline for the management of Candidiasis: 2016 Update by the infectious diseases Society of America. Clin. Infect. Dis. 62 (4), e1–e50. doi:10.1093/cid/civ933
Qi, X., Wang, S., Fang, C., Jia, J., Lin, L., and Yuan, T. (2025). Machine learning and SHAP value interpretation for predicting comorbidity of cardiovascular disease and cancer with dietary antioxidants. Redox Biol. 79, 103470. doi:10.1016/j.redox.2024.103470
Resztak, M., Sobiak, J., and Czyrski, A. (2021). Recent Advances in therapeutic drug monitoring of voriconazole, Mycophenolic acid, and vancomycin: a literature review of pediatric studies. Pharmaceutics 13 (12), 1991. doi:10.3390/pharmaceutics13121991
Schwartz, G. J., Feld, L. G., and Langford, D. J. (1984). A simple estimate of glomerular filtration rate in full-term infants during the first year of life. J. Pediatr. 104 (6), 849–854. doi:10.1016/s0022-3476(84)80479-5
Schwartz, G. J., Munoz, A., Schneider, M. F., Mak, R. H., Kaskel, F., Warady, B. A., et al. (2009). New equations to estimate GFR in children with CKD. J. Am. Soc. Nephrol. 20 (3), 629–637. doi:10.1681/ASN.2008030287
Shen, L., Wu, J., Lan, J., Chen, C., Wang, Y., and Li, Z. (2024). Interpretable machine learning-based prediction of 28-day mortality in ICU patients with sepsis: a multicenter retrospective study. Front. Cell Infect. Microbiol. 14, 1500326. doi:10.3389/fcimb.2024.1500326
Shima, H., Miharu, M., Osumi, T., Takahashi, T., and Shimada, H. (2010). Differences in voriconazole trough plasma concentrations per oral dosages between children younger and older than 3 years of age. Pediatr. Blood Cancer 54 (7), 1050–1052. doi:10.1002/pbc.22451
Soler-Palacin, P., Frick, M. A., Martin-Nalda, A., Lanaspa, M., Pou, L., Roselló, E., et al. (2012). Voriconazole drug monitoring in the management of invasive fungal infection in immunocompromised children: a prospective study. J. Antimicrob. Chemother. 67 (3), 700–706. doi:10.1093/jac/dkr517
Strutz, S., Liang, H., Carey, K., Bashiri, F., Jani, P., Gilbert, E., et al. (2025). Machine learning for predicting critical Events among hospitalized children. JAMA Netw. Open 8 (5), e2513149. doi:10.1001/jamanetworkopen.2025.13149
Takahashi, T., Smith, A. R., Jacobson, P. A., Fisher, J., Rubin, N. T., and Kirstein, M. N. (2020). Impact of Obesity on voriconazole pharmacokinetics among pediatric hematopoietic cell transplant recipients. Antimicrob. Agents Chemother. 64 (12), e00653-20. doi:10.1128/AAC.00653-20
Tilen, R., Paioni, P., Goetschi, A. N., Goers, R., Seibert, I., Müller, D., et al. (2022). Pharmacogenetic analysis of voriconazole treatment in children. Pharmaceutics 14 (6), 1289. doi:10.3390/pharmaceutics14061289
Touw, D. J., and van den Anker, J. N. (2022). Therapeutic drug monitoring of Antimicrobial drugs in Neonates: an Opinion article. Ther. Drug Monit. 44 (1), 65–74. doi:10.1097/FTD.0000000000000919
van den Born, D. A., Martson, A. G., Veringa, A., Punt, N. C., van der Werf, T. S., Alffenaar, J. W. C., et al. (2023). Voriconazole exposure is influenced by inflammation: a population pharmacokinetic model. Int. J. Antimicrob. Agents 61 (4), 106750. doi:10.1016/j.ijantimicag.2023.106750
Wang, T., Chen, S., Sun, J., Cai, J., Cheng, X., Dong, H., et al. (2014). Identification of factors influencing the pharmacokinetics of voriconazole and the optimization of dosage regimens based on Monte Carlo simulation in patients with invasive fungal infections. J. Antimicrob. Chemother. 69 (2), 463–470. doi:10.1093/jac/dkt369
Wang, T., Yan, M., Tang, D., Xue, L., Zhang, T., Dong, Y., et al. (2018). A retrospective, multicenter study of voriconazole trough concentrations and safety in patients with Child-Pugh class C cirrhosis. J. Clin. Pharm. Ther. 43 (6), 849–854. doi:10.1111/jcpt.12724
Wang, H., Shen, Y., Luo, X., Jin, L., Zhu, H., and Wang, J. (2025). Population pharmacokinetics and dose optimization of voriconazole in patients with COVID-19-associated pulmonary aspergillosis. Front. Pharmacol. 16, 1554370. doi:10.3389/fphar.2025.1554370
Weiss, J., Ten Hoevel, M. M., Burhenne, J., Walter-Sack, I., Hoffmann, M. M., Rengelshausen, J., et al. (2009). CYP2C19 genotype is a major factor contributing to the highly variable pharmacokinetics of voriconazole. J. Clin. Pharmacol. 49 (2), 196–204. doi:10.1177/0091270008327537
Keywords: voriconazole, children, under 2 years, machine learning, SHAP analysis
Citation: Shen L, Hu M, Xu X, Zhou Y, Wu W, Ge X, Wang G, Wang Y and Li Z (2025) Precision dosing of voriconazole in immunocompromised children under 2 years: integrated machine learning and population pharmacokinetic modeling. Front. Pharmacol. 16:1671652. doi: 10.3389/fphar.2025.1671652
Received: 23 July 2025; Accepted: 01 September 2025;
Published: 15 September 2025.
Edited by:
Karel Allegaert, Faculty of Medicine, KU Leuven, BelgiumReviewed by:
Catherine M. T. Sherwin, University of Western Australia, AustraliaHui Xie, First Affiliated Hospital of Guangzhou Medical University, China
Copyright © 2025 Shen, Hu, Xu, Zhou, Wu, Ge, Wang, Wang and Li. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Zhiping Li, enBsaUBmdWRhbi5lZHUuY24=; Yi Wang, eWl3YW5nQHNobXUuZWR1LmNu
†These authors share first authorship