In total, 58 publications were identified through the PubMed search. After contacting the corresponding or senior authors by e-mail, we obtained individual-level data from 13 identified publications [7,8,9, 11,12,13,14,15, 32, 45,46,47,48] to generate a pooled dataset which covers a broad range of patients ranging from premature neonates [8] to very elderly people [11, 45, 46, 48], from underweight [11, 13, 46,47,48] to obese adults [7, 11, 13,14,15, 32, 45, 47, 48] and from critically ill patients with sepsis [11, 13, 15, 32, 46, 47] to febrile neutropenic patients with haematological malignancy [7]. The included PK data and patient characteristics are summarised in Table 1 and their distributions are shown in Figure S1 of the ESM1. In total, 3798 PIP concentrations and 1948 TAZ concentrations in different types of samples derived from 415 patients were included in this population PK analysis. The PIP observations comprised 2855 total plasma concentration (CTOT) observations, 888 CUNB observations and 55 observations of concentration in dried blood spot samples (CDBS). However, for TAZ there were only CTOT observations (n = 1893) and CDBS observations (n = 55). The included patients consisted of 32 newborns, 127 children/adolescents, 74 elderly patients, 19 very elderly patients and 163 adults including 42 obese and 12 underweight. For 191 individuals, both PIP and TAZ concentrations were available. For 20 patients both CTOT and CUNB observations for PIP were available. For 28 patients both CTOT and CDBS observations for PIP and TAZ were available.

The developing hierarchy of single-drug models for PIP and TAZ is shown in Table 2. Generally, these two single-drug models were established through the same procedures. Both of them are 2-compartment models with linear elimination and identified covariates including TBW, PMA and SCR. Besides, PK alterations were observed for PIP in two component datasets. In the Sime et al. [7] study, which included patients with haematological malignancies, we found an elevated CLPIP (+ 72.1%) and a lower V2_PIP (− 48.3%). In the Sukarnjanaset et al. [11] study, a high fUNB_PIP (100%) was observed. In the single-drug model for PIP, an additional set of proportional and additive residual errors is used for CUNB observations, which is independent of that for CTOT and CDBS observations. For both drugs, non-linear elimination was tested but model fit was not significantly improved. We also tested non-linear protein binding of PIP using an Emax function. The estimated dissociation constant was high (854 mg L−1) and the fit of the model to the data was not significantly improved (ΔAIC = − 0.12). Therefore, non-linear protein binding was not included in the single-drug model for PIP. More details about the model development process of single-drug models are described in the ESM2.

The combined population PK model was established based on the final single-drug models for PIP and TAZ. We considered IIV terms, covariate fixed effects and residual error terms for combination across models. No apparent deterioration in model performance was caused by combining the θSCR of PIP and that of TAZ into one (ΔAIC = − 1.995), indicating that SCR influences CLPIP in the same way as CLTAZ. For every 0.20 mg dL−1 rise in SCR, CLPIP and CLTAZ both decrease by 6.7% according to the final combined model.

The maturation functions for the two drugs could be merged without significant changes to model fit (ΔAIC = − 0.38) with CLPIP and CLTAZ reaching 50% maturation at 54.2 weeks PMA (MAT50). The decline functions could not be merged (ΔAIC = + 13.7) and in the final model, CLTAZ declines by 50% at 61.6 years PMA (DEC50_TAZ) whereas for CLPIP this is 89.1 years PMA (DEC50_PIP). The typical-for-PMA standardised CLPIP (L h−1 70 kg−1) for all included patients is shown in Fig. 1A (solid line). For comparison, the maturation-decline function for PIP (dashed line), extracted from the pooled analysis by Lonsdale et al. [6] are also shown. As for TAZ, the typical-for-PMA standardised CLTAZ (L h−1 70 kg−1) for all the included subjects with TAZ observations is shown in Fig. 1B.

Standardised clearance [CLstd] (L h−1 70 kg−1) of piperacillin (PIP, A) and tazobactam (TAZ, B) throughout life. The solid lines represent the typical CLstd of PIP and TAZ according to our final combined pharmacokinetic (PK) model. The solid grey circles represent the post hoc CLstd values for all patients included in this study. The region between the 5% and 95% percentile of those post hoc CLstd values is shown with grey shadow. The maturation-decline function for PIP according to the Lonsdale et al. [6] study is shown by the dashed line

Merging the ratio of CDBS to CTOT for PIP (fDBS_PIP) and the one for TAZ (fDBS_TAZ) together resulted in a worse fit (ΔAIC = + 8.82), revealing that significant differences existed between the estimates for this pair of parallel parameters. As estimated in the final combined PK model, concentrations of PIP and TAZ in dried blood spot samples are 63.2% and 55.2% lower, respectively, than those in plasma.

As for the other fixed-effect parameters of covariate-based and study-specific corrections, their estimates were slightly influenced although they were not enrolled into the correlation test. The fUNB_PIP was estimated as 64.5% in the final combined PK model. In the patients derived from the study of Sukarnjanaset et al. [11], a higher fUNB_PIP was observed (fUNB_Sukarnjanaset = 100%). In addition, an increase by 73.2% and a decrease by 48.8% were sequentially identified for CL and V2 in the Sime et al study, which included patients with haematological malignancies [7].

Later, we tested the correlation within every pair of parallel IIV terms. Ultimately, we could combine the IIV for PIP and TAZ for V1, V2 and Q2, respectively, which produced the lowest AIC (ΔAIC = − 370.844) without interfering with the interpretation of the fixed-effect parameters. Combining the IIV terms of CLPIP and CLTAZ reduced the AIC (ΔAIC = − 6.57) but caused an increase in the estimate for DEC50_TAZ (+ 19.1 years). When the approach of covariance estimation was performed (instead of combining eta terms as in the final model) fewer pairs of IIV terms could be implemented (V1 and V2, but not CL or Q2) before numerical issues occurred. The resulting model fitted the data more poorly than the final model (ΔAIC = + 95.2). Similarly, when estimating different eta magnitudes, fewer pairs of IIV terms could be implemented (V1 and V2, but not CL or Q2) before numerical issues occurred and the model fitted the data more poorly than the final model (ΔAIC = 91.2).

Merging each pair of parallel residual error magnitude resulted in a significant increase of AIC (ΔAIC = + 7.39 to + 9.133), indicating that the residual-error distributions for PIP are significantly different from those for TAZ. We also considered including off-diagonal elements in $SIGMA but the complexity of the error model with proportional-additive (PIP total/DBS, and unbound) and proportional (TAZ) components made it difficult to construct a model capturing this intent.

After completing the evaluation of correlations between the PK of PIP and TAZ, we tested a model consisting of a parallel linear and non-linear elimination pathway for PIP (for details the reader is referred to ESM2 and in particular Eqs S1–S3). The inclusion of this parallel linear and non-linear elimination pathway did not improve the fit of the model to the data (ΔAIC = + 2.3).

The final combined PK model for PIP/TAZ was thereby obtained, which is shown in Table 3 and Eqs. 6–25.

$$\begin_} \left(\text\right)=_}\times _}\times ^_},\end$$

(6)

$$\begin}_} \left(}^\right)=_\_\text}\times _}}^\times _}\times _\_\text}\times _}\times _\_\text}\times e}^_}\end,$$

(7)

$$\begin_} \left(\text\right)=_2\_\text}\times _}\times _2\_\text}\times ^_},\end$$

(8)

$$\begin_} \left(}^\right)=_2\_\text}\times _}}_}}\right)}^\times ^_},\end$$

(9)

$$\begin}_} \left(}^\right)=\frac_}}_}}\times _\_\text}\times _\_\text},\end$$

(10)

$$\begin_} \left(\text\right)=_}\times _}\times ^_},\end$$

(11)

$$\begin}_} \left(}^\right)=_\_\text}\times _}}^\times _}\times _\_\text}\times _}\times e}^_}\end,$$

(12)

$$\begin_} \left(\text\right)=_}\times _}\times ^_},\end$$

(13)

$$\begin_} \left(}^\right)=_}\times _}}_}}\right)}^\times ^_},\end$$

(14)

$$\begin}_} \left(}^\right)=\frac_}}_}}\times _\_\text},\end$$

(15)

$$\begin_}=\frac\left(\text\right)},\end$$

(16)

$$\begin_}=\frac\left(\text\right)}^_}}\left(\text\right)}^_}+}_}^_}},\end$$

(17)

$$\begin_\_\text}=1-\frac\left(\text\right)}^_}}\left(\text\right)}^_}+}_\_\text}}^_}},\end$$

(18)

$$\begin_\_\text}=1-\frac\left(\text\right)}^_}}\left(\text\right)}^_}+}_\_\text}}^_}},\end$$

(19)

$$\begin}_}=}^\left(\text\left(\text\right)/100\right)\right)}\left(\text\right)/100}}\right]},\end$$

(20)

$$\begin_\_\text}=\left\1.73,\quad for\, the\, study\, by \,Sime \,et\, al. \,\left[7\right]\\ 1,\quad for\, other\, included \,studies\,\end\right.,\end$$

(21)

$$\begin_2\_\text}=\left\0.512,\quad for\, the\, study\, by\, Sime\, et \,al. \,\left[7\right]\\ 1,\quad for\, other\, included\, studies\end\right.,\end$$

(22)

$$\begin_\_\text}=\left\0.645,\quad for \,_}\, observations\, of\, PIP \,which \,are\, not\, from\\ the \,Sukarnjanaset \,et \,al.\, \left[11\right]\, study \\ \\ 1, \quad for\, other \,observations\qquad\qquad\qquad\qquad\qquad\end,\right.\end$$

(23)

$$\begin_\_\text}=\left\0.368,\quad for\, _}\, observations\, of \,PIP\\ 1,\quad for\, other\, observations\quad\ \end\right.,\end$$

(24)

$$\begin_\_\text}=\left\0.448,\quad for\, _}\, observations \,of \,TAZ\\ 1,\quad for\, other\, observations\quad\ \ \end\right.,\end$$

(25)

In above equations, the parameters with subscript PIP are only used to describe the PK of piperacillin while those with subscript TAZ are exclusively used to describe the PK of tazobactam. V1 and V2 are the central and peripheral volume of distribution; CL and Q2 denote the elimination and inter-compartment clearance. Size-related changes, PMA-induced maturation, PMA-induced decline, and SCR-related changes in the PK of PIP/TAZ are described by FSIZE, FMAT, FDEC and FSCR, respectively. The θCL_Sime and θV2_Sime represent the elevated CLPIP and decreased V2_PIP in patients with haematological malignancies [7]. ηi (i = 1–5), with variances of ωi, represent IIV of typical PK parameters. The fUNB and the fDBS are the fraction unbound and the ratio of CDBS to CTOT. IPRED represent individual predictions of all kinds of observations. A1 denotes predicted amounts in the central compartments.

Backwards elimination of FSIZE, FMAT, FDEC or FSCR led to significant OFV increases and difficulties in convergences. Goodness-of-fit and pvcVPC plots for the final combined PK model for PIP/TAZ are shown in Figs. 2, 3, 4. In addition, Figures S3–S4 in the ESM1 show the GOF and pvcVPC plots stratified by observation type and patient subgroup. Our combined model shows acceptable performance across these diagnostics, despite the apparent underprediction of PIP in the first 1–2 hours after stopping the infusion (Fig. 4). The underprediction does not seem to be specific to any subgroup (Fig. S4). We were unable to remove this underprediction in model development. Model code of the final PIP/TAZ population PK model is available in ESM3.

Goodness-of-fit plots for the final combined population pharmacokinetic model for piperacillin concentrations. Scatterplots show the distributions of observed piperacillin concentrations versus population and individual predictions, conditionally weighted residuals (CWRES) versus population predictions and time after the end of the dose, as well as normalised prediction distribution errors (NPDE) versus population predictions and time after the end of the dose. Circles, triangles and crosses denote total plasma concentrations, unbound plasma concentrations, and concentrations in dried blood spots, respectively. Solid black lines represent lines of unity or zero lines. Red dashed lines are non-parametric smoothers of those distributions. Negative time points mean that observations were collected during the infusion, while positive time points denote observations taken after the infusion is finished

Goodness-of-fit plots for the final combined population pharmacokinetic model for tazobactam concentrations. Scatterplots show the distributions of observed tazobactam concentrations versus population and individual predictions, conditionally weighted residuals (CWRES) versus population predictions and time after the end of the dose, as well as normalised prediction distribution errors (NPDE) versus population predictions and time after the end of the dose. Circles and crosses denote total plasma concentrations and concentrations in dried blood spots, respectively. Solid black lines represent lines of unity or zero lines. Red dashed lines are non-parametric smoothers of those distributions. Concentrations with negative times were collected during the infusion, while positive times denote observations taken after the infusion was finished

Prediction- and variability-corrected visual predictive check plots for the final combined population pharmacokinetic model for total plasma concentrations of piperacillin (A), unbound plasma concentrations of piperacillin (B), piperacillin concentrations in dried blood spots (C), total plasma concentrations of tazobactam (D) and tazobactam concentrations in dried blood spots (E), respectively. Solid red lines are the 50th percentiles of observations corrected by prediction and variance, while dashed red lines denote the 10th and 90th percentiles. Red dashed lines are non-parametric smoothers of those distributions. Negative time points mean that observations were collected during the infusion, while positive time points denote observations taken after the infusion is finished. Grey shaded rectangles represent the 95% confidence intervals for the simulated 10th, 50th and 90th percentiles of the prediction- and variance-corrected observations

Simulated dosing regimens and characteristics of the virtual patients were summarised in Table S3 of the ESM4. The SmPC label does not provide specific dosing recommendations for neonates and infants. In the simulations for the SmPC dosing recommendations, the lowest weight-based dose from the SmPC label (70/8.75 mg kg−1 PIP/TAZ every 8 h) was applied for both age groups. Similarly, the lowest weight-based dose from the FDA label (80/10 mg kg−1 PIP/TAZ every 8 h) was used for neonates in the simulations for the FDA dosing recommendations.

The simulated steady-state PTA versus MIC profiles resulting from the dosing recommendations in the SmPC label are shown in Fig. 5.

Simulated probability of target attainment (PTA) of piperacillin (PIP) at steady state for different age groups of virtual patients with the highest dose according to the summary of product characteristics (SmPC) label. The PTAs of the target that unbound plasma concentration of PIP remain above the minimum inhibitory concentration (MIC) for the whole dosing interval (PTA fT>MIC 100%) and the PTAs of the target that the unbound plasma concentration of PIP remain above four times the MIC for the whole dosing interval (PTA fT>4*MIC 100%) are shown in the top and the bottom rows, respectively. PTAs under intermittent, extended, and continuous infusion are shown in the left, the middle and the right columns, respectively. PTA for different age groups are shown in different colors. The shadow areas represent the 95% confidence intervals of the PTA versus MIC curves. The dashed lines denote the PTA of 90%

Simulations for both the SmPC (Fig. 5) and the FDA dosing recommendations (Fig. S5) show that PTA versus MIC profiles are considerably different across age groups, suggesting that current dosing recommendations do not result in consistent PTAs across age groups. The highest PTAs are found in neonates and the lowest in infants, even though dosing recommendations for these groups are weight adjusted. A similar but smaller difference is apparent between elderly and adults, as well as adolescents and children, with PTAs being higher in elderly compared to adults and higher in adolescents compared to children.

The PTA versus MIC profiles of the SmPC dosing recommendations (Fig. 5) and the FDA dosing recommendations (Fig. S5) shift to the right for longer duration infusions indicating higher PTAs are obtained. The PTAs are highest and have lowest variability for continuous infusions, indicating that most patients receive an effective treatment.

In addition, Cm_TAZ was calculated across different age groups to evaluate the steady-state TAZ exposure achieved by the recommended dosing regimens. The median unbound TAZ concentration (considering fUNB_TAZ = 70% [49]) across age groups and simulated scenarios was 3.69 mg L−1 (ranging from 1.64 to 6.77 mg L−1), with the lowest exposure occurring in the infants.