The data used for analysis have been previously published (13). Briefly, the data consisted of sCr observations in (near)term neonates with moderate-to-severe HIE, who underwent TH during the first 3 days after birth. The data covered the first 10 days of life. The available covariates were birth weight (BWT), gestational age (GA), survival (DEATH, neonatal death, day 1–28), and acute kidney injury (AKI). AKI detection was based on the KDIGO definition (any AKI): sCr↑≥0.3 mg/dL within 48 h or sCr↑≥1.5 fold versus the lowest prior sCr within 7 days, irrespective of urine output, as these data were unavailable (13). Information on the severity (moderate or severe) is not available in the pooled sCr dataset. The reference population comprised of neonates with a GA of at least 24 weeks with at least one measurement of sCr during 42 days after birth and sCr value ≤ 2 mg/dL who were included in the model development dataset of (12). The reference dataset consisted of 1080 subjects with 7977 sCr observations. The descriptive statistics of covariates for TH and reference populations are presented in Table I. The sCr values obtained by the Jaffe assay were converted to values equivalent to ones obtained by an isotope dilution mass spectrometry (IDMS) traceable enzymatic assay using the following Eq. (1):

$$}_}=1.003\bullet }_}+0.057$$

(1)

The sCr values of obtained by an unidentified assay were excluded from the analysis. To facilitate application of a population model that was developed for the reference population, the data were further reduced to include subjects with sCr value ≤ 2 mg/dL and GA ≤ 42 weeks not exceeding the reference ranges. The final dataset consisted of 975 subjects with 4636 sCr observations (initial dataset had 1136 subjects and 4724 sCr observations). A spaghetti plot of individual sCr time courses is shown in Fig. 1. SCr data were reported in mg/dL (conversion, 1 mg/dL = 88.4 µmol/L). A histogram of the GA distribution is shown in Fig. 1S.

Time courses of individual serum creatinine (sCr) measurements for 975 in (near)term neonates with moderate-to-severe hypoxic-ischemic encephalopathy who underwent therapeutic hypothermia (TH) during the first 3 days of life. The bold (red) line is the mean of observed values

We adopted our previously introduced population model of sCr kinetics in neonates (see Fig. 2) (12):

$$s\mathrm=\left(1+\frac)}}\right)\frac_}}(\mathrm)}$$

(2)

where \(_}\) is the rate sCr is secreted from the muscle, \(Q(\mathrm)\) is a back-flow of creatinine from the renal tubules to the plasma, and \(\mathrm(\mathrm)\) is the glomerular filtration rate at the postnatal age \(PNA\). Q(PNA) was modeled by the truncated Gaussian distribution function designed to describe peak and a gradual disappearance of the back-flow:

$$Q\left(\mathrm\right)=\frac}_}}}\bullet \mathrm\left(-\frac}^}\bullet -}_}\right)}^\right)$$

(3)

where \(_}\) is the peak of the back-flow that occurs at time PNAp. The parameter \(K\) is the time scale factor that controls the width of Q (1/K is equal to the standard deviation of the Gaussian distribution). \(\mathrm\left(\mathrm\right)\) was described by the sigmoidal Emax model:

$$\mathrm\left(\mathrm\right)=}_+\frac}_}-}_)\bullet }^}}_^+}^}$$

(4)

where \(}_\) is the GFR at birth, \(}_\) is the time after birth at which 50% of GFRss is reached, and \(\gamma\) is the Hill coefficient.

Schematic diagram of pharmacokinetic model of creatinine. Serum creatinine (sCr) is synthesized from the muscle at a zero-order rate (ksyn) and cleared by the kidney glomerular filtration (GFR(t)) into the renal tubules (rtCr). The creatinine from the renal tubules is excreted into the urine by the urinary flow (UF). Transiently, rtCr is also secreted into the plasma by the back-flow (Q(t))

All identifiable model parameters (P) were assumed to be log-normally distributed among subjects:

$$P=_}\mathrm\left(_}\right)\;\mathrm\;_}\sim N(0,_}^)$$

(5)

where \(_}/\mathrm\),\(}_}, K, }_}/}_, _}/}_, }_, \gamma \}\). The \(\mathrm\) was identified as the only significant covariate (alternative to the body weight) and influencing the model parameters \(_}/\mathrm\),\(}_}, K, }_}/}_, _}/}_, }_,\}\) through the linear relationship

$$P=_}\left(1+_\_\mathrm}\left(\mathrm-}_}\right)\right)\mathrm\left(_}\right)$$

(6)

where \(}_}=34.2\) weeks.

The additional covariates in addition to \(\mathrm\) for the current analysis were acute kidney injury (\(\mathrm)\) and patient’s death, defined as neonatal death (day 1–28) (\(\mathrm\)). Both covariates were dichotomous with the value 1 indicating presence and 0 absence of the covariate:

$$_}=_,\mathrm}\bullet \mathrm+_,\mathrm}\bullet (1-\mathrm)$$

(7)

The constant residual error model was applied to the log-transformed sCr data:

$$\mathrm\left(}_}\right)=\mathrm\left((\mathrm}_})\right)+_} \;\mathrm\;_}\sim N(0,^)$$

(8)

where \(}_}\) are sCr observations for subject \(i\) at age \(}_}\).

The performance of the population model was evaluated using the diagnostic plots observed vs. predicted values, conditional individual weighted residuals (CWRES), conditional individual weighted residuals (CIWRES), and visual predictive check (VPC). To obtain VPC, 500 simulations of the data were performed with parameter estimates from the final model. Simulated 5th, median, and 95th percentiles and their 95% CIs were compared with observed values. The plots were generated in RStudio.

Individual time courses of \(\mathrm\) and \(\mathrm/}_\) for TH and a reference neonate (\(\mathrm=39\) weeks, \(\mathrm=0\), and \(\mathrm=0\)) were simulated using estimated model parameters for both populations. For each population, 700 individual time courses with no residual error were generated at PNA increments 0.1 day over the PNA range of 10 days. At each PNA, the median, 5th, and 95th percentiles of \(\mathrm\) and \(\mathrm/}_\) were calculated and plotted against PNA. Additionally, for TH neonates, analogous plots were generated of \(\mathrm=39\) weeks, \(\mathrm=1\), and \(\mathrm=0\).

All models were implemented in NONMEM 7.4 (ICON Development Solutions, Ellicott City, MD). Population model parameters were estimated using the importance sampling (IMP) method with CTYPE = 1 convergence criterion. The model performance diagnostic plots were obtained using R packages: ggplot2, lattice, and vpc (14) implemented in RStudio version 2022.07.2 (RStudio Inc., Boston, MA).