This was a multicenter, retrospective, and observational study of hospitalized adult patients receiving linezolid for confirmed or suspected multiresistant Gram-positive bacterial infections from January 2018 to December 2019 at the First Affiliated Hospital of Wenzhou Medical University and Ruijin Hospital, Shanghai Jiaotong University School of Medicine. The inclusion criteria were as follows: 1) patients ≥18 years with confirmed or suspected multiresistant Gram-positive bacterial infection; 2) patients who received intravenous and/or oral linezolid for at least 10 days; and 3) at least one steady-state concentration of linezolid was collected. The exclusion criteria were as follows: 1) patients who received renal replacement therapy or extracorporeal membrane oxygenation and 2) patients who died within 24 h after being treated with linezolid. The patients who were included for toxicity analysis were further excluded if 1) baseline PLT <75 × 10⁹ cells//L, 2) baseline hemoglobin <6.8 g/dl for males or 6 g/dl for female individuals, 3) baseline absolute neutrophil count <500 cells/μL, and 4) baseline total bilirubin > 5-times the upper limit of normal. The baseline was defined as the initiation of linezolid therapy.

Patient data including demographics, comorbidities, medication therapy, laboratory values, physiologic parameters, and indication for linezolid therapy were obtained from the electronic medical record. Renal function was assessed by serum creatinine and creatinine clearance (CrCL) estimated by the Cockcroft–Gault formula. Data organization and visualization were performed using R (version 3.6.0) and R Studio (version 1.2.1335).

The study was designed in accordance with legal requirements and the Declaration of Helsinki and was approved by the Ethical Committee of the First Affiliated Hospital of Wenzhou Medical University, China ([2021]034) and the Ethical Committee of Ruijin Hospital (KY2020-68). The study has been registered at the Chinese Clinical Trial Registry (ChiCTR2100047882). The informed consent was passed by the ethics committee in clinical research.

Routine clinical TDM data of patients treated with linezolid were retrospectively obtained from a database maintained at the Department of Pharmacy. The decision to administer linezolid and its dosing regimens (dose amount, dosing interval, duration of intravenous administration, and duration of therapy) was made by the attending physician. An opportunistic sampling strategy was performed when a steady-state concentration (attained after at least five continuous doses) of linezolid has been achieved. Dates and the exact time of linezolid treatment and TDM sampling were able to be indexed.

Plasma samples were separated by centrifugation for 5 min at 15,000 rpm immediately after sampling and stored at −80°C. Samples were determined within 24 h after sampling. The quantification of plasma concentration of linezolid was performed using a validated high-performance liquid chromatography–tandem mass spectrometry assay (Dai et al., 2021). The intra- and inter-day assay coefficients of variation were <10%, and the lower limit of quantification was 0.1 mg/L.

Toxicity was defined as follows: 1) thrombocytopenia: platelet count <125 × 10⁹ cells/L and a decrease of platelet count ≥25% in comparison with baseline levels; 2) anemia: a reduction of ≥25% of hemoglobin level compared to the baseline level; 3) leukopenia, white blood cell (WBC) count <4.0 × 10⁹ cells/L; and 4) hyperlactacidemia: serum pH < 7.35 and serum lactate >4 mmol/L (Dai et al., 2021). The baseline levels were defined at the initiation of linezolid therapy.

PopPK analysis was performed using non-linear mixed-effects modeling program NONMEM (version 7.4, Icon Development Solutions, Ellicott City, MD, United States) and Pirana (version 2.9.7). R (version 3.6.0) and Xpose (version 4.3.2) software packages were applied to generate diagnostic plots. The first-order conditional estimation method with inter- and intra-subject variability was used throughout the model development procedure.

One- and two-compartment structural models with first-order elimination were explored for the linezolid plasma concentration–time profiles. Between-subject variability (BSV) was modeled using exponential function. Residual variability was assessed using additive, proportional, and combined (additive plus proportional) error models. The base model was selected based on the visual inspection of diagnostic plots and various goodness-of-fit criteria, including precision and plausibility of parameter estimation, improvement of the objective function value (OFV), and the Akaike Information Criteria and Bayesian information criterion.

Exploratory analysis was performed before covariate modeling to examine the distribution of covariates in the population and the correlation between covariates of interest. Primary covariates included gender, age, height, body weight, white blood cell count, hemoglobin, platelet count, total bilirubin, albumin (ALB), alanine aminotransferase, aspartate aminotransferase, serum creatinine, and CrCL. Relationships between individual empirical Bayesian estimates of PK parameters and patient covariates were examined visually. Covariates were included using a stepwise forward selection process, with a threshold decrease in the OFV of 3.84 [p < 0.05, 1 degree of freedom (df)] until no further decrease in OFV was observed. All the significant covariates were then incorporated into the basic model to construct a full model. In backward elimination, the covariate was retained in the final model with a threshold increase in the OFV of 6.63 (p < 0.01, 1 df); otherwise, it was eliminated from the model. The additional criterion for retaining the covariate in the final model was a decrease in the unexplained BSV and an increase in PK parameter estimate precision.

Goodness-of-fit plots, non-parametric bootstrap, and visual predictive check were performed to evaluate the final model and parameter estimates. Goodness-of-fit plots include observed concentrations versus individual prediction, observed concentrations versus population prediction, conditional weighted residuals versus population prediction, and conditional weighted residuals versus time after the last dose. A non-parametric bootstrap procedure was conducted to assess the performance and stability of the final model. Random sampling with replacement was utilized to generate 1,000 replicate datasets using the individual as the sampling unit. The median and 95% confidence intervals of the resulting parameters were calculated and compared with the final parameter estimates obtained using NONMEM program. To evaluate the predictive performance, the statistics of the observed and simulated time–concentration profiles were compared using prediction- and variability-corrected visual predictive check. The dataset was simulated 1,000 times using the SIMULATION block in NONMEM for prediction- and variability-corrected visual predictive check. The 95% CI for the 5^th, 50^th, and 95^th percentiles of the simulated concentrations were calculated, plotted against time after the last dose, and compared with the observed concentrations.

Statistical analyses were performed using SPSS version 21.0 (IBM Corp). All study variables were summarized by descriptive statistics. The PK parameters such as the AUC_0-24 and C_min at day 1 (AUC_{0-24, d1} and C_{min, d1}), day 2 (AUC_{0-24, d2} and C_{min, d2}), and steady state (AUC_{0-24, ss} and C_{min, ss}) of each patient were predicted via the maximum a posteriori probability (MAP) Bayesian function of NONMEM using the final PK model as the Bayesian prior. The Kolmogorov–Smirnov test was performed to assess whether the data were normally or non-normally distributed. The relationship between the linezolid concentration–time profile and linezolid-associated toxicity was examined. Categorical variables were compared by using the chi-squared or Fisher exact test, and the continuous variables were compared using the Student t test or Mann–Whitney U test, as appropriate, comparing toxicity versus non-toxicity and dichotomous linezolid exposure thresholds.

Thresholds in the distribution of the linezolid exposure variables where the incidence of toxicity was most disproportionate were derived using the classification and regression tree (CART) analysis (Zhang, 1999). The predictive performance of linezolid exposure variables, including the CART-derived and other prior defined exposure thresholds, was evaluated using receiver operating characteristic (ROC) curves along with negative predictive value and positive predictive value (NPV and PPV, respectively). Multivariable logistic analyses were performed to determine the independent association between linezolid exposure variables and toxicity, while adjusting for confounding variables. Each linezolid exposure variable was assessed in a separate logistic regression model. All baseline variables associated with toxicity in the bivariable model analyses at a p < 0.1 were included in the explanatory multivariable model at model entry, and variables were excluded from the model using a backward stepwise approach. All statistical tests were two-sided. p values < 0.05 were considered statistically significant.

In total, 152 hospitalized adult patients with 270 plasma concentrations were included in the population PK analysis set; 141 patients were evaluable for hematological toxicity analysis. The demographic data and clinical baseline characteristics of these patients are summarized in Table 2, and the demographic profiles were similar across different analysis sets. The overall toxicity rate, including thrombocytopenia, anemia, leukopenia, and lactic acidosis, in this study was 43.26%. The median daily dose of linezolid selected by the physician was 18.83 mg/kg/d. The median [IQR] time of the sampling was 7.0 [4.0–10.0] days.

A total of 270 linezolid concentrations from 152 patients with a range of 0.25–24.79 mg/L were obtained for PopPK modeling. The linezolid concentration versus time after the last dose profile is shown in Figure 1.

The PK characteristics of linezolid could be well-illustrated by a 1-compartment model, with linear elimination showing the best fit of the observed concentration–time data based on the reduction in OFV and residual variability compared to the 2-compartment model. BSV was successfully estimated on both CL and V but not the absorption rate constant (Ka) and bioavailability (F) in the base model. The proportional error model was selected to evaluate the residual variability. Parameter estimates and diagnostic plots from the base model are provided in Supplementary Table S1, Figure S1.

The stepwise forward selection and backward elimination process are provided in Supplementary Tables S2, S3. The covariate model building identified CrCL and WBC as covariates of linezolid CL and ALB as a covariate of V. These covariates were not strongly correlated with one another in the population (|r|≤0.4). The final PopPK model is represented as follows: F = 80.6 % (1) K A ( h − 1 ) = 0.753 (2) C L ( L h ) = 2.93 + 2.33 × C r C L 73 + 0.685 × W B C 9.63 (3) V ( L ) = 97.5 − 1.42 × A L B (4) where CL is the individual clearance, V is the individual volume of distribution, CrCL is the estimated creatinine clearance, WBC is the white cell count in ×10⁹/L, and ALB is the serum albumin in g/L. The parameter estimates of the final model are displayed in Table 3.

Mean ± SD individual empirical Bayesian estimates of CL and V were 6.36 ± 2.50 L/h and 53.98 ± 14.59 L, respectively, across all patients with F estimated at 0.81 and KA estimated at 0.75 h⁻¹ in the population. Furthermore, patients included for toxicity analysis were divided into the group with linezolid-associated toxicity and the group without. The results of the Student t test indicated that the estimated CL values were significantly lower in patients with linezolid-associated toxicity than patients who did not experience toxicity (7.08 ± 2.65 L/h vs 5.86 ± 2.23 L/h). But no statistically significant difference in the estimated V value was observed between the two groups (55.88 ± 17.40 L vs. 52.63 ± 12.09 L).

The diagnostic goodness-of-fit plots of the final model are shown in Figure 2. The scatterplots of population prediction and individual prediction revealed an improvement in the final model compared with that of the basic model. The conditional weighted residuals vs the population prediction of the final model showed a stochastic distribution around zero, and most residuals were within an acceptable range (−2–2). The median with 95% CI parameter estimates obtained from a 1,000-run bootstrap analysis are given in Supplementary Table S4. The parameter estimates of the final PopPK model lie within the 95% CIs resulted from the non-parametric bootstrap procedure, and the biases between the final model estimates and bootstrapped median parameter estimates were < ±10% for all parameters, which demonstrated the good stability of the final model. The prediction- and variability-corrected visual predictive check of concentrations versus time after the last dose reflected a good fit between the observations and simulations (Supplementary Figure S2). Overall, the linezolid PopPK model evaluation results revealed that the final model provided an adequate description of the data and a good prediction of the individual PK parameters.

A total of 141 patients met the criteria for toxicity analyses. The incidence rate of thrombocytopenia, anemia, leukopenia, and hyperlactacidemia was 35.46, 19.15, 6.57, and 9.22% respectively. The Kaplan–Meier plot revealed that the median time from the initiation of therapy to the development of linezolid-induced toxicity was 12 days (Figure 3). Bivariate comparisons between patients who experienced and did not experience linezolid-related toxicity are listed in Table 4. Patients who experienced toxicity had significantly lower CrCL and longer treatment duration. In addition, patients with toxicity had significantly higher linezolid exposure quantified by AUC_{0-24, d1}, AUC_{0-24, d2}, AUC_0-24,ss, C_{min, d1}, and C_{min, d2}, or C_{min, ss}.

In the CART analysis, the rate of occurring toxicity was significantly higher among patients with AUC_{0-24, d1} ≥163 mg h/L, AUC_{0-24, d2} ≥207 mg h/L, AUC_{0-24, ss} ≥210 mg h/L, and C_min,d2 ≥6.9 mg/L, C_min,ss ≥6.9 mg/L, while no threshold was discovered for C_{min, d1}. The predictive performance of the CART-derived linezolid thresholds and other candidate thresholds are listed in Table 5. The CART-derived thresholds for AUC_{0-24, d1}, AUC_{0-24, d2}, AUC_{0-24, ss}, C_{min, d2}, and C_{min, ss} were the most predictive in receiver operating characteristic (ROC) curve analysis. Compared with other thresholds, the AUC_{0-24, d1} threshold of ≥163 mg h/L had lower predictive value. In addition, it was noted that the predictive performance of AUC_0-24 and C_min at day 2 was close to that at the steady state.

In the final logistic regression models, linezolid AUC_{0-24, d1} ≥ 163 mg h/L, AUC_{0-24, d2} ≥ 207 mg h/L, AUC_0-24,ss ≥ 210 mg h/L, C_min,d2 ≥ 6.9 mg/L, and C_min,ss ≥ 6.9 mg/L were still independently associated with linezolid-related toxicity (Table 6). In addition, the treatment duration of linezolid was also included as another independent risk factor in these final models. It was noted that CrCL was included in the final model when the exposure variable was AUC_{0-24, d1} ≥163 mg h/L.