A Multi-Center Study of Population Pharmacokinetics of Polymyxin B in Critically Ill Patients

Introduction

Polymyxin B (PB), a peptide antibiotic, is frequently used to combat infections caused by gram-negative bacteria that have developed resistance to multiple drugs, particularly Pseudomonas aeruginosa, Acinetobacter baumannii, and Klebsiella pneumoniae.¹ The global escalation of bacterial resistance, exacerbated by the widespread use of broad-spectrum antibiotics and suboptimal infection control practices, has necessitated the reintroduction of PB into clinical practice.^2–6 This is particularly evident in intensive care units (ICUs), where PB is commonly used to treat multidrug-resistant gram-negative bacterial infections, and optimizing its pharmacokinetics (PK) can significantly enhance therapeutic outcomes through precision dosing.^7–12

PB is administered as a sulfate salt and exhibits considerable plasma protein binding, approximately 50%–60% in general, increasing to over 75% in critically ill patients.¹³ It also demonstrates limited penetration across the blood-brain barrier.¹⁴ As a concentration-dependent antibiotic,^15,16 its efficacy is best predicted by the ratio of the area under the concentration-time curve over 24 hours (AUC₂₄) to the minimum inhibitory concentration (MIC).¹⁷ However, the elimination pathways of PB are not fully elucidated. Although non-renal clearance (CL) is considered predominant, the influence of renal function on PB PK remains controversial.^13,18–20 Previous studies identified body weight and creatinine clearance (CrCL) as significant covariates affecting PB CL.^3,21–25 Notably, existing pharmacokinetic studies are predominantly single-center and involve limited sample sizes, highlighting the need for multicenter population pharmacokinetic (PPK) studies to better characterize the high interpatient variability (often >30%).²⁶

A recent multicenter study by Hanafin et al developed a PPK model for PB in 142 hospitalized patients, identifying body weight as a covariate for volume of distribution and CrCL as a statistically significant covariate but clinically modest for CL.²⁷ Their findings suggested that fixed dosing may be preferable to weight-based dosing for patients weighing 45–90 kg and reported significantly increased CL in patients receiving continuous veno-venous hemodiafiltration (CVVHDF).

In contrast to the broader inpatient population studied by Hanafin et al, our study specifically focuses on critically ill patients in the ICU, a cohort with distinct pathophysiological and management challenges. Furthermore, we employ the CKD-EPI equation to estimate the glomerular filtration rate (eGFR) as a potential covariate, exploring its influence alongside other factors in a two-center design. This study aimed to develop a PPK model for PB in this specific population, identify covariates significantly influencing drug exposure, and propose individualized dosing regimens through Monte Carlo simulations. The goal is to provide refined, clinically applicable dosing guidance tailored to critically ill patients with varying renal function and pathogen MICs.

Materials and Methods

Study Design

The study included critically ill adults (≥18 years old) receiving intravenous polymyxin B sulfate (Shanghai First Biochemical Pharmaceutical Company Limited, Polymyxin B for Injection, Shanghai, China) following at least the third dose. Blood samples were collected between August 2020 and October 2022 at the Guangdong Provincial People’s Hospital and Maoming People’s Hospital. The experimental procedures adhered to the relevant standard operating protocols. The patients or their guardians were fully informed about the research details, and written informed consent was obtained. Patients were excluded if they were pregnant, had allergies or intolerance to PB, or were missing essential data, such as weight or renal function indices.

Basic information (patient ID, sex, age, weight, etc)., PB therapy (sampling time, administration time, dosage, etc)., and laboratory data were gathered from electronic medical records. Laboratory data included total protein (TP), albumin (ALB), total bilirubin (TBIL), direct bilirubin (DBIL), alanine aminotransferase (ALT), aspartate aminotransferase (AST), white blood cell count (WBC), platelet count (PLT), blood urea nitrogen (BUN), serum creatinine (Scr), and uric acid levels. Additionally, the use and modality of continuous renal replacement therapy (CRRT) and the use of extracorporeal membrane oxygenation (ECMO) were recorded. CrCL was calculated using the Cockcroft-Gault formula, and eGFR was determined using the CKD-EPI formula.²⁸ Body mass index (BMI) was calculated as weight in kilograms divided by the square of height in meters (kg/m²).

Polymyxin B Administration and Sample Collection

Enrolled patients received intravenous infusions of PB sulfate. The prescribing practice generally followed the product label recommendation of 1.5–2.5 mg/kg per day (1 mg = 10,000 IU), divided into two doses. In our cohort, specific regimens varied based on the treating clinician’s assessment of the individual’s clinical status, as dose adjustments were not based on CRRT status or renal function. This included the use or omission of a loading dose, and a range of fixed maintenance doses. For instance, observed regimens included a 100 mg loading dose followed by a 50 mg maintenance dose every 12 hours, or maintenance doses ranging from 50 mg to 100 mg every 12 hours without a preceding load. Using EDTA-K₂ tubes, blood samples were collected at seven time points: pre-dose (10 min before administration); at 5 min, 1 h, 2 h, 4 h, and 8 h after the end of the infusion; and pre-next dose (10 min before the subsequent dose). The collected blood samples were then centrifuged at 996 × g for 10 min at 4 °C, after which the separated plasma was frozen at –80 °C for subsequent analysis. The plasma concentration of PB was quantified using high-performance liquid chromatography-mass spectrometry (HPLC-MS/MS). The analytical method employed for this purpose was developed by our team and has been previously described.²⁹

Population Pharmacokinetics Modeling

A PPK model was developed using nonlinear mixed-effects modeling (NONMEM) software (version 7.3.0, Icon Inc, PA, USA). As previous studies indicate that PB disposition has been described by both one- and two-compartment models,³⁰ we established both one- and two-compartment base models and assessed their fit using goodness-of-fit (GOF) plots and parameter evaluation. The model with the best fit was selected as the base model. Model parameters were estimated using the first-order conditional estimation method (FOCE-I), accounting for both inter-individual and residual error. Interindividual variability in parameters was modelled using an exponential error model (Eq.1):

(1)

Where represents the individual parameter value, is the population typical value of the parameter, and is the interindividual random effect, which follows a normal distribution with a mean of 0 and a variance of ω². The residual error was modelled using a mixed proportional and additive error model (Eq.2):

(2)

Where is the observed concentration, is the model-predicted concentration, and and are the proportional and additive residual random errors, respectively, each normally distributed with a mean of zero and variances of and .

Following the development of the base model, the impact of covariates was evaluated by comparing the decrease in the objective function value (OFV) between the base model and models incorporating covariates. The influence of potential covariates was examined through forward inclusion and backward elimination processes. Continuous covariates included age, weight, BMI, ALT, AST, TP, ALB, TBIL, DBIL, BUN, eGFR, WBC, PLT, uric acid, Scr, and CrCL. Categorical covariates were sex, CRRT status, and ECMO status. The significance level for covariate inclusion in the forward selection was set at p < 0.05 (ΔOFV of −3.84, df = 1), and for backward elimination, it was set at p < 0.001 (ΔOFV of −10.828, df = 1). The effect of continuous and categorical covariates on parameters was tested using the following equations (Eq.3 and Eq.4, respectively):

(3)

Where P_ij is the j-th pharmacokinetic parameter of the i-th patient. P_{tv, j} is the population typical value, and θ_j is the fixed effect. η_j represents inter-individual variability and COV_median is the median value of covariate. CRRT status was recorded prospectively at the time of each blood sample collection. To ensure accuracy, these records were subsequently cross-verified against the patients’ electronic medical records and nursing charts, which documented precise CRRT start and stop timestamps. This allowed us to confirm the CRRT status (on/off) at each pharmacokinetic sampling time and to calculate the total duration of CRRT exposure. In the NONMEM dataset, a binary covariate (CRRT = 1 if CRRT was ongoing at the exact sampling time; CRRT = 0 otherwise) was created and time-aligned with each pharmacokinetic sample. This covariate was tested in the model using a categorical covariate model structure (Eq.4).

Model Evaluation and Validation

The model was evaluated using bootstrap, GOF, prediction-corrected visual predictive checks (pcVPC), and normalized prediction distribution errors (NPDE).

A bootstrap analysis was performed with 1000 replicates using Perl Speaks NONMEM (PSN).³¹ Model parameters were re-estimated for each bootstrap dataset using NONMEM. Model robustness was assessed by the success rate of parameter estimation, which was expected to exceed 90%. The median and 95% confidence intervals (CI, 2.5%-97.5% percentiles) from the bootstrap analysis were compared to the final model’s parameter estimates, which should fall within the 95% CI.

GOF was assessed by plotting observed concentrations (DV) against population predictions (PRED) and individual predictions (IPRED). The identity line (y = x) and a locally weighted regression (LOESS) curve were superimposed on these plots to visualize trends. Plots of conditional weighted residuals (CWRESI) versus population predictions (CWRESI vs PRED) and versus time after the first dose (CWRESI vs Time) were examined to detect any systematic bias.

Additionally, NPDE and pcVPC (based on 1000 simulations) were used to evaluate the normality of the prediction distribution errors and the overall model fit, respectively. Together, these methods provided a comprehensive evaluation of the model’s predictive performance for PB pharmacokinetics in the studied population.

Simulations

For simulation purposes, the eGFR covariate was categorized into five levels to reflect varying renal function: <15, 15–30, 30–60, 60–90, and 90–130 mL/min/1.73 m². Pharmacokinetic simulations were performed to evaluate steady-state drug exposure on the third day. Different dosing regimens were established, and their probability of target attainment (PTA) was assessed. PTA was calculated as the percentage of simulated patients achieving the target value. The threshold for efficacy was set at PTA ≥ 80%. We conducted Monte Carlo simulations for a range of dosing regimens administered every 12 hours (q12h), incorporating both first dose (FD) and maintenance dose (MD) starting at a minimum of 50 mg. The dosage was then systematically increased in increments up to 150 mg. If the target PTA was not achieved with a 150 mg maintenance dose, the dose was further increased or the dosing interval was shortened. Monte Carlo Simulations based on 1000 replicates per eGFR stratum ensure stability of PTA estimates. According to reported studies, antimicrobial efficacy was characterized using the AUC₂₄/MIC ratio, with a target value of >50 established as the efficacy threshold.^13,21,32 The pharmacodynamic index of AUC₂₄/MIC ≥ 50 was applied to evaluate the simulation results, and MICs were from 0.125 to 4 mg/L.

Results

Population PK Model

A total of 56 patients contributed 350 blood samples for the analysis. Patient demographics and clinical characteristics are summarized in Table 1 and Supplementary Table 1. Two patients received ECMO support. Twenty patients (35.7%) received CRRT, exclusively via the continuous veno-venous hemofiltration (CVVH) modality, with a blood flow rate of 150–200 mL/min and a replacement fluid flow rate of 1500–2000 mL/h. The effluent rate ranged from 26.5 to 49.8 mL/kg/h. An exploratory analysis found no significant correlation between effluent rate and drug clearance in patients receiving CVVH (absolute r = 0.158).

Table 1 Characteristics of Patients

Comparative analysis of the model fit demonstrated superior GOF for the two-compartment model compared to the one-compartment model. During forward inclusion of covariate screening for CL, ECMO (ΔOFV −7.258), ALT (ΔOFV −12.057), BUN (ΔOFV −9.301), CrCL (ΔOFV −5.705), and eGFR (ΔOFV −12.658) were identified as statistically significant covariates (p < 0.05). The model incorporating eGFR, which resulted in the largest reduction in OFV, was selected as the base model for subsequent covariate testing. Additional covariates were iteratively included, resulting in a full model containing three covariates: eGFR, ECMO, and ALT. During backward elimination, ECMO and ALT were sequentially removed. The final model retained eGFR as the sole significant covariate influencing CL, with final parameter estimates provided in Table 2. The elimination CL, volume of the central compartment (V₁), volume of the peripheral compartment (V₂), and the intercompartmental clearance (Q) were determined to be 1.68 L/h, 14.30 L, 48.84 L, and 4.67 L/h, respectively. The shrinkage of the inter-individual variation in CL was 4%. Weight (ΔOFV −0.466, p > 0.05) and CVVH (ΔOFV −0.137, p > 0.05) did not demonstrate a statistically significant impact on the model.

Table 2 The Population Pharmacokinetic Parameter Values and Bootstrap of the Final Model

The final model was as follows:

Model Evaluation

GOF plots (Figure 1) and pcVPC (Figure 2) indicated that the final model adequately described the observed data. In Figure 1A and B, data points were tightly distributed around the y=x reference line without systematic deviation, indicating no significant concentration-related bias between predicted and observed values. And in Panels C and D, CWRESI values were randomly scattered around the 0-baseline dashed line, with no trends of offset as PRED changes or time elapses. This confirmed the model exhibits neither concentration-dependent nor time-dependent prediction bias. The bootstrap analysis (Table 2) showed a success rate of 93.5%. NPDE analysis (Figure 3) indicated no significant prediction bias in the final model. Most errors were distributed within the CI, and the 5th, 50th, and 95th percentiles of the observed values were within the 95% prediction interval.

Figure 1 Diagnostic plots of final model. (A) Observed polymyxin B concentrations versus population predicted (PRED) concentrations; (B) Observed polymyxin B concentrations versus individual predicted (IPRED) concentrations; (C) conditional weighted residuals with inter- and intra-subject variability interaction (CWRESI) versus population predicted concentrations; (D) conditional weighted residuals with inter- and intra-subject variability interaction (CWRESI) versus time. The solid black line represents the line of unity ((A and B): y=x; (C and D): y=0) and the red dashed line represents a locally weighted smoothing (loess) curve.

Figure 2 Prediction-corrected visual predictive check (pcVPC) for the final model based on 1000 replicates of Monte Carlo simulations. Blue circles represent the observed concentrations. The red solid line represents the median of the observed values, and the red dashed lines represent the 5th and 95th percentiles. The red shaded band represents the 95% confidence interval for the median of the model predictions, and the blue shaded bands represent the 95% confidence intervals for the 5th and 95th percentiles of the predictions.

Figure 3 Normalized prediction distribution errors (NPDE) for the final model. (A) Quantile-quantile plot of NPDE. (B) Distribution of NPDE. (C) NPDE versus time after last dose. (D) NPDE versus predicted concentrations. The observed concentrations are shown as filled circles, and solid lines represent the 5th, 50th, and 95th percentiles of observed data. Red- or blue-shaded areas represent the 95% prediction interval.

Simulation

Simulations (Figure 4) identified the dosing regimens required to achieve the target PTA (≥80%) across a range of MIC values and renal function levels. Simulations output steady-state plasma concentrations corresponding to the third day of therapy. As expected, PTA demonstrated a dose-dependent increase across all simulated scenarios. The analysis indicated that higher loading and maintenance doses were essential to overcome higher MICs and the increased drug clearance associated with better renal function. Based on this analysis, dosing recommendations stratified by eGFR and MIC are provided in Table 3. The key trends can be summarized as follows: for any given eGFR category, successively higher FD/MD combinations were required as the MIC increased. Conversely, for a fixed MIC, patients with higher eGFR (≥60 mL/min/1.73m²) consistently required more intensive dosing regimens to achieve the target PTA compared to those with impaired renal function. The recommended regimens primarily use a 12-hour dosing interval and aim to achieve the pharmacodynamic target with the lowest effective dose to minimize cumulative drug exposure.

Table 3 Recommended Polymyxin B Dosing Regimens Stratified by eGFR and MIC to Achieve PTA ≥80%. Simulations Based on 1000 Replicates per eGFR Stratum

Figure 4 Probability of target attainment (PTA) for achieving AUCss/MIC ≥ 50 at steady state on day 3 after a one-hour infusion, stratified by estimated glomerular filtration rate (eGFR) and minimum inhibitory concentration (MIC). The results for the recommended dosing regimens are shown, with the efficacy threshold set at PTA ≥ 80% (indicated by the dashed horizontal line). PTA represents the proportion (%) of simulated patients achieving the target. Simulations were performed with 1000 replicates per regimen.

Discussion

In this study, we developed a two-compartment PPK model for PB in critically ill patients and identified eGFR as the sole covariate significantly influencing PB CL. Our findings offer specific dosing recommendations tailored for patients with varying degrees of renal function and different pathogen MICs.

Consistent with previous reports, our results supported the use of a two-compartment model to describe PB PK.³⁰ Although earlier studies identified body weight and CrCL as covariates,^22–24,33 our analysis revealed that eGFR, estimated via the CKD-EPI equation, provided a superior fit (ΔOFV −12.658, p < 0.001) compared to CrCL (ΔOFV −5.705, p < 0.05).

Our findings corroborate and substantially extend the work of Hanafin et al.²⁷ Both studies support a two-compartment model for PB and confirm that renal function is a statistically significant covariate for clearance, whereas body weight is not. However, our work provides critical advancements. While Hanafin et al concluded that the influence of CrCL was statistically significant but clinically modest, our analysis demonstrates that the incorporation of the CKD-EPI-derived eGFR is clinically consequential. While their study was conducted mainly in Western centers, it enrolled a cohort with a relatively low proportion of patients receiving extracorporeal support. In contrast, our cohort was derived from ICUs of two Chinese hospitals and included a higher proportion of CVVH cases (20/56, 35.71%) as well as some patients on ECMO. Furthermore, Hanafin et al reported significantly increased clearance in patients receiving CVVHDF, whereas our study demonstrated that CVVH had no significant effect on drug clearance. This clear discrepancy underscores that the specific modality of renal replacement therapy is a critical determinant that cannot be generalized. Unlike Hanafin et al who identified body weight as a covariate for the volume of distribution, our final model did not.

ECMO and ALT showed statistical significance during the initial covariate screening and were excluded during backward elimination. Studies have reported elevated PB exposure in ECMO patients, suggesting that ECMO may alter its PK through mechanisms such as drug adsorption and changes in volume of distribution and clearance, making it a potential covariate in PB PK.³⁴ However, in our study, although ECMO showed significance during preliminary screening, it was not retained in the final model during backward elimination. This was likely due to the statistical bias introduced by the very limited number of ECMO patients (n=2) in our cohort. Therefore, this study could not robustly demonstrate an independent effect of ECMO on PB CL, warranting validation in a larger population of ECMO patients. During the covariate screening for factors influencing V₁, no statistically significant covariate was identified, including ALB (ΔOFV −1.581, p > 0.05).

Among the reported simulation schemes, most have designed dosing regimens based on body weight,^22,35 with a smaller subset dosing according to CrCL.^24,33 The dosing regimens proposed in our simulations are grounded in achieving a PK/PD target of AUC₂₄/MIC ≥ 50. Our results demonstrated that the required doses increased with increasing MIC and better renal function (higher eGFR), necessitating personalized dosing strategies. While Hanafin et al suggested fixed dosing may be appropriate for many patients (weight 45–90 kg) with non-pulmonary infections and MIC ≤ 2 mg/L, our study provided a more granular, eGFR-stratified dosing matrix for the critically ill population, which may be particularly valuable in settings with high MIC pathogens or fluctuating renal function.²⁷

Several limitations warrant acknowledgment. First, the sample size (n=56), although reasonable for a PPK study, was moderate, and the inclusion of only two centers may limit generalizability. Furthermore, the study was not powered to adequately assess the influence of rare conditions, such as ECMO or specific organ dysfunctions.

Conclusion

We developed a PPK model for PB in a Chinese critically ill population, identifying eGFR as a key covariate for clearance. This model enabled the creation of a dosing framework stratified by renal function and MIC, delivering actionable guidance for individualized therapy.