## ABSTRACT

The pharmacokinetic (PK) properties of intravenous (i.v.) benzylpenicillin in term neonates undergoing moderate hypothermia after perinatal asphyxia were evaluated, as they have been unknown until now. A system-specific modeling approach was applied, in which our recently developed covariate model describing developmental and temperature-induced changes in amoxicillin clearance (CL) in the same patient study population was incorporated into a population PK model of benzylpenicillin with *a priori* birthweight (BW)-based allometric scaling. Pediatric population covariate models describing the developmental changes in drug elimination may constitute system-specific information and may therefore be incorporated into PK models of drugs cleared through the same pathway. The performance of this system-specific model was compared to that of a reference model. Furthermore, Monte-Carlo simulations were performed to evaluate the optimal dose. The system-specific model performed as well as the reference model. Significant correlations were found between CL and postnatal age (PNA), gestational age (GA), body temperature (TEMP), urine output (UO; system-specific model), and multiorgan failure (reference model). For a typical patient with a GA of 40 weeks, BW of 3,000 g, PNA of 2 days (TEMP, 33.5°C), and normal UO (2 ml/kg/h), benzylpenicillin CL was 0.48 liter/h (interindividual variability [IIV] of 49%) and the volume of distribution of the central compartment was 0.62 liter/kg (IIV of 53%) in the system-specific model. Based on simulations, we advise a benzylpenicillin i.v. dose regimen of 75,000 IU/kg/day every 8 h (q8h), 150,000 IU/kg/day q8h, and 200,000 IU/kg/day q6h for patients with GAs of 36 to 37 weeks, 38 to 41 weeks, and ≥42 weeks, respectively. The system-specific model may be used for other drugs cleared through the same pathway accelerating model development.

## INTRODUCTION

Therapeutic whole-body hypothermia for (near) term neonates suffering from hypoxic-ischemic encephalopathy (HIE) after perinatal asphyxia has been the standard of care in The Netherlands and Belgium since 2008 (1–9). Several drugs, such as the antibiotic benzylpenicillin, are administered to these patients. Although the effect of hypothermia on the pharmacokinetic (PK) properties of some drugs has already been described, this is not the case for benzylpenicillin (10). One of the goals of the PharmaCool Study, a prospective observational cohort study with neonates undergoing therapeutic whole-body hypothermia, is to evaluate and describe these properties (11). To do so, developing population PK models based on clinical data and drug concentrations is imperative. Subsequently, with these PK models, existing dosage regimens can be evaluated and, if necessary, optimized. However, as performing a systematic PK analysis for each new or existing drug is time-consuming, other approaches are welcomed (12). Pediatric population covariate models describing the developmental changes in drug elimination pathways may constitute system-specific rather than drug-specific information and may be incorporated into PK models of drugs cleared through the same pathway (13–15). This system-specific modeling approach can accelerate the modeling process and is helpful if data are limited.

Previously we have developed and validated a PK model for intravenous (i.v.) amoxicillin in cooled neonates (*n* = 125) (16). Amoxicillin and benzylpenicillin have similar physicochemical properties and in neonates are almost entirely eliminated through glomerular filtration (GF), rather than tubular secretion (17, 18). The objective of this study was to incorporate the covariate model of amoxicillin into a population PK model of benzylpenicillin in the same patient study population (i.e., system-specific model) (16). To evaluate the performance of this system-specific model, a reference PK model of benzylpenicillin in these patients was also built and the two models were compared. With the developed system-specific model, the current dosage regimens of benzylpenicillin have been evaluated and, if needed, optimized.

## RESULTS

Patients.In Table 1 the demographics, clinical characteristics, and samples drawn during the study period are presented for the patients receiving benzylpenicillin. The amoxicillin patient population was similar to the current benzylpenicillin patient population with respect to demographics (16). Ten (24%) patients died during the study period due to withdrawal of medical intensive care treatment. Twenty-two samples (6%) were below the limit of quantification and were not included in the data analysis (19). In Fig. S1 in the supplemental material the observed plasma concentrations versus the time after dose are shown. No serious adverse drug events were recorded.

Pharmacokinetic model building.Figure 1 summarizes all the model-building steps. For the structural model, a two-compartment model parameterized in terms of clearance (CL), volume of distribution (*V*) of the central (*V*_{c}) and peripheral (*V*_{p}) compartments, and intercompartmental clearance (*Q*) best described the benzylpenicillin concentrations in time. Parameters were normalized *a priori* to a body weight of 70 kg using the [3/4] rule (20). Estimation of the allometric exponents resulted in unstable models (nonconvergence in the iterative process), probably due to the small range of weights in the population. An additive-error model described the residual variability of the log-transformed data. Although interoccasion variability (IOV) decreased the objective function value (OFV) significantly, η-shrinkage was 30%, reducing the power of individual predictions (IPRED) and individual conditional weighted residuals (IWRES) plots to detect model and residual model misspecifications, respectively. Also, including covariates to the structural model with IOV increased the additive error. Therefore, IOV was not included in the model. The structural model had an OFV of −75.7.

Then two models were developed: (i) a system-specific model and (ii) a reference model. First, for the system-specific model, the previously developed covariate model for i.v. amoxicillin was directly incorporated into the structural model for benzylpenicillin (16). The amoxicillin covariate model is considered system specific with respect to the effects of postnatal age (PNA), gestational age (GA), core body temperature (TEMP), and urine output (UO) on drug clearance. These system-specific associations were incorporated into the benzylpenicillin model, whereas drug-specific population values for CL, *V*_{c}, *Q*, and *V*_{p} were estimated. The equations for the system-specific model and the drug-specific final parameter estimates can be found in the footnotes to Table 2. For CL and *V*_{c}, interindividual variability (IIV) could be estimated and was correlated (*r* = 0.41). The OFV was −278.3. The η- and ε-shrinkages were ≤20% (21). Second, the reference model was developed by performing a systematic covariate analysis. In the univariate covariate analysis, PNA, TEMP, UO, GA, and multiorgan failure (MOF) were identified as significant covariates for CL, with ΔOFV of −183.9, −147.6, −20.0, −9.3, and −6.5 points, respectively. Associations between aspartate aminotransferase (ASAT), alanine aminotransferase (ALAT), serum creatinine (SCr), sex, urea, and Thompson score versus the PK parameters were not significant. Addition of maturation models describing the association between both PNA and postmenstrual age (PMA) and CL lead to unstable models (nonconvergence in the iterative process), probably due to the small range of PNAs in the population. SCr and ALAT were significant covariates for *V*_{c}, with ΔOFV of −66.8 and −4.8, respectively. After the forward inclusion and backward deletion, UO on CL and all the covariates on *V*_{c} were discarded, while the other above-mentioned significant covariates were retained. The final parameter estimates of the reference model and the equation describing the relationship between the covariates and CL are given in Table 2. For CL and *V*_{c}, IIV could be estimated and was correlated (*r* = 0.82). The OVF was −305.3. The η- and ε-shrinkages were ≤21% (21).

Based on the system-specific population model, CL increased from 0.48 liter/h at PNA of 2 days (i.e., TEMP, 33.5°C) to 0.75 liter/h at 5 days PNA (i.e., TEMP, 37.0°C) for a typical patient (GA, 40 weeks; birthweight [BW], 3,000 g; UO, normal [2 ml/kg/h]).

Model evaluation and validation.Table 2 shows the results of the bootstrap analysis of the final system-specific model and final reference model. The bootstrap results were in agreement with those of the respective final models, and therefore, the final model parameter estimates were considered reliable. The goodness-of-fit (GOF) plots of both models are shown in Fig. 2. The plots of population predictions (PRED) and IPRED versus observed benzylpenicillin concentrations (observations) show uniformly distributed points on both sides of the lines of identity, without substantial and systematic deviations from the lines of identity. Differences between PRED and observations are caused by IIV, while differences between IPRED and observations are caused by residual variability. The conditional weighted residuals (CWRES) versus time after dose and PRED show that most points fall between −2 and 2 and are distributed uniformly around 0. These points indicate that there is no major bias in the population component of the final model and that an appropriate structural model is found for most individuals. The normalized prediction distribution error (NPDE) plots of both models are shown in Fig. 3. The quantile-quantile plot compares the distribution of the NPDE with a normal *N(0*, *1)* distribution. The points are close to the line *y = x*, which suggests the absence of (major) bias. Furthermore, in this figure the histogram of the NPDE is shown, which also shows that the NPDE follows the normal distribution, indicating a good fit of the model to the individual data. No ill conditioning was found in the system-specific model (eigenvalue [EV] of 48) or reference model (EV of 80). The visual predictive checks (VPCs; see Fig. S2 in the supplemental material) indicate that the final models describe the observed data well.

Model comparison.The reference model had a 27-point-lower OFV than the system-specific model at a 4-point difference in degrees of freedom (−305.3 and −278.3, respectively) and was therefore statistically better in describing the benzylpenicillin data. A 4-point difference in degrees of freedom was chosen because the system-specific model had 4 fixed covariates and the reference model had 4 estimated covariates. Visual examination of the GOF plots (Fig. 2) indicates that the system-specific and the reference models describe the observed concentrations equally well and that the difference in performance of the two is negligible. In Fig. S3 in the supplemental material the individual *post hoc* CL (CL_{i}) and population predicted CL (TVCL) versus the most predictive covariate (PNA) are shown for both models. Similar CL_{i}s and TVCLs are estimated by the system-specific and reference models. Furthermore, the two models estimate fixed and random effects with acceptable precision (Table 2; coefficients of variation of <50%) and perform equally well in terms of predictive performance (Fig. 2 and 3).

Simulations.Simulations of various dose regimens were performed using the system-specific model (Table 3; see also Table S1 in the supplemental material) Dosage regimens were evaluated for patients with GAs of 36 to 42 weeks stratified by UO. For the simulations, the time that the non-protein-bound plasma benzylpenicillin concentration exceeds the MIC was taken as a pharmacodynamic target. With an MIC of 0.125 mg/liter, all the simulated dosage regimens suffice in attaining a non-protein-bound plasma concentration exceeding the MIC for at least 50% of the dosing interval for minimally 90% of simulated patients, except for patients with a GA of 42 weeks and UO of 2 or 6 ml/kg/h receiving 75,000 IU/kg/day every 8 h (q8h). However, with an MIC of 1 mg/liter, the target was not met in patients with a GA of 38 weeks and UO of ≥2 ml/kg/h and GA of >38 weeks both receiving 75,000 IU/kg/day q8h, nor was it met in patients with GAs of 42 weeks and UO of ≥2 ml/kg/h receiving 150,000 IU/kg/day q8h. On the other hand, dosages of 150,000 IU/kg/day q8h for patients with GAs of 36 weeks and 200,000 IU/kg/day q6h for patients with GAs of 36 or 38 weeks produced peak concentrations that were too high (>100 mg/liter) (Table S1). In the simulations, a dose regimen was acceptable, if the maximum peak concentration was less than 100 mg/liter, as higher concentrations are potentially associated with toxicity (22). Based on these results, neonates undergoing whole-body cooling after HIE and born with GAs of 36 to 37 weeks, 38 to 41 weeks, and ≥42 weeks should receive 75,000 IU/kg q8h, 150,000 IU/kg/day q8h, and 200,000 IU/kg/day q6h, respectively.

## DISCUSSION

Recently, we have developed and validated a population PK model of amoxicillin administered intravenously to neonates undergoing therapeutic hypothermia for HIE (16). In this study, we have successfully applied a system-specific modeling approach to describe the PK of benzylpenicillin in the same patient study population by incorporating the amoxicillin covariate model in a newly developed population PK model with *a priori* BW-based allometric scaling for benzylpenicillin.

The amoxicillin covariate model describes the time dependency of CL on the basis of associations with PNA, GA, TEMP, and UO (16). UO is an indicator for possible renal damage following resuscitation, and its assessment together with other laboratory parameters (i.e., SCr, blood urea nitrogen, and electrolytes) can give an estimation of true renal function. In the system-specific model, polyuria (UO of 6 ml/kg/h) is associated with a 22% increase of CL compared to that with oliguria (UO of 0.5 ml/kg/h). This can be explained by the fact that in neonates tubular secretion is inefficient and glomerular filtration is believed to be the major pathway for elimination (18, 23). Furthermore, in the system-specific model, CL is higher in patients with higher GAs and BWs, which are known predictors of glomerular filtration rate (GFR) maturation and drug CL (18, 24, 25). Finally, separate increases of 27% and 23% in CL are seen if the TEMP increases from 33.5 to 37.0°C and PNA increases from 2 to 5 days, respectively. Although no significant difference in the incidence of renal impairment during cooling compared with standard care of asphyxiated neonates in neonatal intensive care units (NICUs) was seen in a meta-analysis of several randomized controlled trials, animal studies have shown that mild hypothermia may reduce the GFR in the immature kidney (7, 26). Furthermore, we have previously shown that in neonates, moderate hypothermia delayed normalization of CL of gentamicin, an antibiotic also predominantly eliminated by glomerular filtration (27).

The system-specific modeling approach allows for the extrapolation of PK properties between drugs. Amoxicillin and benzylpenicillin exhibit similar PK properties in neonates, as GF, rather than tubular secretion, is the main route of elimination of both drugs (17, 18, 23). Furthermore, their physicochemical properties are alike, with molecular masses of 365.4 (amoxicillin) and 334.4 (benzylpenicillin) g/mol and acid dissociation constants of 2.4, 7.4, and 9.6 (amoxicillin) and 2.74 (benzylpenicillin).

In the present study, a reference model was developed with a systematic covariate analysis based on the benzylpenicillin data set (*n* = 41). Significant covariates on CL were PNA, GA, TEMP, and MOF. The system-specific model differs from the reference model only by including UO as a significant covariate on CL instead of MOF. MOF was defined as renal and/or liver failure on top of the existing encephalopathy due to perinatal asphyxia and thus comprises ASAT, ALAT, SCr, and UO. During the multivariate covariate analyses, UO was more significant than MOF for the amoxicillin PK model and MOF was not retained. For the reference model, MOF was most significant and UO was not retained.

The descriptive and predictive performances of the system-specific model were equal to those of the reference model. Consequently, the covariate model of amoxicillin contains system-specific information on the developmental changes in GF and the effect of TEMP on this process. Krekels et al. have also successfully extrapolated the covariate model for the glucuronidation of morphine in (pre)term neonates to children up to 3 years old to describe the developmental changes in the glucuronidation of zidovudine in term neonates and infants (15). De Cock et al. demonstrated that models using an amikacin covariate model performed similarly to independent reference models of other aminoglycosides and vancomycin (14).

Based on the system-specific model, benzylpenicillin CL was 0.48 liter/h (IIV, 49%) for a typical patient (GA, 40 weeks; BW, 3,000 g; PNA, 2 days [TEMP, 33.5°C]; UO, normal [2 ml/kg/h]), rising to 0.75 liter/h at a PNA of 5 days (TEMP, 37°C). The corresponding CL values expressed per kg of weight (linear) are 0.16 and 0.25 liter/h/kg, respectively. For the same typical patient, *V*_{c} was 0.62 liter/kg (IIV of 53%). There is no other literature on the PK of benzylpenicillin in cooled neonates; such information is available only for (pre)term noncooled patients. Muller et al. found a CL of 0.10 liter/h at a PNA of 3 days (GA, <32 weeks). The *V* at steady state (*V*_{ss}) was 0.54 liter (28). In another study with two groups of very-low-BW (VLBW) neonates, the median CLs were 1.2 ml/min/kg and 1.5 ml/min/kg, respectively, and the median *V*_{ss}s were 0.41 and 0.64 liter/kg, respectively (29). These studies were conducted with preterm and VLBW neonates, explaining the lower CL than in our study. This is not surprising, as PNA and GA are known predictors of GFR maturation and drug CL (24). In neonates, tubular secretion is very inefficient and GF (also low initially after birth) is believed to be the major pathway for elimination (17, 18, 23). However, our results are similar to the findings of two other studies with neonates born after GAs of 27 to 40 weeks (mean CL of 0.12 liter/h/kg and a *V*_{c} of 0.61 liter/kg) (30) and with patients with PNAs of <8 days and BWs of >2,000 g (CL of 3.12 liters/h/1.73 m^{2}) (17), probably due to the more comparable GAs.

The lack of compliance to the Dutch Pediatric Formulary (31) was striking (Table 4); in general, higher dosage regimens were used in the studied population. Simulations were performed on basis of the system-specific model to evaluate both the dosage regimens applied in this study and those recommend by the Dutch Pediatric Formulary. Doses were simulated for patients with GAs ranging from 36 to 42 weeks stratified by UO. Noteworthy was that dosages of ≥150,000 IU/kg/day led to high (>100 mg/liter) peak concentrations in patients with GAs of ≤38 weeks. Neurologic adverse events (including seizure, encephalopathy, crystalluria, and tremor) have been reported during treatment with penicillins due to excessive accumulation. (Pre)term neonates are an at-risk population due to immature renal function resulting in higher blood concentrations. As HIE represents the most common cause of both seizures and strokelike events in the term infant, the additional role of epileptogenic drugs, such as penicillins, may be underestimated (32–37). Based on our results, we advise 75,000 IU/kg/day q8h, 150,000 IU/kg/day q8h, and 200,000 IU/kg/day q6h for patients with GAs of <38 weeks, 38 to <42 weeks, and ≥42 weeks, respectively, regardless of variation in the significant covariates on CL (GA, PNA, TEMP, and UO). This is due to the wide therapeutic range of benzylpenicillin.

As there might be some concern regarding the validity of extrapolating a covariate model of one drug to another, the following requirements were met in this study. First, the amoxicillin covariate model was built using a large amount of data (125 patients and 1,280 measured plasma concentrations) and from patients with the same clinical characteristics and disease status (i.e., asphyxiated newborns undergoing moderate hypothermia admitted to the NICU). Second, the covariate model was internally and externally validated (16). Furthermore, the system-specific model performed as well as the reference model (see “Model comparison” in Results), which confirms the applicability of this approach.

Our study had some limitations. First, we did not measure the free concentration of benzylpenicillin, nor did we evaluate the effect of hypothermia on protein binding (PB) of benzylpenicillin. Based on literature for adults, we have assumed a PB of 50% in our patient population (38–40).

Second, we did not include patients with GAs of <36 weeks in this study, as at the time of conducting this research preterm neonates were not treated with moderate hypothermia in The Netherlands or Belgium (11).

In conclusion, this is the first prospective study evaluating the PK characteristics of benzylpenicillin in a cohort of (near) term neonates undergoing moderate hypothermia for HIE due to perinatal asphyxia. We have shown that a system-specific population model, in which an i.v. amoxicillin covariate model was incorporated, performed similarly well in describing and predicting the PK properties of benzylpenicillin in (near) term neonates undergoing moderate hypothermia. This suggests that the covariate model contains information on the developmental changes in CL and the effect of TEMP on this process. In the future, the system-specific covariate model can be used for other drugs cleared through the same pathway, accelerating model development and optimizing sparse data analysis.

## MATERIALS AND METHODS

Patients.Data were collected from 10 Dutch and 2 Belgian NICUs participating in a multicenter prospective observational cohort study (PharmaCool study, November 2010 until October 2014) (11). The current analysis is based on data from 41 patients receiving benzylpenicillin during moderate hypothermia out of a total of 187 newborns included in the PharmaCool study (11). According to national protocol, (near) term newborns meeting the criteria of perinatal asphyxia and ensuing moderate to severe encephalopathy who were at a PNA of <6 h were cooled to a TEMP of 33.5°C for 72 h and then rewarmed (0.4°C/h) to normothermia (37.0°C) (11). All newborns were eligible for inclusion if undergoing moderate hypothermia. Exclusion criteria were the absence of parental informed consent, serious congenital abnormalities or malformations, or absence of central venous or arterial access. The study was approved by the Institutional Review Board of each participating center.

Data and sample collection.Data on GA, BW, sex, Thompson score (a numeric scoring system for the assessment of hypoxic ischemic encephalopathy during the neonatal period) (41), comedication, mean daily UO (in milliliters per kilogram per hour; collected via indwelling urine catheter), SCr, urea, ASAT, and ALAT were collected. Benzylpenicillin was administered as an i.v. bolus in dosages according to local practice. Dosing information and sampling times were recorded in a digital case report form. The dosage regimens are shown in Table 4. Local dosing guidelines could vary among study centers with respect to total daily dose and dosing frequency. This explains the ranges of total daily dosages in our study population. As most patients were born in referring regional hospitals, the first benzylpenicillin dose was largely unknown and in those cases a standard dose of 25,000 IU/kg was presumed.

For the PK analysis, 6, 3, and 6 blood samples were drawn on days 2, 3 (hypothermia), and day 5 (normothermia) of life, respectively. If available, residual blood samples were also collected on day 4 (rewarming). The benzylpenicillin concentration was analyzed in these samples by zwitterion hydrophilic interaction chromatography coupled to tandem mass spectrometry in only 10 μl of plasma. The range of the method was 0.5 to 40 mg/liter. Accuracy and imprecision at the lowest, middle, and upper limits of quantification were 102% and 14%, 103% and 7.7%, and 103% and 7.5%, respectively. Samples containing >40 mg/liter of benzylpenicillin were diluted 10 times and reanalyzed. As concentrations were measured in milligrams per liter, and dosages were given in international units, a multiplication factor of 1,670 was applied for conversion from international units to milligrams for the PK analysis (Sandoz BV, The Netherlands, personal communication).

Pharmacokinetic modeling.Two benzylpenicillin population PK models were developed. First, a system-specific model in which an extensively validated covariate model for i.v. amoxicillin, built from data of 125 neonates enrolled in the PharmaCool study, was directly incorporated (16). Second, a reference model was developed by performing a systematic pharmacokinetic covariate analysis based on data of the 41 patients receiving benzylpenicillin. The performances of both models were evaluated and then compared with each other.

Model development.Data were analyzed using the first-order conditional estimation (FOCE) method with interaction option in the nonlinear mixed-effects modeling software NONMEM version 7.2 (Globomax LLC, Hanover, MD). Tools like R (https://www.r-project.org/; open-source, S-based statistical software, version 0.98.945), XPose (42), and PsN (43) were used to visualize and evaluate the models. The model building process was performed in a stepwise fashion: (i) choice of structural model, (ii) choice of error model, (iii) choice of covariate model, and (iv) model evaluation. The first two steps were equal for the system-specific model and reference model.

Structural model.First, one-, two-, and three-compartment models were fitted to the log-transformed data. PK parameters were estimated for CL, *V*_{c}, *V*_{p}, and *Q*. IIV and covariance were tested assuming a log-normal distribution for all parameters (44). BW was included *a priori* in all model parameters using an allometric power model in which the parameter values were standardized to a body weight of 70 kg to account for variability in PK parameters due to the various sizes of individual children. The allometric scaling parameter (PWR) was fixed at values of 0.75 and 1.0 for CL and *V*, respectively (20). PWR was also estimated to evaluate if this would result in a better fit. Residual variability was estimated using an additive-error model on the log-transformed data (i.e., a proportional error on a linear scale). IOV was estimated for the log-transformed data with an occasion being a 24-h period. The likelihood ratio test was used to evaluate statistical significance between nested models where a reduction in OFV of ≥3.9 points was considered statistically significant (*P* < 0.05 based on χ^{2} distribution; df = 1). Also, GOF plots (PRED or IPRED versus observations, conditional weighted residuals [CWRES] versus time, and individual CWRES [IWRES] versus IPRED), the total number of parameters, visual improvement of individual plots, correlation matrix, confidence intervals of parameter estimates, ill conditioning, and η- and ε-shrinkages were assessed (21, 44).

Covariate model.The covariate model for i.v.-administered amoxicillin in neonates undergoing moderate hypothermia was directly incorporated into the developed structural model for benzylpenicillin (16). The covariate model is depicted below:
_{i} is CL predicted for individual i and θ_{1} is the typical value of CL for the population. The covariates were normalized to their median values in the amoxicillin database. Exponents of the power functions were fixed to their estimated values. The covariate model describes the developmental changes in CL through GF, whereas the population value of CL (θ_{1}) is still estimated by NONMEM since it is considered a drug-specific property.

For the reference model, a systematic covariate analysis was performed in which covariates were tested for significance. Continuous covariates were tested using a power function equation in which the covariates were normalized to their median value in the database (44). Continuous covariates tested were PNA, PMA, GA, TEMP, ASAT, ALAT, SCr, urea, Thompson score, and UO. PMA was calculated by adding GA (in days) to PNA (in days).

Categorical covariates were implemented in the model according to the following equation:
_{1} is the typical value of the parameter P for the population and θ_{2} is the fractional difference in P between categories. For categorical dichotomous data (sex, MOF [yes/no], and inotropic comedication [yes/no]) the value of the covariate was set to 0 for the reference classification and 1 for the other classification. MOF was considered to be present if a patient had renal and/or liver function failure on top of the existing encephalopathy due to perinatal asphyxia, with renal failure being either SCr of >125 μmol/liter or anuria or oliguria (UO <1 ml/kg/h) during the first 24 h after birth. Liver failure was defined as an increase in ASAT or ALAT of >100 U/liter (45).

A maturation model (sigmoid Emax maturation function [Fmat]) was tested to evaluate the effect of maturation of organ function on the parameter estimates (20). This function is depicted below:

Covariates were tested through forward addition/backward deletion with *P* values of <0.05 (a decrease in OFV of at least 3.9 points; df = 1) and <0.001 (a decrease in the OFV of at least 10.3 points; df = 1), respectively (44).

Model evaluation and model comparison.Parameter precision and model stability for both the system-specific model and the reference model were evaluated by performing nonstratified nonparametric bootstrap analyses using the PsN Toolkit (43) in which the benzylpenicillin data set was resampled 1,000 times to produce a new data set with the size of the original but with a different combination of individuals. The parameter estimates obtained with the bootstraps (median values and 95% confidence intervals [CI]) were compared to the parameter estimates of the referring model.

The NPDE method was used to evaluate the predictive properties of both models by simulating the benzylpenicillin data set 2,000 times in NONMEM using Monte-Carlo simulations in which the random effects were included. If the model adequately describes the simulated data, the NPDE is expected to follow a *N(0*,*1)* distribution. A histogram and quantile-quantile plot of the NPDE distribution were used to evaluate both models (46). VPCs were performed to evaluate how well the models describe the central tendency and variability of the observed data (simulation properties). The simulated median and 95% prediction intervals, derived from 1,000 data sets simulated from the developed model using the original data set as a template, were compared to the ones from the observed data. The simulated and observed medians and 95th percentiles correspond closely if the model adequately describes the data.

Several diagnostics were evaluated to compare the performances of the system-specific and reference models (14, 15): (i) although not nested, both models are based on the same patient study population and data, and therefore, the −2 log likelihoods (OFV) of both models were used to statistically compare their descriptions of the benzylpenicillin data; (ii) the GOF plots were compared to visually evaluate the descriptive performance; (iii) the individual and population predicted benzylpenicillin CL values simulated by both models were plotted against the most predictive covariate (PNA) to evaluate whether the individual predicted parameters were equally distributed around the population parameters; (iv) the results of the bootstrap analysis as well as ill conditioning and shrinkage of both models were assessed; and (v) the predictive performance of the models was evaluated by comparison of the NPDE results.

Simulations.For simulations, the system-specific model was chosen above the reference model, as the covariate model of the former was based on more patients. Monte-Carlo simulations (*n* = 1,000) were performed in NONMEM to examine the optimal dosing regimen on a selection of patients from the original benzylpenicillin database (GAs of 36 weeks [BW, 2,500 g], 38 weeks [BW, 3,130 g], 40 weeks [BW, 3,500 g], and 42 weeks [BW, 3,800 g]). The i.v. benzylpenicillin dosage regimens most frequently used in the patient population (150,000 IU/kg/day q8h and 200,000 IU/kg/day q6h) and the dosage advised by the Dutch Pediatric Formulary (31) (75,000 IU/kg/day q8h) were simulated. For graphical display, UO was fixed at 0.5 ml/kg/h (oliguria), 2 ml/kg/h (normal UO), or 6 ml/kg/h (polyuria), although in practice UO can vary within patients. During the first 72 h of life, TEMP was set at 33.5°C, after which it increased with 0.4°C/h to normothermia (37°C). The median (including 5th and 95th percentiles) of the simulated trough and peak concentrations were computed.

For penicillins, the time that the non-protein-bound concentration in blood exceeds the MIC (T>MIC) should be at least 40 to 50% in neonates (47). The most important microorganism causing early-onset sepsis is Streptococcus agalactiae (incidence rate, 0.43% [95% CI, 0.37 to 0.49]; MIC, 0.125 mg/liter), while Listeria monocytogenes (MIC, 1.0 mg/liter) should also be considered (48, 49). The plasma PB of benzylpenicillin in neonates is unknown. In adults the PB of benzylpenicillin ranges from 50 to 65% (38, 39). As serum albumin is lower in newborns and gradually increases with postnatal age (40), we used a PB of 50% in our patients in order to calculate the desired non-protein-bound MIC. Therefore, the total benzylpenicillin concentration must be above 2 mg/liter for at least 40 to 50% of the dosing interval to achieve bacteriologic efficacy (38, 39). As a target we have aimed for a non-protein-bound plasma concentration exceeding the MIC for at least 50% of the dosing interval for minimally 90% of simulated patients. Although β-lactam antibiotics are considered safe due to the wide therapeutic range, excessive accumulation can lead to the development of adverse drug events, such as seizures and crystalluria (37). Therefore, we have chosen a maximum peak level of 100 mg/liter (22).

## ACKNOWLEDGMENTS

This study was funded by the Dutch government (ZonMw grant number 40-41500-98-9002).

We thank I. Corten, an employee of the Clinical Research Unit (CRU) Academic Medical Center, for her contributions with regard to the data management.

Y.A.B., T.R.D.H., and R.A.A.M. designed the research, analyzed the data, and wrote the manuscript. D.H.G.M.N., A.H.V.K., J.H.V.D.L., F.G., and R.C.J.D.J. designed and performed the research. P.H.D., A.V.H., K.P.D., H.L.M.V.S., M.R., I.A.Z., F.C., and A.Z. performed the research.

D.H.G.M.N. is the Clinical Research Coordinator for the PharmCool study.

Collaborators of the PharmaCool study group include Mieke J. Brouwer, Department of Neonatology, Wilhelmina Children’s Hospital, University Medical Center Utrecht, Utrecht, The Netherlands; Marcel P. van den Broek, and Carin M. A. Rademaker, Pharmacy Department, University Medical Center Utrecht, Utrecht, The Netherlands; Djien Liem and Katerina Steiner, Department of Neonatology, Radboud University Medical Center, Nijmegen, The Netherlands; Sinno H. P. Simons and Annelies A. Bos, Department of Pediatrics, Division of Neonatology, Erasmus MC-Sophia Children's Hospital, Rotterdam, The Netherlands; S. M. Mulder-de Tollenaer and L. J. M. Groot Jebbink-Akkerman, Department of Neonatology, Isala Clinics, Zwolle, The Netherlands; and Michel Sonnaert and Fleur Anne Camfferman, Department of Neonatology, Vrije Universiteit Brussel, Brussels, Belgium.

## FOOTNOTES

- Received 13 November 2017.
- Returned for modification 9 December 2017.
- Accepted 21 January 2018.
- Accepted manuscript posted online 29 January 2018.
Supplemental material for this article may be found at https://doi.org/10.1128/AAC.02311-17.

- Copyright © 2018 American Society for Microbiology.