Model-based clinical dose optimization for phenobarbital in neonates: An illustration of the importance of data sharing and external validation

most of all available data. In this study, we utilize previously collected therapeutic drug monitoring (TDM) data of term and preterm newborns to develop a population pharmacokinetic model for phenobarbital. We externally validate the model using prospective phenobarbital data from an ongoing pharmacokinetic study in preterm neonates. Methods: TDM data from 53 neonates (gestational age (GA): 37 (24 – 42) weeks, bodyweight: 2.7 (0.45 – 4.5) kg; postnatal age (PNA): 4.5 (0 − 22) days) contained information on dosage histories, concentration and covariate data (including birth weight, actual weight, post-natal age (PNA), postmenstrual age, GA, sex, liver and kidney function, APGAR-score). Model development was carried out using NONMEM ® 7.3. After assessment of model ﬁ t, the model was validated using data of 17 neonates included in the DINO (Drug dosage Improvement in NeOnates)-study. Results: Modelling of 229 plasma concentrations, ranging from 3.2 to 75.2 mg/L, resulted in a one compartment model for phenobarbital. Clearance (CL) and volume (V d ) for a child with a birthweight of 2.6 kg at PNA day 4.5 was 0.0091 L/h (9%) and 2.38 L (5%), respectively. Birthweight and PNA were the best predictors for CL maturation, increasing CL by 36.7% per kg birthweight and 5.3% per postnatal day of living, respectively. The best predictor for the increase in V d was actual bodyweight (0.31 L/kg). External validation showed that the model can adequately predict the pharmacokinetics in a prospective study. Conclusion: Data-sharing can help to successfully develop and validate population pharmacokinetic models in neonates. From the results it seems that both PNA and bodyweight are required to guide dosing of phenobarbital in term and preterm neonates.


Introduction
Rational dosing guidelines for drugs in neonates are urgently needed. However, datasets from prospective clinical trials in children are scarce and both the number of included children and the number of samples per child are usually very small . To overcome this problem data-sharing is of utmost importance and can help to make the most out of all available data Ince et al., 2009;. Existing data can be utilized to determine the optimal design of prospective trials but it may also aid dose finding in ongoing trials in case the collected data is not (yet) sufficient to draw valid conclusions (Ince et al., 2009;Krekels et al., 2015).
Besides that, the application of advanced data analysis techniques, namely the population pharmacokinetic modelling approach, allows handling sparse and infrequently collected samples . Moreover, it offers the possibility to quantify interindividual variability and to identify covariates that determine the pharmacokinetics of drugs along the whole pediatric life-span and can thereby be used to optimize drug dosing (Brussee et al., 2016;De Cock et al., 2011).
Phenobarbital remains the traditional first-line treatment for seizures in neonates although evidence to favour one antiepileptic agent over the other is lacking (Blume et al., 2009). Using phenobarbital alone, only around 50% of seizures can be effectively controlled (Painter et al., 1999;Tulloch et al., 2012). As persistence of seizures might cause permanent functional and structural damage to the brain and existing brain damage might be worsened (Wirrell et al., 2001;Rennie and Boylan, 2003), safe and efficacious treatment is of primary importance. The optimal dosage of phenobarbital in term and preterm neonates remains a topic of discussion. The therapeutic effect is dose dependent with a suggested therapeutic range between 15 and 40 mg/L (Gilman et al., 1989). At higher concentrations sedation and feeding difficulties might occur (Ouvrier and Goldsmith, 1982). The Dutch National Children's Formulary recommends a loading dose of 20 mg/kg and a maintenance dose of 2.5-5 mg/kg/day (Dutch National Children's Formulary, 2016). However, pharmacokinetic data is sparse in term (Ouvrier and Goldsmith, 1982;Donn et al., 1985;Heimann and Gladtke, 1977;Marsot et al., 2014;Yukawa et al., 2011) and preterm (Marsot et al., 2014;Yukawa et al., 2011;Oztekin et al., 2013) newborns.
In this analysis, we utilize therapeutic drug monitoring (TDM) data collected between 1997 and 2003 in the neonatal intensive care unit of the Maastricht University Medical Centre (study 1) to build a population pharmacokinetic model for phenobarbital in term and preterm newborns. We validate the model with data originating from an ongoing PK study (DINO-study: Drug dosage Improvement in preterm NeOnates, NL47409.078.14) (study 2). This study collects pharmacokinetic data of phenobarbital and eight other frequently used off-label drugs in preterm neonates to increase the knowledge on the pharmacokinetics and pharmacodynamics using sparse sampling and limited sample volumes to minimize the burden on the individual child. Using this approach, we illustrate that data sharing and external validation can lead to model-based clinical dose optimization, particularly in neonates where ethical and practical constraints limit the possibilities to perform studies.

Model development dataset (TDM data, study 1)
TDM data (study 1) were obtained from the database of the Maastricht University Medical Centre between 1997 and 2003 with approval from the medical ethical committee (MEC 02-204.3). Neonates younger than 35 days at the start of phenobarbital treatment were included. A total of 229 samples from 53 neonates (28 male, 25 female) were available for analysis. First dose and consecutive doses (intravenous as well as oral) varied between 4 and 40.7 mg/kg and 1.3 and 20 mg/kg, respectively, and the study period ranged from 4 to 85 days (Table 1). The median first dose was 20 mg/kg, 17 children received a first dose that was higher than 25 mg/kg. The median maintenance dose was 3.9 mg/kg. Covariates were retrieved from the patient's records.

External model validation dataset (prospective data, study 2)
The DINO-study (NL47409.078.14, MEC-2014-067, study 2) prospectively studies a total of nine drugs including phenobarbital used as standard of care in preterm infants born before 32 weeks of gestation aiming at evidence-based individualized dosing regimen and is still ongoing. From September 2014 to September 2016, 61 blood samples from 17 children (7 female, 10 male) containing phenobarbital were evaluable for the analysis (Table 1).

Bioanalytical analysis
Phenobarbital concentrations of the model development dataset (study 1) were determined using a fluorescence polarization assay (Steijns et al., 2002) on the COBAS INTEGRA 700 (Roche Diagnostics; Basel, Switzerland) using COBAS INTEGRA reagent system cassettes at the pharmacy of Maastricht University Medical Centre. Fluoresceinlabelled phenobarbital binds an antibody and the emitted light is polarized due to the reduction in freedom of rotation. In case phenobarbital is present in the patients serum it reduces the extent of fluorescence polarization due to antibody binding (Steijns et al., 2002). The test range of the assay was 0.6-60 mg/L (coefficient of variation (CV) intra-assay: 0.8-2.8%, CV inter-assay: 2.9-6.5%).
Phenobarbital concentrations of the validation dataset (study 2) were determined using a particle-enhanced turbidimetric inhibition immunoassay (Petinia) on the Architect C4000 (Abbott Diagnostics; Hoofddorp, The Netherlands) using Architect C4000 reagent system cassettes at the pharmacy of the ErasmusMC, Rotterdam. The assay is based on competition between drug in the sample and drug coated onto a microparticle for antibody binding sites of the phenobarbital antibody reagent. The test range is 2.0-80 mg/L (CV intra-assay: 0.9-4.8%, CV inter-assay CV: 3.2-4.8%). A recent method comparing different immunoassays for phenobarbital found a Pearson correlation coefficient of 0.989 between a COBAS and an Architect system (Shipkova et al., 2014).

Population pharmacokinetic analysis
The analysis was performed using NONMEM ® version 7.3 (ICON Development Solutions, Ellicott City, MD, USA), supported by Perl- speaks-NONMEM (PsN ® ) version 3.4.2 and Xpose version 4.3.5 (Keizer et al., 2013). Model development started using the intravenous data. After covariate inclusion and assessment of the model fit, oral data was added to the model. The absorption rate constant (k a ) was fixed to a value of 50 h − 1 obtained from literature . The first-order conditional estimation with interaction method was used throughout model development. In case of missing covariate information for renal function or liver function, the last value observed in the subject was carried forward. For actual bodyweight, linear interpolation between available measurements was performed. The objective function value was used to discriminate between nested models. Standard errors obtained from NONMEM ® and the confidence intervals of the bootstrap analysis in PsN ® (n = 1000) were used to evaluate the precision of the parameter estimates.
For the covariate analysis η-values (ETA) and conditional weighted residuals (CWRES) were used. ETA values are defined as the random effect describing the deviation of the individual empirical Bayes estimate of the parameter from the typical population parameter estimate of the corresponding parameter, e.g. CL or central volume of distribution (V), for a given subject. Using the base model, plots of eta values on CL and V versus covariates were used to investigate potential covariates. In a next step potential covariates were evaluated using a stepwise covariate modelling (SCM) procedure (Lindbom et al., 2005). A significance level of p ≤ 0.05 was used for the forward inclusion and a significance level of ≤ 0.01 for the backward elimination.

Model evaluation
Key models as well as the final model were evaluated using goodness-of-fit plots and normalized prediction distribution errors (npde) (Comets et al., 2008) based on 1000 simulations of the model and a bootstrap analysis based on 1000 samples of the data.

External model validation
Phenobarbital plasma concentrations from study 2 were predicted by fixing the parameters to the parameter estimates in the final model using post hoc Bayesian forecasting. To compare predicted and observed values, goodness-of fit was assessed and bias (mean prediction error [MPE]) was calculated using Eq. (1): where n denotes number of observations. Predicted parameter values were population predicted values for each individual.

Evaluation and optimization of dosing regimen
Based on the validated model , the current dosing regimen was critically reviewed for its ability to reach a target concentration between 15 and 40 mg/L in all children. To this end, 1000 simulations of a child with a birthweight of 0.6, 1, 2 and 3 kg were performed upon an intravenous loading dose of 20 mg/kg and an intravenous maintenance dose of 5 mg/kg daily over 60 days. The weight of the children born weighing 0.6, 1 and 2 kg was simulated to develop according to the preterm growth curves published by Anchieta et al. (2003) and the weight of child born weighing 3 kg was assumed to develop as described by Thulier (2016).
In a next step the dosing regimen was modified in order to improve target attainment in the simulated individuals.

Population pharmacokinetic model
Modelling of plasma concentrations resulted in a one compartment model with intra-individual variability (IIV) on CL and V. A proportional error yielded the best description of the residual variability. Using this model without covariates, the plasma concentrations of the smallest children were underpredicted and concentrations of the heaviest children were overpredicted.
In the covariate analysis, birthweight, actual bodyweight, GA and height showed high correlation coefficients with CL, while only weak or no correlations could be seen for Apgar scores, serum creatinine, umbilical artery pH and PNA. V was highly correlated with birthweight, actual bodyweight, GA and height. Other covariates showed weak correlations. The stepwise covariate modelling proposed birthweight or GA and PNA as predictors for CL and actual bodyweight as predictor of V. No other covariates were found. Model fit was comparable for the combination of PNA and GA or PNA and birthweight on CL. In line with the literature (De Cock et al., 2014) the model containing birthweight and PNA on CL (Fig. 1) was carried forward as dosing based on birthweight might lead to a more comprehensible dosing guideline. The ETAs of the final model showed no correlation to parameters relevant for maturation and development indicating that covariates have been correctly implemented (Supplement 1). Covariate inclusion reduced IIV on CL and V from 56% to 30% and from 59% to 24%, respectively.
The final parameter estimates as well as their respective bootstrap estimate and confidence interval (bootstrap convergence rate = 97.3%) are displayed in Table 2. Bootstrap estimates are in good agreement with the final NONMEM ® parameter estimates and show narrow confidence intervals for all structural parameters.

Model evaluation (study 1)
The model was able to sufficiently describe the data of study 1 that were used to build the model (Fig. 2, left panels). The npde-analysis showed no trends towards a model misspecification (Fig. 3).

External model validation (study 2)
Model validation with the external prospective dataset (study 2) shows that the model can adequately predict the data of the included infants without bias (Fig. 2, right panels). Only the peak concentrations of one individual were underpredicted by the model (Fig. 2e, diamonds). The MPE was − 8.4% corresponding to an adequate description of the data.

Evaluation and optimization of dosing regimen
According to the final model, the current dosing regimen consisting of an intravenous loading dose of 20 mg/kg and an intravenous maintenance dose of 5 mg/kg/day resulted in distinct results in children with different birthweights (Fig. 4a). The figure shows that the mean maximal plasma concentration (C max ) after loading dose is lower in smaller children. This is caused by a higher V [L/kg] in lighter and thereby less mature children, leading to a lower C max following an equal dose per kg bodyweight. After the loading dose the mean plasma concentration of phenobarbital keeps increasing during the first 10 days of life. Thereafter a quasi-steady state concentration is reached for about 5 days after which concentrations decrease over time due to the model predicted maturation of CL with PNA ( Fig. 1). At the same time the percentage of simulations below the threshold of 15 mg/L increases. Lower initial C max values in the smallest newborn result in the highest proportion of subtherapeutic levels at the end of the observed period. The current dosing regimen did not result in any concentrations above the 40 mg/L within the 95% confidence interval of 1000 simulations. The highest concentrations were observed at a PNA of 10-15 days in the simulations of a biggest child (birthweight = 3 kg). The lowest concentrations were observed in the smallest child (birthweight = 0.6 kg) after 6-8 weeks.
When increasing the loading dose from 20 mg/kg to 30 mg/kg (Fig. 4b), the therapeutic window is immediately reached for all birthweights. If the maintenance dose is increased from 5 mg/kg/day to 6 mg/kg/d at a PNA of 15 days, the simulated median plasma concentration stays well in the desired target range (Fig. 4c).

Discussion
In this study, we successfully developed a population PK model for phenobarbital in neonates on the basis of existing TDM data and validated it with data from a prospective study. These results are of particular importance as typically studies in (preterm) neonates are complicated due to ethical and practical constraints which results in frequent off-label use. The results of this study show that CL increases with both birthweight and PNA (Fig. 1) while V is determined by actual bodyweight (Table 2).
Besides providing new insights into the PK of phenobarbital in term and preterm newborns this analysis highlights the advantages of datasharing. Only due to the utilization of TDM data from a previous study (study 1) it was possible to develop a population pharmacokinetic model for term and preterm newborns. The data from the DINO-study (study 2) would have been insufficient to build such a model as it only contained data in infants born before 32 weeks of gestation. Covering the whole age range helps at identifying maturational processes and eventually aids dose finding. Furthermore, data sharing provided the possibility of external validation with an independent dataset, however, only in the very preterm age range. Due to the excellent fit of the validations data (Fig. 2) we do believe that our model is applicable for the whole age range.
So far, only few population PK models of phenobarbital in neonates have been published. Yukawa et al. (2011) studied 70 neonates and infants (GA: 24.1-43 weeks, PNA:1-73 days) and found that PNA and total bodyweight were the dominant predictors of phenobarbital PK. Furthermore, the authors observed and accounted for non-linearity in the PK of phenobarbital at concentrations above 50 mg/L (Yukawa et al., 2011). Data in dogs also indicate that auto-induction might be present at very high doses resulting in concentrations outside the therapeutic window (Abramson, 1988). As a possible mechanism leading to non-linearity, auto-induction of phenobarbital metabolism has been discussed in literature, even though it is anticipated to play a minor role in adults (Nelson et al., 1982). Our analysis does not indicate non-linearity, maybe due to the low number of concentrations above 50 mg/L that might have complicated discovering this trend. On the other hand, it might as such not be relevant as non-linearity occurs outside the therapeutic range of 15-40 mg/L. Marsot et al. (2014) modelled phenobarbital data of 48 newborns born at 27-42 weeks of gestation. Their covariate model is based on allometric scaling with an exponent of three quarters on CL and a linear exponent on V. Additional maturation of CL with PNA was not observed. Based on simulations of the first 12 days of life the authors recommend an intravenous loading dose of 20 mg/kg in all neonates and a weight-based maintenance dose. According to their advice, the recommended maintenance dose decreases with increasing bodyweight (Marsot et al., 2014). Our model as well as the model of Yukawa et al. (2011) indicate that CL increases with PNA which would necessitate an increased maintenance dose over the treatment course (Fig. 4). Based on our simulation period of 60 days we propose that the maintenance dose could be increased by 1 mg/kg at a PNA of 15 days to account for this change.
Although information on the maturation of phenobarbital metabolism in neonates is sparse, Boreus et al. (1978) showed that the excretion of conjugated phenobarbital metabolite is significantly reduced in newborns when compared to adults. In addition, phenobarbital is metabolized by Cytochrome P450, mainly through CYP2C9, with minor metabolism by CYP2C19 and CYP2E1 (Pacifici, 2016). Peak concentrations of one individual were slightly underpredicted by the model. The C-reactive protein (CRP) in this child was elevated up to 100 mg/L. Literature suggests that inflammation as represented by high CRP values might lead to a reduced function of Cytochrome P450 (CYP) enzymes (Aitken et al., 2006). Reduced CYP activity might therefore have resulted in elevated (peak) phenobarbital levels.
The CYP enzymes responsible for phenobarbital metabolism may be subject to maturation processes (Bjorkman, 2006) and thereby explain our findings. Ward et al. showed that the ontogeny of CYP2C19 has a profound effect on the weight-normalized CL of orally administered pantoprazole. In their analysis, the CL of pantoprazole, another CYP2C19 substrate, was positively correlated with PNA in 33 neonates aged 1 to 19 days (Ward et al., 2010).
In line with these findings, Heimann et al. found a phenobarbital half-life of 118.6 ± 16.1 h in term newborns (0-28 days) and a halflife of 62.9 h ± 5.2 h in infants (1-12 months) using a classical two stage approach (Heimann and Gladtke, 1977). Oztekin et al. observe rapidly increasing phenobarbital trough levels during the first four days of life in very low birth weight infants (< 1500 g) (Oztekin et al., 2013), which is in agreement with our results (Fig. 4). We can explain these findings because of a very low CL, upon which steady-state is reached after a long treatment period. Thus, instead of CL, V mainly determines the PK profile of phenobarbital during the first week of life.
In view of these observations, we found high volumes of distribution per kg bodyweight in newborns and thereby a high percentage of initial plasma concentrations below the therapeutic range (Fig. 4). Our model thus implies that neonates require higher loading doses than the currently applied 20 mg/kg. Donn et al. investigated a higher loading dose of 30 mg/kg in asphyxiated term newborns . The mean observed plasma concentration 2 h after loading dose was 30.0 ± 3.2 mg/L and was not associated with an increased number of side effects. Safety data on higher loading doses in preterm newborns is however lacking. Therefore, safety aspects should be evaluated when administering higher loading doses in this vulnerable population.
Phenobarbital is the most frequently used treatment modality in term newborns following asphyxia and whole body hypothermia (Pokorna et al., 2015). In a recent evaluation Pokorna et al. (2015) showed that the combination of asphyxia and hypothermia might have an effect on clearance of phenobarbital (personal communication). Therefore, it might also be desirable to develop models for phenobarbital that can lead to optimized dosing in these conditions.
The results of the presented analysis indicate that the use of the currently applied bodyweight based loading and maintenance dose of phenobarbital in preterm newborns may have its limitations. Recent observations from the population pharmacokinetic field show that bodyweight is often not the only determinant of the PK Calvier et al., 2016). Our model as well as the model of Yukawa et al. (2011) predict that maturation of CL is not related to bodyweight alone but also to PNA. Therefore, preterm infants might require adjusted dosing regimen correcting for this maturation.
The therapeutic window for phenobarbital remains a topic of discussion. It might be possible that the therapeutic window in preterm newborns differs from that in term newborns, children or adults. To be able to detect a potential age-dependent difference, phenobarbital concentrations in the central nervous system would need to be linked to effect measures such as EEG or seizure frequency. Furthermore, the effect of concomitant medication was not tested in our study. It has, for example, been demonstrated that co-medication with phenytoin or valproate might lead to increased phenobarbital serum concentrations (Pacifici, 2016). Thus, for subjects receiving these drugs other doses than the ones recommended here might be suitable.
In conclusion, we strongly encourage data-sharing in the pediatric and neonatal PK setting in order to use the data in an optimal fashion and subsequently to minimize the burden on newborns. We show that data-sharing can help to successfully develop and validate population pharmacokinetic models in neonates. From the results it seems that for phenobarbital both PNA and bodyweight are required to guide dosing of phenobarbital in term and preterm neonates. Based on the presented analysis the loading dose should be adjusted to 30 mg/kg in order to immediately reach the therapeutic window and maintenance dose should be increased by 1 mg/kg/day due at a PNA of 15 days due to a maturation of CL with PNA.

Acknowledgments
The DINO study and all accompanying research were funded by the Netherlands Organisation for Health Research and Development ZonMw (Grant number: 80-83600-98-10190). All authors have no conflicts of interest to declare. Paula Pokorna is supported by the General University Hospital RV-project (64-165/2012) and an unrestricted research grant of the Intensive Care of the ErasmusMC-Sophia Childrens Hospital. Fig. 4. Phenobarbital concentration time profiles (n = 1000) with blue line representing median and grey shaded areas 95% confidence intervals for children with a birthweight of 0.6 kg, 1 kg, 2 kg and 3 kg (dashed red lines = therapeutic window) according to different dosing regimens. a: current dosing regimen (i.v. loading dose: 20 mg/kg; i.v. maintenance dose: 5 mg/kg/day); b: increased i.v. loading dose (30 mg/mg), maintenance dose: 5 mg/kg/day; c: increased i.v. loading dose (30 mg/kg) and increased i.v. maintenance dose (6 mg/kg/day vs. 5 mg/kg/day) starting at a PNA of 15 days. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)