This was a retrospective, observational study of hospitalized critically ill 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. The inclusion criteria were as follows: 1) patients receiving intensive care aged ≥18 years with confirmed or suspected Gram-positive bacterial infections; 2) patients who received intravenous linezolid for at least 3 days; and 3) patients for whom 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; 2) patients who died within 24 h after being treated with linezolid; 3) patients with baseline PLT <75 × 10⁹ cells//L; 4) patients with baseline hemoglobin <6.8 g/dl for males or <6 g/dl for females; 5) patients with baseline absolute neutrophil count <500 cells/μl; and 6) patients with baseline TBIL > 5 times the upper limit of normal. The baseline was defined at the initiation of linezolid therapy.

The study was designed in accordance with legal requirements and the Declaration of Helsinki and was approved by the Ethical Committees of the First Affiliated Hospital of Wenzhou Medical University, China [(2021)034], registered at the Chinese Clinical Trial Registry (ChiCTR2100047882). The informed consent was free passed by the Ethics Committee in Clinical Research (ECCR) of the First Affiliated Hospital of Wenzhou Medical University.

Routine therapeutic drug monitoring (TDM) data of 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) were made by the attending physician. An opportunistic sampling strategy was performed when a steady-state concentration of linezolid has been achieved (at least 3 days from the start of treatment). Most of the TDM samples collected were peak or trough concentrations. The concentration–time profile of linezolid for the included patients is shown in Supplementary Figure S1. 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 and were determined within 24 h after sampling. Quantification of the plasma concentration of linezolid was performed using a validated high-performance liquid chromatography–tandem mass spectrometry (LC-MS/MS) assay (Phillips et al., 2001). Intra- and inter-day assay coefficients of variation were <10%, and the lower limit of quantification was 0.1 mg/L. Individual laboratory parameters and demographic data of patients, including gender, age, height, body weight (WT), white blood cell count (WBC), hemoglobin (Hb), platelet count (PLT), total bilirubin (TBIL), serum albumin (ALB), alanine aminotransferase (ALT), aspartate aminotransferase (AST), and serum creatinine (SCr), were collected. CrCL was calculated using the Cockcroft–Gault equation. Data organization and visualization were performed using R (version 3.6.1) and R Studio (version 1.3.1093).

Myelosuppression was defined as follows: 1) thrombocytopenia: PLT <125 × 10⁹ cells/L and a decrease of ≥25% PLT in comparison with the baseline levels; 2) anemia: a reduction of ≥25% of Hb compared with the baseline. The baseline levels were defined as the hematological parameters at the initiation of linezolid therapy.

Population PK analysis was performed using the nonlinear 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), Akaike information criteria (AIC), and Bayesian information criterion (BIC).

Covariates were included using a stepwise forward selection process with a threshold decrease in OFV of 3.84 (p < 0.05, 1 degree of freedom [df]) until no further decrease in OFV was observed. All of 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 OFV of 6.63 (p < 0.01, 1 df); otherwise, it was eliminated from the model. The additional criteria for retaining the covariate in the final model were a decrease in the unexplained BSV and an increase in PK parameter estimate precision.

For internal validation of the final model, goodness-of-fit plots, including observed concentrations (DV) versus individual prediction (IPRED), DV versus population prediction (PRED), conditional weighted residuals (CWRES) versus PRED, and CWRES versus time after last dose (TAD), were performed. A nonparametric bootstrap procedure was conducted to assess the performance and stability of the final model (Ette and Onyiah, 2002; Ette et al., 2003). 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 (CIs) of the resulting parameters were calculated and compared with the final parameter estimates obtained using the 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 (pvcVPC) (Bergstrand et al., 2011). The dataset was simulated 1,000 times using the $SIMULATION block in NONMEM^® for pvcVPC. The 90% CIs for the 5th, 50th, and 95th percentiles of the simulated concentrations were calculated, plotted against time after the last dose, and compared with the observed concentrations.

Monte Carlo simulations were performed using the final model to identify the pragmatic dose adjustment of linezolid. Simulations were performed with the covariate sets from individual patients included in model building, serving as a template for 1,000 simulated subjects. The probability of target attainment (PTA) of the C_min target was calculated for each dose (daily doses of 300–1,800 mg/day with increments between 150 and 300 mg) among simulated subjects with CrCL values of <30 ml/min, 30–59 ml/min, 60–89 ml/min, and ≥90 ml/min.

Statistical analyses were performed using SPSS version 21.0 (IBM Corp). All study variables were summarized by descriptive statistics. The C_min and AUC_0-24 at steady state 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. AUC_0-24 was calculated by the linear-log trapezoidal rule using the concentrations at continuous time (10 min-interval) predicted via Bayesian estimation. Specifically, the linear trapezoidal approach was used during the ascending phase and the log-linear method was used during the descending phase. Logistic regression modeling for the toxicity of linezolid was conducted to determine whether the trough concentration or AUC_0-24 of linezolid was a significant predictor of myelosuppression during treatment. The time from the initiation of linezolid treatment to the development of myelosuppression was estimated using the Kaplan–Meier curve analysis. All statistical tests were two-sided. p values <0.05 were considered statistically significant.

In total, 83 critically ill patients with a mean ± S.D. age of 60.57 ± 14.64 years and body weight of 64.11 ± 11.21 kg were included in the present study. The demographic data of the included patients are summarized in Table 1. The overall rate of myelosuppression occurrence, including thrombocytopenia and anemia, in this study was 40.97%. The incidence rates of thrombocytopenia and anemia were 36.14 and 16.87%, respectively. The median daily dose of linezolid selected by the physician was 20.17 mg/kg/d, and most patients received the drug every 12 h.

A total of 127 linezolid concentrations derived from 83 patients at a range of 0.25–21.61 mg/L were obtained for population PK modeling. The PK characteristics of linezolid illustrated by a 1-compartment model with linear elimination showed the best fit of the observed concentration–time data based on reduction in OFV and residual variability. The intra-individual variability for V was fixed as 0 because it was very small (ŋ<0.0025); it might be because of the limited data of peak concentrations. The proportional error model was used to evaluate the residual variability. Parameter estimates and diagnostic plots of the base model are provided in Supplementary Table S1 and Supplementary Figure S2. Covariate model building identified CrCL as the only covariate on linezolid CL that significantly explained interindividual variability. The final PK model is represented as follows: C L ( L / h ) = 3.66 + 2.18 × C r C L 65 , (1) V ( L ) = 54 , (2) where CL is the individual clearance, V is the individual volume of distribution, and CrCL is the estimated creatinine clearance. The parameter estimates of the final model are displayed in Table 2.

The mean ± SD individual empirical Bayesian estimate of CL was 6.35 ± 2.32 L/h, across all patients with V estimated at 54 L. Furthermore, dividing the included patients into the group with linezolid-induced myelosuppression and the group without, the estimated CL was significantly lower in patients with linezolid-induced myelosuppression compared with that in the patients without (4.91 ± 1.51 L/h vs. 7.57 ± 2.35 L/h, p = 0.027).

The diagnostic goodness-of-fit plots of the final model are shown in Figure 1. The CWRES versus PRED of the final model showed a stochastic distribution around zero, and most residuals were within an acceptable range (−2 to 2). The median with 95% CI parameter estimates obtained from a 1,000-run bootstrap analysis is given in Supplementary Table S2. The parameter estimates of the final PK model lay within the 95% CIs from the nonparametric bootstrap procedure, and the biases between the final model estimates and the bootstrapped median parameter estimates were within ±10% for all parameters, which demonstrated the good stability of the final model. The pvcVPC of concentrations versus time after the last dose reflected a good fit between the observations and simulations (Supplementary Figure S3). Overall, the linezolid population PK model evaluation results revealed that the final model provided an adequate description of the data and a good prediction of individual PK parameters.

A strong correlation (r2 = 0.9969) was found between the C_min and AUC0-24. The linear regression line was AUC0-24 = 26.354×C_min+91.607. The ratio of the AUC to MIC was determined to be the most important PK/PD index for linezolid against strains of S. aureus and strains of S. pneumoniae with varying antibacterial susceptibilities, and the antibacterial activity of linezolid may be maximized when AUC0-24/MIC ratio ≥80, which was derived from a study conducted in critically ill patients (Rayner et al., 2003). Given the MIC90 of linezolid for most staphylococci is 2 mg/L (Jones et al., 2007), to achieve the PK/PD index (AUC0-24/MIC ≥80) for the efficacy of linezolid, an AUC0-24 of 160 mg h/L was required when the MIC was 2 mg/L. Using the linear regression equation, the C_min for AUC0-24 = 160 mg h/L was 2.6 mg/L. Thus, the C_min required for sufficient efficacy was ≥2.6 mg/L.

Linezolid-induced myelosuppression was observed in 34 (40.97%) of the 83 patients. Figure 2 shows the relationship between the C_min and the myelosuppression (absence, 0; presence, 1). The C_min of linezolid was a significant predictor of myelosuppression in critically patients according to the following equation: probability of myelosuppression = 1/[1 + exp (3.767–0.481*C_min)]. The threshold for the C_min of linezolid that caused myelosuppression with 50% probability was 7.8 mg/L. In general, the target C_min range based on the aforementioned results was considered to be 2.6–7.8 mg/L. In addition, the Kaplan–Meier plot revealed that the overall median time from the initiation of therapy to the development of myelosuppression was 12 days (Figure 3).

The simulated PTA of linezolid therapeutic range (C_min 2.6–7.8 mg/L) with various dosage regimens in patients with various CrCLs were quantified in Table 3 and Figure 4. The target for PTA was defined as > 80%. The low probability of attaining therapeutic C_min values in patients with renal impairment (CrCL <60 ml/min) with a standard linezolid dose of 600 mg every 12 h is primarily due to an increase in the probability of attaining supratherapeutic C_min (>7.8 mg/L). An empirical dose reduction to 600 mg every 24 h was the optimal to balance safety and efficacy in patients with CrCL of 30–59 ml/min, whereas 450 mg every 24 h was the alternative for patients with CrCL <30 ml/min.