Prediction of Drug Exposure in Critically Ill Encephalopathic Neonates Treated With Therapeutic Hypothermia Based on a Pooled Population Pharmacokinetic Analysis of Seven Drugs and Five Metabolites.

Drug dosing in encephalopathic neonates treated with therapeutic hypothermia is challenging; exposure is dependent on body size and maturation but can also be influenced by factors related to disease and treatment. A better understanding of underlying pharmacokinetic principles is essential to guide drug dosing in this population. The prospective multicenter cohort study PharmaCool was designed to investigate the pharmacokinetics of commonly used drugs in neonatal encephalopathy. In the present study, all data obtained in the PharmaCool study were combined to study the structural system specific effects of body size, maturation, recovery of organ function, and temperature on drug clearance using nonlinear mixed effects modeling. Data collected during the first 5 days of life from 192 neonates treated with therapeutic hypothermia were included. An integrated population pharmacokinetic model of seven drugs (morphine, midazolam, lidocaine, phenobarbital, amoxicillin, gentamicin, and benzylpenicillin) and five metabolites (morphine-3-glucuronide, morphine-6-glucuronide, 1-hydroxymidazolam, hydroxymidazolam glucuronide, and monoethylglycylxylidide) was successfully developed based on previously developed models for the individual drugs. For all compounds, body size was related to clearance using allometric relationships and maturation was described with gestational age in a fixed sigmoidal Hill equation. Organ recovery after birth was incorporated using postnatal age. Clearance increased by 1.23%/hours of life (95% confidence interval (CI) 1.03-1.43) and by 0.54%/hours of life (95% CI 0.371-0.750) for high and intermediate clearance compounds, respectively. Therapeutic hypothermia reduced clearance of intermediate clearance compounds only, by 6.83%/°C (95% CI 5.16%/°C-8.34%/°C). This integrated model can be used to facilitate drug dosing and future pharmacokinetic studies in this population.

WHAT IS THE CURRENT KNOWLEDGE ON THE TOPIC?

☑ Pharmacokinetics and thereby drug exposure in encephalopathic neonates is highly variable due to a multitude of patient, disease and treatment related factors. Population pharmacokinetic models and dosing guidelines have been developed for several individual drugs. It can be anticipated that structural system‐specific effects exist that can be applied to pharmacokinetics in general.

WHAT QUESTION DID THIS STUDY ADDRESS?

☑ Can drug exposure in asphyxiated neonates undergoing therapeutic hypothermia be predicted based on current pharmacokinetic knowledge?

WHAT DOES THIS STUDY ADD TO OUR KNOWLEDGE?

☑ Data on 12 drugs and metabolites from the prospective multicenter cohort study PharmaCool were successfully combined in an integrated population pharmacokinetic model. Clearance increased over time with 1.23%/hours of life and 0.54%/hours of life for high‐clearance and intermediate‐clearance drugs, respectively; hypothermia reduced clearance of intermediate‐clearance drugs only by 6.83%/°C.

HOW MIGHT THIS CHANGE CLINICAL PHARMACOLOGY OR TRANSLATIONAL SCIENCE?

☑ This integrated model can be used to predict the impact of illness, body size, maturation, and hypothermia on drug exposure and, thus, facilitates treatment and informs future pharmacokinetic studies in this fragile population.

Dosing of drugs in neonates is challenging due to limited evidence from clinical studies, body size related differences in pharmacokinetics, and ongoing maturation of organ function. For neonates admitted to a neonatal intensive care unit (NICU), dose optimization is even more challenging due to additional variability associated with illness severity and recovery thereof. These neonates are exposed to a multitude of drugs, such as sedatives, analgesics, antibiotics, and anti‐epileptic drugs (AEDs) for which the optimal dose is often unknown.

Perinatal asphyxia is one of the leading causes of neonatal morbidity and mortality worldwide. It can lead to hypoxic‐ischemic organ damage throughout the body, potentially resulting in multi‐organ failure. , Central nervous system dysfunction leading to neonatal encephalopathy (NE) is a major concern as the brain is least likely to recover. Therefore, (near‐)term neonates with moderate or severe NE after perinatal asphyxia are routinely treated with therapeutic hypothermia (TH) to reduce the incidence of death and long‐term developmental disability. , According to the Dutch National Clinical Protocol, TH is intended for neonates with moderate or severe NE, a gestational age (GA) of at least 36 weeks, and has to be started within 6 hours after birth. During TH, the body temperature of the neonate is reduced to 33.5°C for 72 hours, after which the neonate is gradually rewarmed. Appropriate dosing of essential drugs in this vulnerable population is challenging due to the multitude of complicating factors, such as organ failure and recovery, changes in body temperature, infections, and multiple concomitant medications.

To investigate the pharmacokinetics of frequently used drugs in neonates with NE during and after treatment with TH, the prospective multicenter observational cohort PharmaCool study has been conducted in 12 level III NICUs in the Netherlands and Belgium.

The resulting publications from the PharmaCool study group have separately described the pharmacokinetics of antibiotics (amoxicillin, gentamicin, and benzylpenicillin), sedatives (morphine and midazolam), and AEDs (phenobarbital and lidocaine) in this population and population pharmacokinetic models have been developed for each individual drug. , , , , For all drugs, body size related changes were adequately described with standard allometric relationships. However, the effects of TH, maturation, and recovery of organ function after asphyxia varied between the different drugs and metabolites. This can partly be explained by the different number of patients available for each individual compound with sometimes relatively small numbers, albeit that this was the largest study performed thus far in this population.

It can be anticipated that the effects regarding body size, maturation, recovery of organ function, and body temperature transcend individual drugs and reflect underlying physiological processes that can be applied to pharmacokinetics in general. Therefore, we hypothesized that for drugs with similar clinical pharmacokinetic characteristics (e.g., hepatic or renal clearance) the impact of these effects is similar. Amoxicillin, gentamicin, and benzylpenicillin are all renally cleared drugs. , , Morphine is metabolized by glucuronosyltransferase (UGT) 2B7 into morphine‐3‐glucuronide (M3G) and morphine‐6‐glucuronide (M6G). Midazolam is metabolized by cytochrome P 450 (CYP) 3A into 1‐hydroxymidazolam (OHM) and subsequently by UGT into hydroxymidazolam glucuronide (HMG). Phenobarbital undergoes slow hepatic metabolism through several CYP enzymes and lidocaine is metabolized into monoethylglycinexylidide (MEGX), mainly through CYP1A2 and CYP3A. , M3G, M6G, HMG, and MEGX are excreted renally. In the current study, all data available for analysis from the PharmaCool study were integrated to examine the structural system specific effects of body size, maturation, recovery of organ function, and temperature on the pharmacokinetics of commonly used drugs in NE. Our aim is to develop a framework for the prediction of clearance of other drugs in this fragile population based on current pharmacokinetic knowledge, to facilitate treatment, and to guide future pharmacokinetic studies.

METHODS

Setting, study design, and study population

The prospective multicenter observational cohort PharmaCool study (www.trialregister.nl, NTR2529) was conducted in 12 tertiary NICUs in the Netherlands and Belgium. Neonates undergoing TH for NE were eligible for inclusion. Exclusion criteria were GA < 36.0 weeks, severe congenital malformations, encephalopathy due to other causes than perinatal asphyxia, and the absence of central venous or arterial access for noninvasive blood sampling. Written parental informed consent was obtained from each included neonate. The study was approved by the institutional review boards of all participating centers.

Pharmacokinetic data were collected from seven drugs: gentamicin, amoxicillin, benzylpenicillin, morphine, midazolam, phenobarbital, and lidocaine. For none of the drugs, initial dosing or choice of therapy was influenced by the study protocol. Any dose adjustment was based on clinical care or therapeutic drug monitoring according to local clinical protocol. Full dosing information was recorded in the online Case Report Forms of the neonates participating in this study.

Pharmacokinetic sampling and bioanalyses

For the pharmacokinetic analyses, blood samples were obtained on days 2–5 after birth, both during and after TH. All plasma concentrations, including metabolites M3G and M6G (morphine), OHM and HMG (midazolam), and MEGX (lidocaine), were determined using validated liquid chromatography‐tandem mass spectrometry assays. The sampling schedule for each drug as well as details regarding the bioanalyses have been described previously. , , , , , ,

Body size

All physiological processes are related to body size, and body weight (BW) is the most common descriptor for body size. It is well documented that many physiological processes and organ sizes exhibit an allometric relationship with BW, not only for humans but across different species. Because pharmacokinetic parameters are dependent on these physiological processes, parameters, such as clearance and volume of distribution, can also be described using allometric equations relative to BW. , In this study, an exponent on BW of 0.75 for clearance and an exponent of 1 for volume of distribution were used.

Maturation

In neonates, the function of organs responsible for drug clearance is immature and, as a result, pharmacokinetic parameters in neonates differ from older children and adults. Neonatal renal function at birth is underdeveloped compared with adults and undergoes maturation over the first weeks to months of life until it reaches body size adjusted adult values between 8 and 12 months after birth. Hepatic clearance in neonates is also attenuated compared with adults and the enzymes responsible for drug metabolism, such as CYP and UGT, mature at different rates. , , Rhodin et al. showed that maturation of renal clearance across the entire pediatric population was well described using postmenstrual age (PMA) with a sigmoidal Hill equation. The TM50, the PMA at which clearance is 50% of the mature value, was estimated at 55.4 weeks and the Hill coefficient describing the slope of the sigmoidal curve at 3.33. Knøsgaard et al. found that maturation of morphine clearance was also related to PMA with a sigmoidal Hill equation and with similar values for TM50 and Hill coefficient (54.2 weeks and 3.92, respectively).

In our study, PMA was almost fully determined by GA as data was only collected until 5 days after birth. GA was, therefore, introduced in our model as a covariate describing maturation using the following equation:

Our population consisted of neonates with a GA between 36 and 42 weeks, which is at least 12 weeks before reaching the TM50 for both models describing maturation. Both published models on maturation behave similarly at these early time points after birth and, therefore, we choose to fix the Hill coefficient to 3.92 and TM50 to 54.2 weeks.

To test the validity of these assumptions for BW and maturation, the population conditional weighted residuals for each compound were plotted against BW and GA.

Recovery of organ function

In the individual models, a strong increase in clearance after birth was identified, which is much larger than can be explained by maturation or the influence of TH. During the hypoxic‐ischemic event, both the kidneys and liver are deprived of oxygen, resulting in possible functional nephron and hepatocyte damage. , , , , Additionally, cardiac output might also be hampered due to a loss in myocardial function. After resuscitation and stabilization at the NICU, gradual recovery of organ function after birth was anticipated. In our model, this was described using postnatal age (PNA) as a separate covariate on clearance. The effect was tested separately for the drugs and metabolites fully renally cleared (gentamicin, amoxicillin, benzylpenicillin, M3G, M6G, HMG, and MEGX) vs. hepatically cleared drugs and metabolites (morphine, phenobarbital, midazolam, OHM, and lidocaine). Subsequently, it was tested whether differences in effect could be identified for hepatically cleared drugs with high clearance (lidocaine) vs. intermediate (morphine, midazolam, and OHM) and vs. low clearance drugs (phenobarbital). Because we hypothesized that recovery of organ function could differ between neonates due to differences in asphyxia severity, interindividual variability of the effect of PNA on clearance was also included.

Body temperature

Alterations in body temperature (TEMP) may be of influence for drug clearance. TH decreases heart rate and cardiac output, which will subsequently reduce kidney and liver perfusion. , Additionally, hepatic clearance might also be affected by altered activity of liver enzymes. Because most enzymatic processes exhibit temperature dependency, a lower body temperature might result in reduced enzyme activity and, thus, reduced clearance. , Upon rewarming, processes hampered by the hypothermic state can recover, resulting in an increase in drug clearance.

TEMP was tested as a continuous variable using a dynamic model of temperature over time. For all neonates, the reported start and end times of TH were used to determine the period of TH treatment. TEMP during TH was set at 33.5°C, with consecutive rewarming at 0.4°C/hour (i.e., rewarming time 7.5 hours) until 36.5°C, after which body temperature was set to 36.5°C for the remainder of the study time. Similar to the effect of recovery of organ function, the effect of TEMP on clearance was tested separately for the renally vs. hepatically cleared drugs.

Correlation in clearance

It was expected that the clearance of the different drugs and metabolites were correlated within one individual, especially between drugs eliminated via the same organ system. Interindividual variability in clearance not explained by BW, GA, PNA, or TEMP was introduced using a log‐normal distribution according to the following equation:

In which CL , is the clearance of compound n in individual i, TVCL is the typical value of clearance of compound n, and η , describes the interindividual variability assuming that all values of η have a normal distribution with mean 0 and SD ω. In the first step, a full OMEGA matrix was considered to study the correlation between interindividual variability of the different clearance components. However, the resulting OMEGA matrix contained 12 diagonal and 66 off‐diagonal parameters to be estimated. Therefore, the model was simplified by introducing one common η value on all clearance parameters describing the correlation between the different components according to the following equations:

In which CL , is the clearance of compound n in individual i, TVCL is the typical value of clearance of compound n, and η , describes the interindividual variability assuming that all values of η have a normal distribution with mean 0 and SD ω. In the first step, a full OMEGA matrix was considered to study the correlation between interindividual variability of the different clearance components. However, the resulting OMEGA matrix contained 12 diagonal and 66 off‐diagonal parameters to be estimated. Therefore, the model was simplified by introducing one common η value on all clearance parameters describing the correlation between the different components according to the following equations:

In which CL ,1 and CL ,2 denote the clearance of compound 1 and 2 for individual i, respectively. TVCL1 and TVCL2 denote the typical values of clearance for compound 1 and 2, respectively, and η ,1 and η ,2 describe interindividual variability of compound 1 and 2, respectively. Correlation is described by η,common, which is the part of interindividual variability that is common for both clearance components, and θ 1 and θ 2 denote the scaling factor for this common variability component. For morphine, this scaling factor was fixed to 1 and all other scaling factors were, therefore, relative to morphine clearance. Interindividual variability in clearance and correlation between different clearance components were subsequently calculated using:

In which ω 1,total and ω 2,total denote the total interindividual variability of compound 1 and 2, respectively; ω 1 2 and ω 2 2 denote the compound‐specific interindividual variability for compound 1 and 2, respectively, and ω common denotes the common interindividual variability for all clearance components. R1,2 describes the correlation in clearance between compounds 1 and 2.

Population pharmacokinetic analysis

A population pharmacokinetic model was developed using the nonlinear mixed effect modeling program NONMEM (version 7.3, Icon Development Solutions) with R (version 3.4.1), Xpose (version 4) for data visualization, and Piraña for run management. Datasets from the individual drugs were restructured and combined so that all data could be fitted simultaneously. Each neonate was included once in the final combined dataset but could contain data from multiple drugs. The previously developed structural models were used with regard to the number of compartments and structural pharmacokinetic parameters estimated. To reduce model complexity, all volumes of distribution and the associated variability were fixed to the parameter estimates from the individual models. As a separate population pharmacokinetic model for lidocaine and MEGX from PharmaCool, data were not available, the volume of distribution for lidocaine and MEGX were fixed to values obtained in a previous population pharmacokinetic study in the same population.

Separate proportional error models for all compounds were used to model residual unexplained variability. For midazolam, OHM, HMG, lidocaine, and MEGX, measurements below the lower limit of quantification (LLOQ) were fixed to LLOQ/2. Therefore, for these substances, an additive error fixed on LLOQ/2 was also included. For the other compounds, no data below LLOQ were present.

Parameter precision was assessed with sampling importance resampling. Both graphical (e.g., goodness‐of‐fit plots) and statistical model evaluation procedures were used to assess model adequacy.

RESULTS

Patient characteristics

In total, 192 neonates from the PharmaCool study were included in this analysis. Patient characteristics are displayed in Table .

Table 1

Patient characteristics

Parameter	Patients (n = 192)
Gestational age; weeks, mean ± SD	39.7 ± 1.66
Birth weight; kg, mean ± SD	3.38 ± 0.617
Male, n (%)	118 (61.5%)

Body size and maturation

No systematic deviation among the conditional weighted residuals for each compound and BW and GA were observed, indicating the appropriateness of these assumptions (Supplementary Materials). The relative influence of BW and GA on drug clearance compared to a BW of 3,500 g and a GA of 40 weeks are depicted in Figure .

Figure 1

The relative influence of birth weight (left) and gestational age (right) on clearance of all drugs in the final pharmacokinetic model.

PNA was identified as a covariate on clearance for all compounds except phenobarbital. After grouping these compounds into renal clearance (amoxicillin, benzylpenicillin, gentamicin, M3G, M6G, HMG, and MEGX), hepatic high clearance (lidocaine) and hepatic intermediate clearance (morphine, midazolam, and OHM), it was found that the effect of PNA on renally cleared compounds was similar to the effect of PNA on lidocaine clearance. Therefore, it was decided to separate the covariate PNA into two groups: high‐clearance compounds (renally cleared compounds and lidocaine) and intermediate‐clearance compounds (morphine, midazolam, and OHM). In the high‐clearance group, the relative effect of PNA on clearance was 1.23%/hours of life (95% confidence interval (CI) 1.03–1.43); in the intermediate‐clearance group, this was 0.54%/hours of life (95% CI 0.371–0.750). More complex forms of the relationship between clearance and PNA than the currently used linear relationship were investigated but proved unidentifiable because only data up to 120 hours after birth were available.

As expected, the effect of PNA by far exceeded the effects of maturation. According to the sigmoidal Hill equation used to describe maturation, clearance would increase with ~ 0.05%/hours of life in the first 5 days after birth in this population. Furthermore, large interindividual variability was found for these effects: 71.6% for the high‐clearance compounds and 54.8% for the intermediate‐clearance compounds. Interindividual variabilities in both effects showed a very high correlation, which was subsequently fixed to 100% to reduce model complexity (Supplementary Material).

The influence of TEMP on clearance was only significant for the intermediate‐clearance drugs (morphine, midazolam, and OHM). During TH, clearance was decreased by 20.5% (6.83%/°C; 95% CI 5.16%/°C–8.34%/°C) compared with normothermia. The influence of both PNA and TEMP on the average clearance for the three identified drug groups are shown in Figure .

Figure 2

The relative influence of postnatal age and body temperature on the clearance for the three identified drug groups in the final pharmacokinetic model. The vertical solid lines indicate therapeutic hypothermia; and the vertical dashed line indicates return to normothermia.

Correlation in clearance was calculated among all compounds except phenobarbital, lidocaine, and MEGX. Phenobarbital was not tested because its clearance visually did not correlate with any other compound. Correlation for lidocaine and MEGX could not be estimated due to the sparseness of the data (only 28 neonates and 77 samples). Correlation in clearance for the remaining compound is presented in Table . All correlations were positive. The highest correlation was between clearance of M3G and M6G (96.2%), as expected. Correlation in clearance in the hepatically cleared compounds was relatively high; the renally cleared compounds did not show a higher correlation within the group than compared with the hepatically cleared compounds.

Table 2

Correlation in clearance

	Morphine	Midazolam	OHM	M3G	M6G	HMG	Amoxicillin	Benzylpenicillin	Gentamicin
Morphine
Midazolam	62.2%
OHM	63.3%	57.0%
M3G	35.4%	31.9%	32.2%
M6G	32.5%	29.3%	29.6%	94.8%
HMG	62.6%	56.4%	57.0%	31.9%	29.3%
Amoxicillin	42.1%	38.0%	38.3%	21.4%	19.7%	37.9%
Benzylpenicillin	58.6%	52.8%	53.4%	31.9%	27.4%	52.8%	35.5%
Gentamicin	45.9%	41.4%	41.8%	23.4%	21.5%	41.4%	27.9%	38.8%

HMG, hydroxymidazolam glucuronide; M3G, morphine‐3‐glucuronide; M6G, morphine‐6‐glucuronide; OHM, 1‐hydroxymidazolam. [Colour version of this table can be viewed at wileyonlinelibrary.com]

Model evaluation demonstrated that the final model was adequate in describing the data. Goodness‐of‐fit plots of observed vs. population and individual predicted concentrations showed no systematic deviation and the weighted residuals were homogeneously scattered vs. predicted values and time for all compounds (Supplementary Materials).

The final structural models from the original publications were used as the structural model for this analysis. , , , , Estimates of the parameters related to clearance from the final model are shown in Table . Full pharmacokinetic parameter estimates from the final model are included in the Supplementary Materials.

Table 3

Pharmacokinetic parameter estimates and SIR results relating to clearance

Compound		Parameter
		Cl, L/h a	PNA on Cl, %/h	TEMP on Cl, %/°C
Morphine	Estimate	0.811	0.540	6.83
	SIR 95% CI b	0.707–0.937	0.371–0.750	5.16–8.34
Midazolam	Estimate	0.511	0.540	6.83
	SIR 95% CI b	0.387–0.620	0.371–0.750	5.16–8.34
OHM c	Estimate	1.72	0.540	6.83
	SIR 95% CI b	1.43–2.05	0.371–0.750	5.16–8.34
M3G c	Estimate	0.241	1.23	NA
	SIR 95% CI b	0.220–0.269	1.03–1.43	NA
M6G c	Estimate	0.765	1.23	NA
	SIR 95% CI b	0.697–0.854	1.03–1.43	NA
HMG c	Estimate	0.111	1.23	NA
	SIR 95% CI b	0.0977–0.126	1.03–1.43	NA
Amoxicillin	Estimate	0.178 d	1.23	NA
	SIR 95% CI b	0.159–0.196	1.03–1.43	NA
Benzylpenicillin	Estimate	0.359 d	1.23	NA
	SIR 95% CI b	0.297–0.423	1.03–1.43	NA
Gentamicin	Estimate	0.108 d	1.23	NA
	SIR 95% CI b	0.0968–0.120	1.03–1.43	NA
Lidocaine	Estimate	0.937	1.23	NA
	SIR 95% CI b	0.783–1.11	1.03–1.43	NA
MEGX c	Estimate	1.51	1.23	NA
	SIR 95% CI b	0.991–2.06	1.03–1.43	NA
Phenobarbital	Estimate	0.00930	NA	NA
	SIR 95% CI b	0.00785–0.0111	NA	NA

BW, body weight; CI, confidence interval; Cl, clearance; GA, gestational age; HMG, hydroxymidazolam glucuronide; M3G, morphine‐3‐glucuronide; M6G, morphine‐6‐glucuronide; MEGX, monoethylglycinexylidide; NA, not applicable; OHM, 1‐hydroxymidazolam; PNA, postnatal age; SIR, sampling importance resampling; TEMP, body temperature.

Estimates for a neonate with BW 3.5 kg, GA 280 days, PNA 0 hours and TEMP 36.5°C.

Six iterations; no. of samples 4,000, 4,000, 4,000, 4,000, 4,000, and 4,000; no. of resamples 250, 250, 250, 250, 250, and 1,000.

All metabolite estimates are relative to their formation fraction.

Estimated clearance in the central compartments. Peripheral compartment estimates are included in the Appendix.

DISCUSSION

This study successfully combined data from 12 compounds in the PharmaCool study into one pharmacokinetic model. Body size related effects were adequately captured with allometric scaling. Previously developed models for maturation were successfully implemented in this population using GA as descriptor. Subsequently, the model was extended to quantify the effects of body temperature and recovery of organ function.

The allometric relationship between BW and clearance is based on physiological observations both within and across species and is the most widely used method to describe size differences. , , , , Although its validity has been debated over the years, especially in neonates, we showed that allometry could describe the relationship with body size in this dataset if also other structural effects, including maturation, were included.

PMA is acknowledged as the most reliable age‐related factor to reflect the biology of clearance maturation. In our study population, PMA was determined by GA for > 98% because data were collected for only a short period after birth. Therefore, we included GA as a descriptor for maturation, which also enabled us to separate the effect of maturation from the effect of organ recovery. In the literature, maturation of both renal and hepatic clearance has been described using PMA with a sigmoid maximum effect (Emax) model with similar values for Hill coefficient and TM50. , We hypothesized that this function would adequately describe maturational differences in clearance for all compounds. Careful model evaluation supported this hypothesis. It should, however, be noted that full evaluation of maturation in this population is hampered by the small range in GA and a relatively short observation period.

After including BW and GA as fixed effects, clearance proved to be strongly dependent on PNA. This effect was much larger than can be expected from maturation alone. Therefore, we hypothesized that this was caused by recovery of organ function after asphyxia. The effect of PNA was largest in the high‐clearance compounds. For these drugs and metabolites, clearance is closely linked to organ perfusion. As damaged liver and kidney cells regenerate, organ perfusion and thereby clearance will increase rapidly. Clearance of intermediate‐clearance compounds is not only dependent on organ perfusion but also on enzyme capacity. Enzymes, such as CYP3A and UGT2B7, might show slower recovery or may even be unaffected by asphyxia, which could explain the smaller effect of PNA in this group compared with the high‐clearance compounds. No effect of PNA was identified on phenobarbital clearance. As phenobarbital is an extremely low clearance drug for which clearance is independent of hepatic perfusion, this finding was not unexpected. Interindividual variability on the effect of PNA on clearance was 71.6% and 54.8% for high and intermediate‐clearance compounds, respectively. As the severity of asphyxia differs between neonates with NE, so will recovery of organ function. A much stronger recovery in less severely ill neonates vs. no recovery at all in the most severely ill is common in this population and might explain the relatively high variability, which supports the interpretation that the large effect of PNA on clearance is most likely caused by recovery of organ function. As recovery of organ function is closely linked with NE, it is unlikely that such a strong dependency of clearance on PNA is also present in neonates not suffering from NE. It should be noted that in the final model a linear relationship between PNA and clearance was assumed based on data up to 5 days after birth. Extension of this empirical relationship beyond the range of these data should be done with great caution, as recovery of organ function will likely reach a plateau later in life.

TEMP could only be identified as covariate on clearance of the intermediate‐clearance compounds. As depicted in Figure , clearance in this group drops after birth due to TH, after birth followed by a hampered increase attributed to PNA during the hypothermic phase, and a steeper increase during rewarming. Although this is somewhat contradictory to the individual models, we believe that in the current model the relatively small decrease in clearance caused by a lower temperature is not identifiable due to the much more pronounced effect of PNA in the high‐clearance group. In the intermediate‐clearance group, the effect of TEMP on clearance could be identified separately from the smaller effect of PNA.

This study is the first to integrate data from some of the most important and most frequently used drug in neonates treated with TH for NE into one integrated population pharmacokinetic model. Previously, the population pharmacokinetic model for amoxicillin was successfully applied to predict clearance of benzylpenicillin in the PharmaCool study population. The present model included compounds with hepatic and renal elimination and has identified covariates on clearance that transcend individual drugs and routes of elimination. This integrated model can be used to predict clearance of other drugs in this population based on data from older children or adults. Levetiracetam is an AED that is used increasingly in neonates treated with TH for NE. Based on its clinical pharmacokinetic profile, levetiracetam can be grouped with the high‐clearance drugs and it can be expected that clearance increases strongly with PNA in this population. , The same prediction can be made for 2‐iminobiotin and allopurinol, high‐clearance drugs that are currently being investigated for additional neuroprotection in combination with TH. , , ,

Drug dosing is highly challenging in neonates in general, and may be even more difficult in critically ill encephalopathic neonates treated with TH. Studies in this populations are difficult to perform and, therefore, it is of importance to elucidate and quantify the processes that influence the pharmacokinetics of drugs used in this population. Although our model cannot be used to fully describe the pharmacokinetics of any (new) drug administered to neonates treated with TH for NE, it can be used to predict changes in drug clearance. Whereas individual pharmacokinetic studies will still be necessary, knowledge obtained from this integrated model can facilitate the design and might reduce the number of patients needed for those studies.

CONCLUSION

Data from seven different drugs (12 compounds) administered to neonates treated with TH for NE were successfully combined into one integrated population pharmacokinetic model. PNA was identified as a covariate on clearance of both high‐clearance and intermediate‐clearance compounds and TEMP was subsequently identified as covariate on intermediate‐clearance compounds. Individual clearance values were positively correlated for nine compounds. This integrated model can be used to facilitate drug dosing and future pharmacokinetic studies for other drugs in this population by predicting (changes in) drug clearance based on the clinical pharmacokinetic properties of that drug.

Author Contributions

L.F. and A.H. wrote the manuscript. L.F., T.d.H., Y.B., C.R., D.N., and F.G. performed the research. L.F., Y.B., R.M., F.G., and A.H. analyzed data. T.d.H., C.R., T.E., and D.N. reviewed the manuscript. T.d.H., C.R., T.E., and F.G. designed the research. Y.B., R.M., and F.G. reviewed the manuscript.

Conflicts of Interest

The authors declared no competing interests for this work.

Funding

The PharmaCool study was funded with a grant from the Netherlands Organization for Health Research and Development (ZonMw) Priority Medicines for Children (Grant number: 40‐41500‐98‐9002).

Fig S1

Click here for additional data file.

Fig S2

Click here for additional data file.

Fig S3

Click here for additional data file.

Table S1

Click here for additional data file.