Population Pharmacokinetics and Limited Sampling Strategy for Therapeutic Drug Monitoring of Polymyxin B in Chinese Patients With Multidrug-Resistant Gram-Negative Bacterial Infections

Polymyxin B is used as a last therapeutic option for the treatment of multidrug-resistant Gram-negative bacterial infections. This study aimed to develop a population pharmacokinetic model and limited sampling strategy, a method to estimate the area under the concentration curve (AUC) by using a limited number of samples, to assist therapeutic drug monitoring of polymyxin B in Chinese patients. Population pharmacokinetic analysis was performed using Phoenix® NLME with data obtained from 46 adult patients at steady state. Various demographic variables were investigated as potential covariates for population pharmacokinetic modeling. The limited sampling strategies based on the Bayesian approach and multiple linear regression were validated using the intraclass correlation coefficient and Bland-Altman analysis. As a result, the data was described by a two-compartment population pharmacokinetic model. Through the modeling, creatinine clearance was found to be a statistically significant covariate influencing polymyxin B clearance. The limited sampling strategies showed the two-point model (C0h and C2h) could predict polymyxin B exposure with good linear relativity (r2 > 0.98), and the four-point model (C1h, C1.5h, C4h, and C8h) performed best in predicting polymyxin B AUC (r2 > 0.99). In conclusion, this study successfully developed a population pharmacokinetic model and limited sampling strategies that could be applied in clinical practice to assist in therapeutic drug monitoring of polymyxin B in Chinese patients.

Polymyxin B is used as a last therapeutic option for the treatment of multidrug-resistant Gram-negative bacterial infections. This study aimed to develop a population pharmacokinetic model and limited sampling strategy, a method to estimate the area under the concentration curve (AUC) by using a limited number of samples, to assist therapeutic drug monitoring of polymyxin B in Chinese patients. Population pharmacokinetic analysis was performed using Phoenix ® NLME with data obtained from 46 adult patients at steady state. Various demographic variables were investigated as potential covariates for population pharmacokinetic modeling. The limited sampling strategies based on the Bayesian approach and multiple linear regression were validated using the intraclass correlation coefficient and Bland-Altman analysis. As a result, the data was described by a two-compartment population pharmacokinetic model. Through the modeling, creatinine clearance was found to be a statistically significant covariate influencing polymyxin B clearance. The limited sampling strategies showed the twopoint model (C 0h and C 2h ) could predict polymyxin B exposure with good linear relativity (r 2 > 0.98), and the four-point model (C 1h , C1 .5h , C 4h , and C 8h ) performed best in predicting polymyxin B AUC (r 2 > 0.99). In conclusion, this study successfully developed a population pharmacokinetic model and limited sampling strategies that could be applied in clinical practice to assist in therapeutic drug monitoring of polymyxin B in Chinese patients.
Keywords: polymyxin B, population pharmacokinetics, limited sampling strategy, therapeutic drug monitoring, multidrug-resistant Gram-negative bacterial infection INTRODUCTION Polymyxin B, a polypeptide antibiotic obtained from the fermentation products of Paenibacillus polymyxa, was first used clinically in the 1950s but was largely abandoned in the 1970s due to nephrotoxicity and neurotoxicity. Recently, with the prevalence of multidrug-resistance (MDR) Gram-negative bacterial infections, polymyxin B has been re-introduced to clinical practice against these challenging infections. It is currently considered as the last drug resort for MDR Gramnegative bacterial infections (Velkov et al., 2013;Nation et al., 2015).
Despite being used clinically for decades, few pharmacokinetic/ pharmacodynamics (PK/PD) data is available for polymyxin B (Tam et al., 2005;Bergen et al., 2012;Tsuji et al., 2016;Lakota et al., 2018;Landersdorfer et al., 2018). These studies found the area under the concentration-time curve to minimum inhibitory concentration ratio (fAUC: MIC), as the PK/PD index, appeared to be a good efficacy predictor for polymyxin B. Accordingly, a global multidisciplinary panel of infectious diseases and pharmacotherapy experts recommended therapeutic drug monitoring (TDM) was required to optimize the clinical use of polymyxin B. An AUC across 24 hours at steady state (AUC ss,24 hr ) target of 50 to 100 mgcḣ/L, corresponding to a C ss,avg of 2 to 4 mg/L, might be acceptable for polymyxin B therapy (Tsuji et al., 2019). However, the guidelines also noted that data were lacking for PK properties and AUC ss,24 hr targets of polymyxin B, and the recommended therapeutic target largely based on in vitro and animal data rather than clinical outcome (Alosaimy et al., 2019). In China, polymyxin B is only commercially available since 2017, there is almost no PK data on polymyxin B from Chinese patients. Therefore, it is necessary to assess the population PK of polymyxin B in Chinese patients.
In addition, to calculate AUC, it requires multiple blood draws during a dosing interval to get a full pharmacokinetic curve. This is not only time-consuming and expensive but also a huge burden to patients, and therefore unfeasible in clinical practice. To solve this problem, the limited sampling strategy (LSS), estimating AUC with one or a few samples, would be suitable in the clinic (David and Johnston, 2000). It has been proposed for TDM of many drugs, such as mycophenolate mofetil and colistin (Zhang et al., 2018;Kim et al., 2019;Van Der Galiën et al., 2020). However, no LSS report is available for polymyxin B at present. Therefore, this study aimed to develop and validate a population PK model, along with clinically feasible LSS for TDM of polymyxin B in Chinese patients with MDR Gramnegative bacterial infections.

Patients Demographics
A single-center clinical trial was conducted between April 2018 and November 2019 at the First Affiliated Hospital of Zhengzhou University. The study protocol was approved by the hospital ethics committee review board, and all participating subjects signed the informed consent (Zhengzhou University Medical Research and Ethics Committee, No. L2018K129).

Polymyxin B Administration and Sample Collection
According to the Chinese package insert, the maintenance dose of polymyxin B was 50 to 100 million units (1 million units equal to 1 mg) twice daily. The loading dose was 100 to 150 million units in clinical practice. The recommended infusion time was one hour but could be lengthened as needed. Polymyxin B treatment including dose, infusion time, and duration of therapy depended on the medical teams.
Blood samples were obtained after 3 days of therapy. On day 4, one blood sample (2 ml, C 0h ) was collected immediately at predose, and five to seven blood samples (mainly C 0.5h , C 1h , C 1.5h , C 2h , C 4h , C 6h , and C 8h ) were collected at pre-next dose into EDTA tubes for each patient. To calculate AUC 0-12h , C 0h was regarded as C 12h for those who hid not obtained C 0h . All samples were centrifuged at 3,500×g for 10 min. The supernatant was collected and stored at −80°C until analysis.

Quantification of Polymyxin B Concentrations
Plasma concentrations of polymyxin B were determined using liquid chromatography-tandem mass spectrometry (LC-MS/MS) in the hospital laboratory. The method was previously described by our group (Wang et al., 2020). In brief, the calibration curves showed acceptable linearity over 0.2 to 10 µg/ml for polymyxin B1and 0.05 to 2.5 µg/ml for polymyxin B2. The upper limit of quantification was extended to 20 µg/ml for polymyxin B1 and 5.0 µg/ml for polymyxin B2 after four-fold dilution. For precision and accuracy, the intra-and inter-day imprecision of polymyxin B1 and B2 was less than 13.93%, and the inaccuracy ranged from −10.87 to 11.13%. Since polymyxin B1 and B2 had similar structures, molecular weight, pharmacological activities, and pharmacokinetic characteristics, the plasma concentration of polymyxin B was summed to derive the total polymyxin B1 and B2 concentrations (Tam et al., 2011;Manchandani et al., 2016a).

Population Pharmacokinetics Analysis
Based on a previously published study (Avedissian et al., 2019), the PK parameters were determined by one-and twocompartmental models performed using Phoenix ® NLME software (v7.0, Pharsight, Mountain View, CA, USA). Initial PK parameters were estimated by the Naive-pooled model. PK models were estimated by the first-order conditional where P i is the PK parameters of the i th patient, q is the population PK pharmacokinetic parameters, and h i is a normally distributed random variable with a mean of 0 and a variance of w 2 . Intra-individual variability (residual error) is described using additive (C obs = C pred + ϵ), proportional [C obs = C pred × (1+ ϵ)] or mixed (additive + proportional) models, where C obs and C pred are the observed and predicted concentrations, and ϵ is an error variable with a mean of 0 and a variance of s 2 . Model assessment criteria included precision of parameter estimates (standard error), goodness-of-fit plots, and likelihood ratio test (−2LL). For the one-compartment model, basic PK parameters were the volume of central compartment distribution (V) and central compartment clearance (Cl). For the two-compartment model, there were two additional parameters as the volume of peripheral compartment distribution (V2) and inter-compartmental clearance (Q).
In the process of model development, age, sex, body weight, ALT, AST, GGT, ALP, urea nitrogen, Scr, serum uric acid, serum proteins, serum albumin, TBIL, DBIL, and CrCL were evaluated as the covariates. The covariates selection was evaluated using a stepwise process (a forward-selection process and then a backward-elimination process). By comparing with initial model, the inclusion criteria for covariates was a drop > 3.84 (P = 0.05) of objective function (OFV; −2LL) during forwardselection. Then each selected covariate was re-evaluated by backward-elimination. An increase of OFV > 6.63 (P = 0.01) was required for confirmation.
After covariate selection, the correlations between population PK parameters were graphically examined in the scatter plot. The relevant parameters were introduced to non-diagonal random effects. A drop of OFV > 6.63 (P = 0.01) was retained to obtain the final model.
Model validation was assessed by using plots of observed concentrations (DV) versus population predicted concentrations (PRED) or individual predicted concentrations (IPRED), and conditional weighted residuals (CWRES) versus time (IVAR) or PRED. Additionally, the model performance was evaluated by a prediction-corrected visual predictive check (VPC). The 5 th , 50 th , and 95 th quantiles and the 90% confidence intervals (CIs) of 5 th , 50 th , and 95 th quantiles (200 replicates) were developed from the final model parameters. The model stability was assessed using bootstrap analysis (1,000 bootstrap samples). The median, standard error (SE), and coefficient of variation (CV%) were calculated from the empirical bootstrap distribution and compared estimates with the original dataset.

Model-Based Simulation
For optimal dose selection, the plasma concentration-time profile of 1,000 individuals was simulated using the final population PK model. The following dose regimens were evaluated: 100 mg loading dose with 50 mg maintenance dose twice daily, 150 mg loading dose with 75 mg maintenance dose twice daily, and 150 mg loading dose with 100 mg maintenance dose twice daily. The infusion rate was set as 50 mg/h. The covariate of CrCL was selected as 31.3, 105.9, and 315.2 ml/min. which were the low, medium, and high values from the sample of patients.

Limited Sampling Strategies
The AUC 0-12h of polymyxin B was calculated by the linear-uplog-down method. LSSs were investigated with Bayesian analysis and linear regression analysis. By using the Bayesian approach, Polymyxin B predicted concentrations were estimated from limited point values with the final population PK model. The observed and predicted AUC 0-12h were analyzed using the Pearson correlation coefficient (r). The linear regression analysis was performed on IBM SPSS Statistics (v24.0, SPSS Inc., Chicago, USA). The stepwise forward multiple regression method was used to investigate the correlation between AUC 0-12h and polymyxin B concentrations at different time points. The determination coefficient (r 2 ) evaluated the regression level of the equation. Good regression equations were selected for model validation.
The model was internally validated by the Jackknife method. One sample was removed from the original sample at a time, then regression analysis of the remaining sample was performed to obtain a new equation and calculate the AUC 0-12h of the removed one sample. Prediction error (PE, Equ. 2) and root mean square error (RMSE, Equ. 3) were used to evaluate the accuracy and precision of the model. The acceptance criteria of RMSE < 15 and r 2 > 0.95 (Van den Elsen et al., 2018).
The consistency between the measured and predicted value of AUC 0-12h was evaluated by intraclass correlation coefficient (ICC), and Bland-Altman (BA) analysis. The BA analysis was expressed by the scatter diagram and limits of agreement. The lower limit of 95% CI < 0.9 and the limits of agreement within ±15% were acceptable for the clinic (Zhang et al., 2018). The best model was selected based on the values of r 2 , PE, RMSE, ICC, and BA analysis.

Patients
A total of 46 patients contributed to 331 plasma samples were enrolled in this study. The patients' demographic and clinical information were listed in Table 1. Among them, pathogenic bacteria cultures results showed that sputum was the most common primary site of infection, along with blood, cerebrospinal fluid, and puncture fluid, etc. MDR Acinetobacter baumannii and Klebsiella pneumoniae were the major causative agent of the infections, and 6 patients were infected with two kinds of MDR Gram-negative bacterial.

Population PK Model
Preliminary analysis of base model showed the OFVs of one-and two-compartmental models were 1062.13 and 694.63, respectively. Based on OFV, CV values, and diagnostic scatter plots, a two-compartment model with a proportional option was chosen as the base model. In the next step, covariate model building identified CrCL as the covariate for Cl (DOFV = 8.26, P < 0.01). Age, gender, and other laboratory data had no significant effect on population PK parameters. Furthermore, a strong correlation between Cl, Cl2, and V was observed, and then incorporated it into non-diagonal random effects (DOFV = 30.60, P < 0.01). The correlations between V-Cl, V-V2, and Cl-V2 were 0.713, 0.667, and 0.571, respectively. Accordingly, the final PK model was shown in Equ. 4 to 7, where 105.9 ml/min was the median of CrCL.
The goodness-of-fit plots for the final model were shown in Figure 1. The observed concentrations were consistent with PRED and IPRED, and the plots of CWRES vs time and PRED were normally distributed. The estimated covariates and bootstrap replicates were shown in Table 2, which indicated the final model had qualified stability. In prediction corrected-VPC (Figure 2), most of the observed plots were distributed within the 90% CIs of predicted corresponding quantiles, which indicated the prediction of simulated data matched the observed plots.

Model-Based Simulation
The median simulated plasma concentration-time profiles of polymyxin B based on different doses and CrCL values were displayed in Figure 3, and the simulated AUC 24h and C ss,avg of polymyxin B on day four was quantified in Table 3.

Limited Sampling Strategy
The AUC 0-12 h of 46 patients was 43.64 ± 27.68 mg·h/L with a range of 8.50 to 122.84 mg·h/L. In this study, the infusion time of 37 patients was one hour, and that of others was half an hour or two hours. To ensure the accuracy of the results, only 37 patients contributed 275 plasma samples were enrolled in LSS analysis.
Ten best-performing strategies using the Bayesian approach were displayed in Table 4. The steady-state trough concentration (C 0h ) and peak concentration (C 1h ) showed a medium correlation with AUC 0-12h (r 2 = 0.862 and 0.846). While C 2h and C 4h both had a good correlation with AUC 0-12h (r 2 = 0.922 and 0.976). For the two, three, and four-point samples, all models showed a good agreement with the AUC 0-12h (r 2 ≥ 0.98). With the increase of sample points, the PE% range was narrowed down and RMSE value became small, which was consistent with the ICC and limits of agreement results. In BA analysis plots of models 6 and 10 (Figure 4), all points were within ±10% of the limits of agreement.
The simple and multiple regression equations for estimating AUC 0-12h with good correlation r 2 were shown in Table 5. C 2h model and all models using two to four samples met the acceptance criteria (Van den Elsen et al., 2018;Zhang et al., 2018). With the increase of sample points, the linear correlations between predicted and measured AUC 0-12h were improved. In BA analysis plots of models 16 and 20 (Figure 4), all points were within ±10% of the limits of agreement. Based on the above analysis, all two, three, and four-point models met the criteria for acceptability in the clinic, and the four-point strategy was the best LSS.

DISCUSSION
To date, few studies have been conducted on population PK of polymyxin B (Kwa et al., 2008;Sandri et al., 2013;Kubin et al., 2018;Manchandani et al., 2018;Miglis et al., 2018). To our best knowledge, this study included the largest samples and most covariates from acutely-ill patients, aimed at evaluating the population PK of polymyxin B and identifying potential factors influencing the PK variability. The result indicated a twocompartment model adequately described population PK of polymyxin B in Chinese patients. Although there were studies performed using the one-compartment model (Kwa et al., 2008;Kubin et al., 2018;Manchandani et al., 2018), which may be due to limited blood points collected at distribution and elimination phases. This study showed that the median polymyxin B CL was 1.786 L/h, which was less than those reported previously (range, 1.87-2.5 L/h). For V estimate, the median was 6.218 L, closed to 6.35 L reported by Sandri et al. (2013), and was quite different from 20.39 to 33.77 L reported by Avedissian et al. (2018) and  Thamlikitkul et al. (2017). With regard to V2 and Q estimates, the values were variant in different literature, probably because variability in population PK parameter estimates was high, with CV% often >30% (20.6-73.3%) (Sandri et al., 2013;Avedissian et al., 2018;Miglis et al., 2018;Avedissian et al., 2019). The impact of clinical variables on population PK models of polymyxin B was always inconsistent, especially the influence of body weight and CrCL on polymyxin B clearance. FDA package inserts suggested polymyxin B dose reduction in patients with renal insufficiency, which was inconsistent with the recommendation provided by international consensus guidelines for the optimal use of polymyxins (Li et al., 2019;Tsuji et al., 2019). With a low correlation coefficient (Supplementary Figure 1), this study found CrCL was a statistical covariate of polymyxin B clearance (Equ. 6), which was in agreement with other studies Manchandani et al., 2018). However, Manchandani et al. (2018)   indicated the effect was clinically insignificant. Furthermore, by comparing the clearance and exposure of polymyxin B in normal renal function patients (CrCL ≥ 80 ml/min) and renal insufficient patients (CrCL < 80 ml/min), Thamlikitkul et al. (2017) reported the exposure of polymyxin B was similar between two groups (P = 0.8) after standardizing AUC for the daily dose. The clearance value of renal insufficient group was lower than that of normal function group (2.0 L/h vs. 2.5 L/h, P = 0.06), but there was no statistical difference. This might be because renal clearance of polymyxin B was only a minor elimination pathway in both critically ill patients and animals (Zavascki et al., 2008;Abdelraouf et al., 2012;Manchandani et al., 2016b). Taken together, given the limited cases, a small number of samples collected from each patient, and different compartment models, most of the studies, including ours, found a borderline effect of CrCL on polymyxin B clearance, the correlation needed further investigation. As for body weight, no correlation between total body weight and clearance was observed in this study. This might be because the body weight range (45-98 kg) was narrow. With a wide range of body weight (41-250 kg), Sandri et al. (2013) demonstrated that total body weight influenced the PK parameters of polymyxin B and total body clearance, and suggested loading and maintenance doses of polymyxin B were best scaled by total body weight. Further research reported the adjusted body weight rather than the total body weight may be a better factor in influencing polymyxin B exposure .
Based on the 50 th percentile of simulated AUC 24h , after administration of three dosage regimens, the AUC 24h for patients with normal renal function and renal hyperfunction was ranged from 37.92 to 108.57 mg·h/L, which was quite different from patients with renal insufficiency (87.20-167.90 mg·h/L). Accordingly, for patients with normal renal function and renal hyperfunction, a regime of 75 to 100 mg maintenance dose twice daily would be sufficient to reach the therapeutic window (AUC ss,24 hr of 50-100 mgcḣour/L; Tsuji et al., 2019). While, as for patients with renal dysfunction, the regime of 50 mg maintenance dose twice daily would be a better option. These observations were in agreement with a previous study, although it did not provide a dosing recommendation for patients with renal dysfunction by using a one-compartment model (Manchandani et al., 2018).
Since obtaining a full concentration-time curve to calculate the AUC is not always feasible in the clinic, LSS offers a practical  approach to estimating the AUC. As shown in Tables 4 and 5, the results of single-time point model by Bayesian method were more accurate than that of linear regression analysis, and the results of multiple-time point models obtained by the two methods were similar. In single time-point models, C 4h showed a better ability to predict polymyxin B AUC than C 0h , C 1h , and C 2h in both methods. It was likely that C 0h , C 1h , and C 2h were obtained in elimination, absorption, and distribution phases; while, C 4h was obtained between distribution and elimination phases, which could better predict polymyxin B exposure. In multiple time-point models, all models (Tables 4 and 5) displayed a good agreement with the measured AUC 0-12 h (r 2 ≥ 0.98), which was acceptable for polymyxin B TDM. The 4point models (C 1h , C 1.5h , C 4h , and C 8h ), which included the absorption, distribution, and elimination phases, were the best predictor of polymyxin B AUC 0-12h . Additionally, it was reported a 2-point model including C 0h , and C 2h was appropriate for colistin TDM (r 2 = 0.98) because the C 0h sample had an association with renal toxicity and the C 2h sample was essential to monitor efficacy (Kim et al., 2019). This result was in agreement with that of model 16 in this study. Accordingly, models 5 to 10 and 15 to 20 all can be recommended for  polymyxin B TDM, especially C 0h -C 2h and C 1h -C 1.5h -C 4h -C 8h models. Our study has several limitations. Firstly, this study enrolled a relatively small number of patients, leading to a lack of external validation of the population PK model and LSS. Future studies should evaluate the PK of polymyxin B with a larger population. Second, since patients with diverse underlying conditions, drugdrug interactions, and co-administration of other antibiotics were not included in the population PK model. Third, individual PK parameters may change during therapy occasions. Through the modeling, ignoring inter-occasion variability may lead to a big bias in parameter estimates, especially for drugs with large intraindividual variability that require TDM (Karlsson and Sheiner, 1993;Abrantes et al., 2019). However, since all data were collected on day four, whether this population PK model is suitable for other occasions is still unknown. Fourth, the LSS results were only applicable to patients who intravenously infused polymyxin B for one hour. Finally, the efficacy and toxicity thresholds based on clinical outcomes were not investigated in this study, which was the next work of our group.

CONCLUSION
In conclusion, a two-compartment population PK model was successfully established to characterize the PK parameters of polymyxin B in Chinese patients with MDR Gram-negative bacterial infections. Furthermore, as far as we know, this is the first study to develop and validate the LSS of polymyxin B. The results suggested 2-point model (C 0h and C 2h ) and 4-point model (C 1h , C 1.5h , C 4h , and C 8h ) performed well in predicting polymyxin B AUC, which could be applied in clinical practice to assist TDM of polymyxin B.

DATA AVAILABILITY STATEMENT
All datasets generated for this study are included in the article/ Supplementary Material.

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by Zhengzhou University Medical Research and Ethics Committee. The patients/participants provided their written informed consent to participate in this study.

AUTHOR CONTRIBUTIONS
PW and JY designed the research. QZ performed the experiments. ZZ analyzed the results. MF, TS, and XZ supervised the research and revised the manuscript. All authors approved the final manuscript.