A single-center prospective study was conducted in the intensive care unit (ICU) of the First Affiliated Hospital of Army Medical University from August 2021 to July 2024. Inclusion criteria: (1) patients aged ≥18 years; (2) patients receiving intravenous PMB for CRO infections confirmed by pathogen testing; (3) patients receiving at least four consecutive doses of intravenous PMB, or a loading dose followed by at least three consecutive doses. Exclusion criteria: (1) patients with missing clinical data; (2) patients hospitalized for fewer than 7 days; (3) patients receiving any form of renal replacement therapy during PMB treatment; (4) pregnant women. The research protocol was approved by the Ethics Committee of our hospital (No. (A) KY2021064).

Data collected from electronic medical records included: (1) demographic characteristics and main diseases; (2) PMB dose, administration route, and concurrent antibacterial agents; (3) routine blood test results; (4) liver and renal function indices, with CrCL calculated using the Cockcroft-Gault equation; (5) blood coagulation parameters; (6) other treatments such as extracorporeal membrane oxygenation (ECMO), mechanical ventilation, as well as relevant parameters such as the duration of these treatments.

After at least 48 h of treatment, two blood samples were collected: one immediately before the infusion and the other immediately after. The exact times of blood sampling and infusion were recorded for each patient. Plasma was separated by low-temperature, low-speed centrifugation and stored at −70°C until analysis.

PMB concentrations were analyzed using a validated ultra-performance liquid chromatography-tandem mass spectrometry (UPLC-MS/MS) method. Briefly, using polymyxin E2 as the internal standard, the assay was linear over 0.2–20.0 mg·L^-1 and 0.05–5 mg·L^-1 (r > 0.995) for PMB1 and PMB2 respectively. Intra-day and inter-day precision tests showed relative standard deviations (RSDs) ≤ 12.06%. The average extraction recovery ranged from 103.04% to 117.44%, and RSDs for the matrix effect and stability tests did not exceed 7.42%. PMB concentration was calculated as follows: total concentration of PMB = [PMB1 concentration/PMB1 molecular + PMB2 concentration/PMB2 molecular]*total PMB molecular (Liu et al., 2023).

The PopPK model was developed using nonlinear mixed-effects modeling software: NONMEM (version 7.5.1, ICON plc, United States), Pirana (version 23.1.2, Certara L.P., United States), and PsN (version 5.3.0, https://uupharmacometrics.github.io/PsN/), with the first-order conditional estimation method including interaction (FOCE-I). The base model was selected based on goodness-of-fit (GOF) diagnostic plots, relative standard error (RSE) of parameters, and the objective function value (OFV). Spearman correlation was used to evaluate the relationships between covariates and individual empirical Bayesian estimates (EBEs) of PK parameters before covariate selection. Covariates included age, weight, APACHE II score, the presence of sepsis, and all laboratory parameters listed in Table 1. A decrease in OFV >3.84 (P < 0.05, χ², df = 1) for forward addition and an increase >6.63 (P < 0.01, χ², df = 1) for backward elimination were the criteria for covariate inclusion. GOF plots were used to assess the model. To assess the stability of the final model and the precision of the PK parameters, a bootstrap method with 1,000 resampling iterations was performed. The final model was then subjected to 1,000 simulations for visual predictive checks (VPC) to evaluate the model’s predictive ability and accuracy. The accuracy and predictive performance of the model were further assessed using normalized prediction distribution errors (NPDE) plots. Statistical validation of the model was conducted using the t-test, Fisher’s test, Shapiro-Wilk test, and the Global test.

The patients were divided into a modeling group and a validation group based on the chronological order of enrollment and in accordance with the inclusion and exclusion criteria. External validation of the final model was performed using the validation group data. The mean prediction error (MPE), mean absolute prediction error (MAPE) (Equations 1, 2), F₂₀ (20% of the absolute value of the prediction error), and F₃₀ were calculated by comparing the predicted values with the observed values. Model performance was considered acceptable if MPE% ≤ ±20%, MAPE% ≤ 30%, F₂₀ ≥ 35%, and F₃₀ ≥ 50% (Mao et al., 2018). MPE % = 1 N ∑ i = 1 N pred i − obs i obs i × 100 (1) MAPE % = 1 N ∑ i = 1 N pred i − obs i obs i × 100 (2) pred_i denotes the i-th predicted value, and obs_i its corresponding observed value.

Using the final PopPK model, 1,000 simulations were conducted for commonly used maintenance doses in ICU patients, specifically 50 mg, 75 mg, 100 mg, and 125 mg every 12 h (q12h) with a 1-hour infusion time. According to clinical research (Tang et al., 2023), PMB has demonstrated better clinical efficacy in treating CRO caused nosocomial pneumonia when AUC_ss,24h/minimum inhibitory concentration (MIC) ≥ 66.9. Given that the majority of patients in this study (65/80, 81.25%) had pulmonary infections, we selected AUC_ss,24h/MIC ≥ 66.9 as the PK/PD target, with a 90% target attainment probability (PTA) was considered to be effective. Furthermore, to mitigate the potential risk of nephrotoxicity associated with excessively high AUC_ss,24h, we defined the therapeutic window for AUC_ss,24h as 50–100 mg·h·L^-1 (Tsuji et al., 2019; Yang et al., 2022). We analyzed the probability of the simulated population achieving this therapeutic window.

Data analysis was performed SPSS (version 26.0, IBM, United States) and GraphPad Prism (version 8.3.0, CA, United States). Variables with a normal distribution were expressed as mean ± standard deviation (SD), while non-normally distributed variables were reported as median (interquartile range, IQR). Categorical data were presented as percentages (%). Comparisons between two groups were conducted using the Mann-Whitney U test for non-normally distributed data and the independent t-test for normally distributed data. The estimation of AUC_ss,24h used the two-point method was presented in Equations 3–5 (Liu et al., 2023), while the estimation of AUC_ss,24h based on the PopPK model was obtained by dividing the total 24-hour drug dose by the individual clearance. Bland-Altman analysis was performed to evaluate the consistency of AUC_ss,24h estimates obtained by the two methods. AUC_ss,24h values were categorized as “within the therapeutic window” (50–100 mg·h·L^-1), and McNemar’s test and Kappa test were applied to assess the consistency of AUC_ss,24h classifications. k e = ln ⁡ C peak , ss − ln ⁡ C trough , ss τ − Infusion (3) C soi ′     = C peak , ss e − k e ⋅ Infusion (4) AUC ss , 24 h = C soi ′ − C trough , ss k e × n (5) Infusion time; τ: dosing interval; Csoi' is the exploratory concentration at the start of dosing based on the one-compartment linear elimination pharmacokinetic assumption; k_e: elimination rate constant; n is the number of doses within 24 h.

The modeling group consisted of 80 patients with 184 PMB blood concentration samples, while the validation group included 15 patients with 30 samples. 11 patients in the modeling group underwent repeated sampling. Clinical characteristics and laboratory parameters are summarized in Table 1, showing no significant differences between the two groups. All patients were infected with CRO, in the modeling group, carbapenem-resistant Acinetobacter baumannii (CRAB) was the most common infection (70 cases), followed by carbapenem-resistant Enterobacterales (CRE; 56 cases, comprising 45 Klebsiella pneumoniae, 4 Enterobacter cloacae, 3 Serratia marcescens, 2 Escherichia coli, and 2 Citrobacter freundii), and carbapenem-resistant Pseudomonas aeruginosa (CRPA; 15 cases). In the validation group, CRE was more prevalent (11 cases) than CRAB (10 cases). In the modeling group, the MICs of PMB against CRO strains were: ≤0.5 mg·L^-1 in 72.25% (57/80), 1 mg·L^-1 in 13.75% (11/80), 2 mg L^-1 in 10% (8/80), and 16 mg·L^-1 in 3.75% (3/80). In the validation group, MIC values were ≤0.5 mg·L^-1 in 60% (9/15), 1 mg·L^-1 in 26.67% (4/15), and 2 mg·L^-1 in 13.33% (2/15), with no statistically significant difference observed (P = 0.21).

A one-compartment model with first-order elimination best fit the population data of PMB in critically ill patients. Inter-individual variability was described using an exponential random effects model, while residual variability was described using both proportional and additive error models. Among the covariates evaluated, CrCL and platelet count (PLT) count were found to significantly influence CL, whereas no covariates had a significant effect on V_d (volume of distribution). The final PK model Equations 6, 7 is as follows: CL = 2.03 × CrCL / CrCL median 0.26 × PLT / PLT median − 0.14 × exp η CL (6) V = 18 L × exp η v (7)

Parameter estimates and GOF plots for the base model were presented in Supplementary Table S1 and Supplementary Figure S1 (Supplementary Material). The GOF plot (Figure 1) of the final model demonstrated strong agreement between observed and predicted values. The bootstrap median (Table 2) closely matched the population estimates of the final model, further confirming the robustness of the population PK model. The observed-versus-predicted plot showed a random scatter, indicating no systematic bias and suggesting that the model accurately describes the concentration data. The VPC plot (Figure 2) demonstrated that the median of the model predictions closely aligns with the median of the observed data, and most observed data points fall within the model’s 95% prediction interval. This indicated that the model has strong predictive ability and reasonable variability. Statistical tests for NPDE results included: t-test (P = 1), Fisher variance test (P = 1), Shapiro-Wilk normality test (P = 0.30), and Global test (P = 0.30). The NPDE histogram and Q-Q plot (Figure 3) showed that prediction errors were close to zero and symmetrically distributed, conforming to the normality assumption. NPDE showed no significant variation over time or across predicted concentrations, indicating that the model’s predictive performance is consistent and stable across various time points and concentration levels. External validation, conducted with 15 patients using the final model, showed an MPE% of 2.69%, MAPE% of 28.45%, F₂₀ of 36.67%, and F₃₀ of 73.33%, confirming acceptable predictive performance.

A total of 106 AUC_ss,24h values were obtained from the modeling group (91 values) and the validation group (15 values). The AUC_ss,24h estimated by the PPK model within the range of 50–100 was 51.89% (55/106), while the AUC_ss,24h estimated using the two-point method was 46.23% (49/106). The Bland-Altman analysis results were shown in Figure 4. The mean difference of AUC_ss,24h between the two methods was 12.98 mg·h·L^-1, indicated a small average difference between the methods. The majority of the sample differences fell within the 95% confidence interval, demonstrated consistency between the two estimation methods in most samples. The Kappa test revealed moderate agreement between the two methods (κ = 0.51, P < 0.001), while McNemar’s test showed no significant difference in classification outcomes (P = 0.33), indicating overall consistency with minor discrepancies.

To evaluate the impact of varying renal function levels, CrCL values of 30, 60, 90, 120, and 150 mL. min^-1 were used in the simulation. Since PLT levels in the modeling data are concentrated in the low and normal PLT ranges, the simulation used PLT at the quartiles of the modeling data: 85·10⁹·L^-1 and 266·10⁹·L^-1. Simulation results were presented in Figure 5 and Table 3. Forty different clinical scenarios were simulated. At MIC ≤ 0.5 mg·L^-1, PTA ≥ 90% was achieved in all groups except the 50 mg q12h maintenance dose group. PTA decreased as CrCL increased for the same maintenance dose and PLT level. Compared with maintenance dose and CrCL levels, PLT levels seem to exert a relatively minor influence on the PTA and proportion of AUC_ss,24h within the therapeutic window. In the 50 mg q12h group, the proportion of AUC_ss,24h within the therapeutic window decreased as CrCL increased, whereas an opposite trend was noted in the 100 mg and 125 mg groups. Overall, the 75 mg and 100 mg q12h groups had a slightly higher (75 mg 58.13%; 100 mg 46.78%; 50 mg 42.40%; 125 mg 14.66%) proportion of AUC_ss,24h within the therapeutic window compared to other groups.