Application of a Physiologically Based Pharmacokinetic Model to Characterize Time-dependent Metabolism of Voriconazole in Children and Support Dose Optimization

Background: Voriconazole is a potent antifungal drug with complex pharmacokinetics caused by time-dependent inhibition and polymorphisms of metabolizing enzymes. It also exhibits different pharmacokinetic characteristics between adults and children. An understanding of these alterations in pharmacokinetics is essential for pediatric dose optimization. Objective: To determine voriconazole plasma exposure in the pediatric population and further investigate optimal dosage regimens. Methods: An adult and pediatric physiologically based pharmacokinetic (PBPK) model of voriconazole, integrating auto-inhibition of cytochrome P450 3A4 (CYP3A4) and CYP2C19 gene polymorphisms, was developed. The model was evaluated with visual predictive checks and quantitative measures of the predicted/observed ratio of the area under the plasma concentration-time curve (AUC) and maximum concentration (Cmax). The validated pediatric PBPK model was used in simulations to optimize pediatric dosage regimens. The probability of reaching a ratio of free drug (unbound drug concentration) AUC during a 24-h period to minimum inhibitory concentration greater than or equal to 25 (fAUC24h/MIC ≥ 25) was assessed as the pharmacokinetic/pharmacodynamic index. Results: The developed PBPK model well represented voriconazole's pharmacokinetic characteristics in adults; 78% of predicted/observed AUC ratios and 85% of Cmax ratios were within the 1.25-fold range. The model maintained satisfactory prediction performance for intravenous administration in pediatric populations after incorporating developmental changes in anatomy/physiology and metabolic enzymes, with all predicted AUC values within 2-fold and 73% of the predicted Cmax within 1.25-fold of the observed values. The simulation results of the PBPK model suggested that different dosage regimens should be administered to children according to their age, CYP2C19 genotype, and infectious fungal genera. Conclusion: The PBPK model integrating CYP3A4 auto-inhibition and CYP2C19 gene polymorphisms successfully predicted voriconazole pharmacokinetics during intravenous administration in children and could further be used to optimize dose strategies. The infectious fungal genera should be considered in clinical settings, and further research with large sample sizes is required to confirm the current findings.


INTRODUCTION
Voriconazole is an essential triazole antifungal agent with in vivo activity against a broad spectrum of yeasts and filamentous fungi, commonly used for the prophylaxis and treatment of various invasive fungal infections (IFI) (Clancy and Nguyen, 1998;Saravolatz et al., 2013;Perfect et al., 2003). However, a high interindividual plasma variability has been observed partially due to its nonlinear and time-dependent pharmacokinetics (Purkins et al., 2003;Pfizer, 2010). It also exhibits differences in clearance and bioavailability between adults and children (Schulz et al., 2019). All these factors complicate the successful therapeutic use of voriconazole in pediatric populations.
Voriconazole is a substrate for cytochrome P450 (CYP) enzymes. CYP2C19 is the primary enzyme that contributes to the main circulating metabolite of voriconazole, voriconazole N-oxide, while CYP3A4 and CYP2C9 are also responsible for its metabolism (Hyland et al., 2003). Only 2% of voriconazole is excreted unmetabolized in the urine (Levêque et al., 2006). The time-dependent inhibition (TDI) of CYP3A4 observed in in vitro studies (Li et al., 2020) may play a role in the elevated exposure to voriconazole at multiple doses, which could not be predicted from single-dose data (Purkins et al., 2003). Genetic polymorphisms of CYP2C19 are also a major determinant of the wide pharmacokinetic (PK) variability in voriconazole. Drug exposures from multiple doses of poor metabolizers (PMs) are almost 3 times higher than those from normal metabolizers (NMs) (Lee et al., 2012). Regarding age-related changes, the total body clearance in children is approximately 2-3-fold higher than that in adults (Levêque et al., 2006). In addition, oral bioavailability in children (45-66%) is only half of that in adults (96%) (Schulz et al., 2019), and this may suggest that primarily gut wall metabolism is also increased in children (Zane and Thakker, 2014). These PK discrepancies may be explained by developmental differences in organs, tissues, and enzymes (Kearns et al., 2003), resulting in an increased ratio of hepatic mass to total body mass and a higher clearance of CYP enzymes in children (Walsh et al., 2010).
The physiologically based pharmacokinetic (PBPK) model combines the knowledge of system-specific factors and drugspecific factors with mathematical modeling methods to quantitatively predict the PK characteristics of drug absorption, distribution, metabolism, and excretion (Maharaj and Edginton, 2014). Previously, an adult and pediatric PBPK model was established with hepatic in vitro data (Zane and Thakker, 2014). However, this model could not predict the nonlinear PKs of voriconazole and alterations in its metabolism over time. In another study, a whole-body PBPK model of voriconazole integrating the TDI of CYP3A4 and genetic polymorphisms of CYP2C19 was constructed (Li et al., 2020). It successfully captured the main PK characteristics of the drug in adults but overpredicted exposure to PMs after multiple doses.
Therefore, the objectives of the present study were to 1) establish an adult PBPK model of voriconazole, focusing on improving predictions of multiple-dose administration, especially in PMs; 2) extrapolate this model to children using age-related scaling methods; and 3) conduct simulations to facilitate the dose-optimization process. Due to the lack of data and high interpatient variability in the PK parameters observed following oral (p.o.) administration of voriconazole in children, the adult model was only extrapolated to intravenous (i.v.) administration in the pediatric model.

Workflow and Software
In this study, a PBPK model of voriconazole was developed and evaluated in the adult population, subsequently extrapolated to the pediatric population, and verified by comparing the simulated plasma exposure with the observed data. The final PBPK model was then used to simulate PK studies for pediatric dose optimization. The workflow of the model development is presented in Figure 1.
The modeling work was conducted in PK-Sim ® (part of Open Systems Pharmacology (OSP) Suite version 8.0, www. open-systems-pharmacology.org). The published data were digitized by GetData Graph Digitizer version 2.26 (getdatagraph-digitizer.com).  Li et al. (2020) fu, fraction of bound drug; logP, log of the partition coefficient between octanol and water; pKa, acid dissociation constant; K m , Michaelis-Menten constant; k cat , in vitro V max per recombinant enzyme; GFR, glomerular filtration rate; k inact , maximum inactivation rate constant; K I , the inhibition concentration when reaching 50% of k inact ; D T,50 , the dissolution time when 50% of the substance dissolved; shape factor, the dissolution shape parameter for Weibull function.
Drug-specific parameters such as physicochemical properties and PK characteristics describing absorption, distribution, metabolism, and excretion were obtained from the literature, as shown in Table 1. Tissue-to-plasma partition coefficients were calculated using the Poulin and Theil method. Cellular permeabilities were predicted using the PK-Sim standard method. The enzymatic clearance process was quantified using Michaelis-Menten kinetics (Michaelis et al., 2011). The initial value of the Michaelis-Menten constant (K m ) for CYP2C19, CYP3A4, and CYP2C9 was 3.5, 11, and 20 μmol/L, and the maximum velocity (V max ) was set to 1.19, 2.3, and 0.0556 pmol/min −1 /pmol, respectively (Hyland et al., 2003;Damle et al., 2011). As the ratio of AUC within the dosing interval after multiple doses to AUC from zero to infinity after a single dose (AUC τ(multipole dose) /AUC inf(single dose) ) is larger than 2 (Purkins et al., 2003), the time-dependent inhibition of CYP3A4 was integrated into the model using Eq. 1.
dE cat /dt describes the turnover of the enzyme, where k deg is turnover rate constant, E 0 is the initial enzyme concentration, I is the inhibitor concentration, k inact is the maximum inactivation rate constant, and K I is the inhibition concentration when reaching 50% of k inact . The initial values of k inact and K I were set to 0.04 min −1 and 9.33 μmol/L, respectively, according to the results of an in vitro inactivity assay (Li et al., 2020). Other parameters related to TDI were set as software default values. The glomerular filtration rate (GFR) fraction was fixed to 1, as there was no evidence for reabsorption and tubular secretion. Initial formulation-related parameters were also obtained from the literature, with the dissolution time when 50% of the substance dissolved (D T,50 ) set to 30 min and dissolution shape parameter for Weibull function (shape factor for tablet) set to 1.29 (Li et al., 2020). The PBPK model was first established based on the initial values indicated above, and the V max of CYP3A4 was optimized based on a single i.v. dose in CYP2C19 PMs, assuming that CYP3A4 contributes to almost 100% of the metabolism in this population. The propriety of CYP2C19-related parameters was then evaluated based on the observed PK parameters and plasma concentration profiles from a single i.v. dose in CYP2C19 RMs/ NMs/IMs. In the next step, the specific intestinal permeability and formulation-related parameters, including D T,50 and shape factor for tablet, were inspected using data from studies following a single p.o. administration in RMs/NMs/IMs/PMs. If the fits were deemed inadequate, these three parameters were optimized based on p.o. data. Data from Scholz et al. (2009) was used for the above parameter optimization. The TDI-related parameter was optimized in the final step due to the lack of multiple i.v. clinical studies with genotype information. If the predicted PK parameters from multiple p.o. studies in PMs were within 2fold of the observed values, the modeling process was complete. Otherwise, k inact for CYP3A4 was optimized and the iterations of the above optimization steps were repeated.

Adult Model Verification
The PK simulation of 100 virtual people for each clinical study was carried out, corresponding to the subject demographics (age range, proportion of male/female, and dosing regimens). The predictive performance was evaluated by visually comparing predicted concentration-time data with the observed data from the literature for initial verification. Ninety percent population prediction intervals covering the observed plasma concentrationtime datasets were considered as a visual criterion for good performance. Next, the quantitative assessment was conducted by calculating the mean fold error (MFE) of PK parameters such as the area under the plasma concentration-time curve (AUC) and maximum concentrations (C max ), expressed as the ratio of predicted to observed mean values. The model was acceptable if it met the 0.5-to 2.0-fold limit, and a more stringent criterion was the 0.8-to 1.25-fold range.

Pediatric Model Development
PK profiles following i.v. administration were searched in PubMed, Embase and the Cochrane Library using the terms "voriconazole", "pharmacokinetic", and children-related items "infant", "child", "children", "pediatric", and "adolescent". Studies with sufficient information for dosage, concentration, and CYP2C19 genotype were included.
Drug-specific parameters defined in the adult PK data were kept constant for the pediatric PBPK model.
Developmental changes in anatomical and physiological parameters such as weight, height, organ volumes, blood flows, organ composition, and plasma protein concentration in PK-Sim ® are based on the population data from previous studies (NHANES III, 1997;ICRP, 2002  CYP2C9 ontogeny information is described in the online PK-Sim Ontogeny Database Open Systems Pharmacology (2018). In summary, the activity of CYP2C19, CYP3A4 and CYP2C9 increases after birth and reaches the adult level over approximately 1, 4 and 1 year, respectively. The model using these default ontogeny profiles overpredicted the exposure in children; therefore, the ontogeny factors of CYP enzymes were estimated based on a recent meta-analysis using in vivo data (Upreti and Wahlstrom., 2016; Supplementary Figure S1).
These data suggested that in childhood, CYP2C19, CYP3A4, and CYP2C9 exhibit maximal activities beyond the levels in adults.

Pediatric Model Verification
The PBPK model performance in children was evaluated using the quantitative verification described in Adult Model Verification.
As some pediatric clinical trials (Walsh et al., 2010) did not clearly show results according to the CYP2C19 genotype as in adults, these data were verified by adding up different genotype results in

Pediatric Dose Optimization
A 100-patients simulation was generated for each subpopulation classified by age (2-6 and 6-12 years) and CYP2C19 genotype (NMs, IMs and PMs). The individual AUC 24h (AUC during a 24-h period) was estimated. The probability of reaching a ratio of free drug (unbound drug concentration) AUC 24h to minimum inhibitory concentration greater than or equal to 25 (fAUC 24h /MIC ≥ 25) was considered to be the PK/pharmacodynamic (PD) index (Wang et al., 2015). The fraction of unbound drug was set to 42% (US FDA, 2005), with the assumption that this value is similar between children and adults (Yanni et al., 2008). Voriconazole MIC distributions for Aspergillus (4 species) and Candida (14 species) infections were obtained from the European Committee on Antimicrobial Susceptibility Testing database (2020); (Supplementary Figure  S2). The probability of target attainment (PTA) at each MIC and the cumulative fraction of response (CFR) for the overall MIC distribution for each species were calculated. CFR values of all species in one genus larger than or equal to 80% was considered an appropriate dosage regimen (Mangal et al., 2018).

Adult Model Verification
Input parameters of voriconazole PBPK model were summarized in Table 1. A proper fit was achieved in the adult model, as shown in Figure 2, and most of the observed concentrations fell within the 5% and 95% concentration-time prediction intervals. All PK parameters were predicted to be within the 2-fold difference, with 78% of predicted/observed AUC ratios and 85% of C max ratios within the 1.25-fold difference (Table 2; Figure 3). The prediction of PMs following multiple doses was satisfactory.

Pediatric Model Verification
The PBPK model well captured the pharmacokinetic features in children after integrating in vivo ontogeny profiles. The

Pediatric Dose Optimization
Twelve pediatric groups were constructed for dose design, and the results are shown in Table 4. IMs with *17 have a similar CYP2C19 abundance to NMs and the recommended dosage for NMs can be a reference for this subpopulation. Therefore, the following IM dosage recommendation is for IMs with *1. Figure 4 shows the contrast of the minimum PTA among different species in two fungal genera at each specific MIC between the following scenarios: one is administered the recommended maintenance dose from the current medication label (8 mg/kg, twice daily (BID)), and the other one is administered the recommended dose determined from the PBPK model. The results suggested that, for the BID dosing regimens, intravenous doses should be adjusted to 12 mg/kg for NMs, 8 mg/kg for IMs, and 5 mg/kg for PMs for 2-6-year-old children infected with Aspergillus spp. As children grow older, the recommended dose decreases, specifically 9 mg/kg for NMs, 6 mg/kg for IMs, and 4 mg/kg for PMs infected with Aspergillus spp. When the infectious fungal genus is Candida spp., approximately half of the above dosages are recommended to attain the appropriate drug exposure using fAUC 24h /MIC as the indicator, due to notable differences in these two fungi in terms of susceptibility to voriconazole.

DISCUSSION
The adult and pediatric PBPK model of voriconazole, which incorporated the TDI of CYP3A4, gene polymorphisms of CYP2C19, and developmental changes in physiology and metabolic enzymes, were able to describe the PKs in both populations. Subsequent simulations revealed that age, CYP2C19 genotype, and infectious fungal genera influence target PK/PD index attainment and should be considered in dose optimization.  Voriconazole is often for long-term use, from weeks to months, in the prophylaxis and treatment of IFI (Benitez and Carver, 2019). Therefore, it is essential to elucidate the timedependent PK characteristics of this medication following multiple doses. To accomplish this goal, the strategy used for model building in this study is slightly different from that used in other studies. TDI-related parameters have always been optimized in the final step in multiple studies. However, the resulting change in these values can affect the previous predictions of single-dose administration. Generally, retrograde clearance can be implemented in the model first, and then sensitivity analysis is conducted to determine the optimal value of uncertain parameters (Wu et al., 2014). However, due to the possible inaccuracies of data digitized from the literature, which may influence the sensitivity analysis, this method was not used. Instead, a "bottom-up" and "top-down" combination strategy incorporating manual iterations was used to facilitate the model building process. The adult PBPK model results showed satisfactory predictive performance for multiple-dose administration in all CYP2C19 genotype populations.
There was an overprediction of exposure in children when default enzyme ontogeny profiles from in vitro experiments were utilized to extrapolate the adult model. Therefore, ontogeny factors based on a meta-analysis of in vivo CYP activity with age-related changes (Upreti and Wahlstrom, 2016) were implemented. In this meta-analysis, the maximal activity of CYP2C19, CYP3A4, and CYP2C9 in children exceeded the corresponding adult values. The higher capacity or expression of CYP enzymes in children was also corroborated from a previous in vitro experiment using human liver microsomes (Yanni et al., 2010). That experiment showed that the apparent V max for voriconazole conversion to the N-oxide metabolite was approximately 3-fold higher in children than in adults. After incorporating the new ontogeny profiles, the predictive performance of the PBPK model was improved. Moreover, the clearance was approximately 2-3-fold higher in children than in adults in the simulation, which is consistent with a previous study (Walsh et al., 2004).
The final extrapolated model presented the PK features of voriconazole in the pediatric population following i.v. administration, except that the exposure of several 6-12-year-old children in the study by Walsh et al. was slightly overpredicted. There may be several reasons. First, the predicted PK parameters were calculated from a cohort of children with a certain percentage of different CYP2C19 metabolic types described in the baseline demographic data. However, a few children discontinued treatment due to adverse events. The metabolic types of these children were not disclosed, which may result in differences between the predicted and observed values. Furthermore, although children who were receiving drugs known to interact strongly with voriconazole, such as terfenadine and pimozide, were excluded, concomitant medications that could potentially interact with voriconazole were permitted, which might have resulted in altered PKs that were not considered in the PBPK model. The results of dose optimization seem to be reasonable. Most of the recommended dosages for NM and IM children aged 2-12 years infected with Aspergillus spp. are similar to the dosages adopted by the European Medicines Agency (EMA) for children aged 2-11 years, specifically a loading i.v. dose of 9 mg/kg BID on day 1, followed by a maintenance dose of 8 mg/kg BID (European Medicines Association, 2012). They are also close to the loading dose of 7 mg/kg BID recommended by the Infectious Disease Society of America (IDSA) (Pappas et al., 2009), except for 12 mg/kg, which is the recommended dosage for 2-6-year-old NMs. Although the tolerance of higher dosage up to 10 mg/kg in children has been confirmed (Sano et al., 2016), further validation of this regimen's exposure and safety is needed before implementing it in clinical practice. In other subpopulations, dose adjustment is proposed based on several factors. First, younger children may have a relatively higher enzyme expression or activity (Upreti and Wahlstrom., 2016), contributing to higher metabolism and lower concentration. Thus, higher recommended doses in younger children, especially <5 years old, seem warranted (Soler-Palacín et al., 2011). Second, as CYP2C19 is the major enzyme in the metabolism of voriconazole (Barbarino et al., 2018), lower doses should be administered to people with alleles associated with enzymatic loss-of-function (Moriyama et al., 2017). Third, invasive fungal infection classification should be considered a key factor for antifungal therapy (Wang et al., 2016). The susceptibility of Candida spp. to voriconazole is higher than that of Aspergillus spp.; therefore, a reduced dose of voriconazole should be sufficient for infections involving the former (Perez-Pitarch et al., 2019).
The recommended dosage in this study should be deemed a reference for the initial dosing regimen because the target CFR value represented the overall distribution of MIC in the population and the clinical data utilized for model validation were mean drug exposure. It cannot completely solve the problem of high PK variability and replace the importance of therapeutic drug monitoring (Gastine et al., 2017). Subsequent dose adjustments should be conducted based on the individual MIC, drug concentration and clinical response (Muto et al., 2015).
One limitation of this study is that the model was established based on some fundamental assumptions. For example, CYP2C19, CYP3A4, and CYP2C9 account for the entire metabolism of voriconazole, and the metabolic pathway is the same in both adults and children. A study of human liver microsomes provided evidence that flavin-containing monooxygenase 3 (FMO3) contributes to higher clearance in children than in adults (Yanni et al., 2010). However, this was not integrated into our model due to the small effect on the main metabolite of voriconazole observed in previous recombinant FMO3 enzyme experiments (Yanni et al., 2008). Moreover, the limited available data on RMs made it impossible to recommend dosage for this subpopulation in children. The established model and the recommended dosage for other subpopulations may also require further verification in carefully designed clinical trials.
In conclusion, the developed PBPK model of voriconazole provides a more accurate description of PK characteristics in adults following single and multiple i.v. and p.o. administrations, especially in PMs. It also predicts exposure in children following i.v. administration with good accuracy. Age, CYP2C19 genotype, and infectious fungal genera were found to significantly influence the attainment of the PK/PD target in the simulation and thus should be considered for clinical dose selection.

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 author.