Ethical approval for this study was obtained from the Hospital Ethics Committee of the Second Affiliated Hospital of Guizhou Medical University (No: 2023-LS-221). Patients were included if they were admitted to the surgical ICU and were receiving meropenem as a continuous infusion. Exclusion criteria included individuals under 18 years of age, as well as those undergoing renal replacement therapy (RRT) or extracorporeal membrane oxygenation (ECMO) during the course of meropenem administration. Meropenem dosing regimens were adjusted according to the measured creatinine clearance (CrCL). A loading dose of 2 g was administered over 30 min, immediately followed by continuous infusion of meropenem. Blood samples were collected from patients admitted to the surgical ICU between 1 March 2023, and 1 October 2024, after a minimum of 6 h of meropenem therapy. In addition to blood samples, data on biochemistry results (albumin, transaminase levels, renal function, etc.), demographic characteristics (age, weight, sex, etc.), SOFA scores, and APACHE II scores at admission were prospectively recorded at baseline and at the time of blood sampling. Detailed information regarding dosing (including date, time, dosage, and duration) and the timing of blood collection were also documented. Specially, for a 24-h continuous infusion, meropenem was therefore reconstituted and prepared three times per day.

The analysis of meropenem was conducted using liquid chromatography-mass spectrometry (LC-MS/MS). Reagents such as water, methanol (LC-MS grade), acetic acid, and ammonium acetate were purchased from Tedia Company Inc. (Beijing, China) and Sigma-Aldrich (Steinheim, Germany), while meropenem-D₆, an internal standard, was acquired from Toronto Research Chemicals (Ontario, CA). The mobile phase used in high-performance liquid chromatography (HPLC) comprised two components: phase A consisted of water containing 2 mM ammonium acetate and 0.1% acetic acid, while phase B consisted of methanol. Sample preparation for plasma and dialysate was performed according to the method described by Koal et al. (2006). To each sample, 10 µL of methanolic meropenem-D6 (target concentration: 1 mg/mL) was added as an internal standard.

Chromatographic separation was achieved using an Acquity BEH C18 column (2.1 mm × 50 mm, 1.7 μm particle size, Waters, United States), coupled with a guard cartridge. The separation gradient was implemented as follows: initially, the mobile phase consisted of 95% phase A and 5% phase B for 0 min. This was gradually changed to 10% phase A and 90% phase B at 4 min. Subsequently, a second gradient was applied, transitioning from 10% phase A and 90% phase B at 4 min to 95% phase A and 5% phase B at 6 min. The system was then equilibrated for 3 min before the next injection. Mass spectrometric analysis was performed using a TSQ Quantum Ultra triple-quadrupole mass spectrometer (Thermo Fisher Scientific Inc., Boston, United States) operating in positive electrospray ionization (ESI) mode. The parameters included a spray voltage of 4.0 kV, capillary temperature of 350 °C, vaporizer temperature of 280 °C, sheath gas pressure of 20 psi, auxiliary gas pressure of 5 psi, collision gas pressure of 2.0 mTorr. The transitions of the precursors to the product ions were m/z 384.1 → m/z 141.2 and m/z 254.2 for meropenem, and m/z 390.1 → m/z 147.2 for meropenem-D₆. The limit of quantification (LOQ) was 0.1 mg/L for plasma, with the assay demonstrating linearity from 0.1 to 100 mg/L. The validation of this bioassay was evaluated, and meet requirements as per guidelines for analysis of biological samples.

PK modeling was performed using Phoenix NLMN™ software (version 7.0, Certara L.P Pharsight, St. Louis, MO, United States), employing a systematic three-step approach: (i) selection of the basic model, (ii) identification of relevant covariates, and (iii) validation of the model. The first-order conditional estimation method with interaction was applied for parameter estimation. Interindividual variability (IIV) was modeled using a log-normal distribution, while residual variance was assessed through additive, proportional, and combined error models. The performance of the model was evaluated using several criteria: (i) the reduction in the objective function value (−2 log likelihood), (ii) graphical comparison of predicted versus observed concentrations, and (iii) examination of changes in the standard error of parameter estimates. These assessments ensured the reliability and robustness of the final model (Li et al., 2022).

In the second step, potential covariates, including both demographic and clinical variables, were evaluated for their inclusion in the basic model. The influence of these covariates on model parameters was assessed using graphical methods and generalized additive models. Continuous covariates were examined as either proportional or power functions, whereas categorical variables were incorporated using an appropriate mathematical equation P_j = PPOP_COV · (1–Covi), where P_j is the PK parameter for patient j, Covi is the covariate, and PPOP is the typical value of the PK parameter. Covariates that improved model fit were retained, based on biological plausibility, graphical concordance, and statistical tests.

To assess the stability and performance of the final model, a nonparametric bootstrap resampling technique was employed (Mathew et al., 2016). The original dataset was resampled at the subject level to produce 1,000 new replications. These resampled datasets were then used to calculate the 2.5th and 97.5th percentiles of the model parameters. If the model was valid, the parameter estimates from the original dataset should closely align with the median and fall within the 2.5th and 97.5th percentiles of the distributions generated through resampling (Li et al., 2022).

The performance of the final model was further assessed using a prediction-corrected visual predictive check (pc-VPC), a technique that integrates virtual predictions and observed data derived from Monte Carlo simulations. In this approach, the percentiles of the simulated data were compared with the corresponding percentiles of the observed data, thereby providing a thorough evaluation of the model’s predictive accuracy.

Monte Carlo simulations were conducted to simulate concentration-time profiles for 1,000 patients per dosing regimen, using the final population PK parameters. Different loading dosages (LD) with four bolus maintenance dosing (MD) regimens were simulated as follows: a. 2000 mg LD+ 1,000 mg every 8 h; 2000 mg LD+ 1,500 mg every 8 h; 2,500 mg LD+ 2000 mg every 8 h; 3,000 mg LD+ 2,500 mg every 8 h b. 750 mg LD+ 500 mg every 6 h; 1,000 mg LD+ 750 mg every 6 h; 1,500 mg LD+ 1,000 mg every 6 h; 1,500 mg LD+ 1,500 mg every 6 h.

The simulations incorporated patient age groups of 20, 40, 60, and 90 years. For each regimen, the proportion of patients achieving the PK/PD target of 40% fT > MIC × 4 was calculated. MIC distributions ranging from 0.008 to ≥64 mg/L were considered as the pharmacodynamic (PD) parameter for meropenem, following the Clinical and Laboratory Standards Institute (CLSI) Executive Standards for Antimicrobial Susceptibility Testing (Humphries et al., 2021). To ensure the convenience of clinical application, MIC valuesof 0.25, 0.5, 1, 2, 4, 8, and 16 mg/L were selected based on a study published previously (Hou et al., 2024). PK/PD target values were set as 40%fT > MIC × 4 and 100%fT > MIC × 4 (Zhang et al., 2024; Roberts et al., 2025), to evaluate both the therapeutic efficacy and the potential for resistance suppression of meropenem. These target values were selected on the basis of previously published criteria for carbapenem efficacy and resistance prevention (Zhang et al., 2024). A clinical susceptibility breakpoint was set to 2 mg/L, consistent with meropenem breakpoints reported by McKinnon et al. (2008). For each dosing regimen, the PTA was determined and subsequently utilized to calculate the CFR against the bacterial population, following approaches described in earlier studies (Koomanachai et al., 2016).

Statistical analysis was performed using SPSS software (version 20 for Macintosh, IBM SPSS Statistics, United States). Categorical variables were expressed as both absolute and relative frequencies, whereas continuous variables were reported as medians with their respective ranges. For comparisons involving normally distributed variables, a two-tailed Student’s t-test was applied, while the Mann-Whitney U test was used for non-normally distributed data. The chi-square test or Fisher’s exact test, depending on the context, was employed for categorical variables. A significance level of p ≤ 0.05 was considered statistically significant for all analyses.

Data from 144 adults and elderly patients (mean (SD) age 60.63 (18.37) years, 58 males and 86 females) in Intensive Care Unit were initially recruited between 1 March 2023 and 1 October 2024, for this study. A summary of the demographic characteristics of the participants is provided in Table 1. The median (IQR) value of trough concentration was 10.40 (3.39,22.53) mg/L. Patients received 2000 mg of meropenem dosing as an initial the loading dose, followed by daily maintenance dose of 500 or 1,000 mg (q8h/q12h).

A total of 144 blood samples were obtained from the enrolled patients. Among these, 5 samples (3.47%) had concentrations below the LOQ. For those observations that were considered normally distributed, values below the LOQ were substituted with random values drawn between negative infinity and the lower LOQ, referred to M3 method, as described previously (Palaparthy et al., 2019; Li et al., 2022). Consequently, all detectable concentrations, including those below the LOQ, were treated as continuous data for population pharmacokinetic (popPK) analysis.

The pharmacokinetic profiles of meropenem was modeled using a one-compartment model with first-order elimination. This model demonstrated a significantly lower objective function value (OFV) of 873.81, compared to the two-compartment model, which yielded an OFV of 928.18. The model described the pharmacokinetics of meropenem in critically ill patients, with the corresponding parameter estimates provided in Table 1. Residual variability was adequately accounted for using an additive error model, as the combined proportional and additive error model did not offer a significant improvement. Furthermore, IIV in CL and Vd decreased from 8.91% to 4.22% in the base model to 4.50% and 1.60% in the final model, respectively. Residual variability showed a minor increase from 0.72 in the base model to 0.74 in the final model. These findings highlight the refinement and improved fit of the model with respect to both individual and residual variabilities.

During the covariate analysis, potential clinical and demographic factors, including age, body weight, sex, serum creatinine, and renal function, were systematically evaluated for their impact on meropenem clearance and volume of distribution. In the forward selection phase, age was initially included in the base model for both Cl and Vd, leading to a significant reduction of 10.36 points in the objective function value (OFV). The inclusion of estimated glucose (GLU) also resulted in a modest, but statistically significant, decrease of 2.16 points in the OFV. However, during the backward elimination step, estimated glomerular filtration rate (eGFR) was excluded from the model. And ultimately, age was the only covariate retained which had the most substantial impact on the meropenem Vd, contributing the largest decrease in OFV among all potential covariates. Compared with the base model, the final PK model led to a significant decrease of Akaike Information Criterion (−1.02%) and Bayesian Information Criterion (−1.33%). Estimates of the shrinkage for random variabilities of Vd and CL were 0.35 and 0.12, respectively. The condition number was 10.93. Model diagnostics, as illustrated in Figure 1, indicated that the final model for meropenem pharmacokinetics demonstrated an acceptable goodness-of-fit. In this final model, the individual pharmacokinetic parameters were optimally described using the following equations: C L L / h = 4.68 × e 0.045 V d L = 4.47 × A g e 63.5 0.29 × e 0.016

The median (range) values of the estimated age-normalized Vd at steady state were 4.41 (4.21–4.62) L. As illustrated in Figure 2, the age-normalized Vd increased with advancing age, indicating a clear trend across different age groups.

The robustness and stability of the final pharmacokinetic model were evaluated using both bootstrapping (Table 2) and prediction-corrected visual predictive checks (Figure 3). In the bootstrapping analysis, the median of the fixed-effects parameter estimates obtained from the resampled datasets was within 5% of the population parameter estimates derived from the original dataset for all parameters, suggesting high reliability and accuracy. Moreover, the comparison between predicted and observed meropenem concentrations demonstrated excellent concordance. Specifically, the 5th, 50th, and 95th percentiles of the predicted concentrations were closely aligned with those observed in the clinical data. This alignment further supports the model’s predictive accuracy. Model diagnostics confirmed that the final meropenem pharmacokinetic model exhibited a satisfactory goodness of fit. Figure 1 shows that there were no discernible trends between the population predicted concentrations (PRED) and the conditional weighted residuals (CWRES), nor was there any correlation between time and CWRES. Additionally, as presented in Table 2, the median values of the parameter estimates obtained through bootstrapping were consistent with those from the final model. This consistency further indicates the stability of the model and reinforces the reliability of the parameter estimates derived from the population pharmacokinetic analysis. Finally, the results from the pcVPC, shown in Figure 3, revealed that the 5th, 50th, and 95th percentiles of the observed concentrations fell within the 95% confidence interval of the predicted concentrations, thereby demonstrating the model’s predictive accuracy and overall adaptability to the data.

Monte Carlo simulations (Figure 4) demonstrated that, in critically ill patients infected with P. aeruginosa, all meropenem dosing regimens achieved a PTA exceeding 90% when the MIC was 0.25 mg/L in both every-8-h (Figures 4a–d) and every-6-h (Figures 4e–h) dosing regimen. At an MIC of 0.5 mg/L, all meropenem dosing regimens achieved a PTA >90% at each age group, while escalation of the meropenem MD dose to 1,500 mg every 8 h was required to maintain adequate exposure in the 20-year-old group (Figure 4a). For an MIC of 1 mg/L, a regimen of LD 3000 mg and MD 2500 mg every 8 h achieved a PTA of 92.20%, 97.20%, 98.20% and 99.40% for the 20-year-old group, the 40-year-old group, the 60-year-old group and the 90-year-old group, respectively. When the MIC increased to 2 mg/L, shortening the dosing interval to 6 h was necessary, yielding PTAs exceeding 90% in each age group. Using an MIC of 2 mg/L as the clinical susceptibility breakpoint, the regimen of 1,500 mg LD followed by 1,000 mg q6h was the most appropriate option for patients aged 20 years, and 1,000 mg LD followed by 750 mg q6h for 40–90 years, as it consistently achieved PTAs above 90%. At MIC values of 4 mg/L, only the regimen of 1,500 mg LD followed by 1,500 mg q6h attained a PTA greater than 90% in patients aged ≥60 years old. At MIC values ≥ 8 mg/L, however, none of the simulated regimens attained a PTA greater than 90%.

When the more stringent pharmacodynamic threshold of 100% fT > MIC × 4 was applied to assess resistance suppression (Figure 5), simulations indicated that administering meropenem at 1,500 mg every 6 h in patients aged over 60 years achieved a PTA of 50.00% at an MIC of 0.25 mg/L. For MICs ≥0.5 mg/L, no regimen achieved PTA values above 50%; among the tested regimens, LD 1500 mg and MD 1500 mg every 6 h provided the highest relative PTA, though it remained below the desired threshold.

When the PK/PD target was defined as 40% fT > MIC×4, the CFR for the regimens consisting of LD of 1,500 mg followed by a MD of either 1,000 mg q6h or 1,500 mg q6h exceeded 90% across the population distribution of MICs (Hou et al., 2024) suggesting that both dosing strategies should be prioritized as empirical treatment options for P. aeruginosa in 20-year-old patients (Table 3). The regimen of LD 2500 mg + MD 2000 mg q8h, LD 3000 mg + MD 2500 mg q8h, LD 1000 mg + MD 750 mg q6h, LD 1500 mg + MD 1000 mg q6h, and LD 1500 mg + MD 1500 mg q6h achieved CFRs above 90% in 40-year-old and 60-year-old patients. With the exception of the regimen of LD 2000 mg + MD 2000 mg q8h, all other dosing strategies achieved CFR values exceeding 90% in 90-year-old patients. When considering a more stringent PK/PD target,100% fT > MIC × 4, only the regimen of LD 1500 mg + MD 1500 mg q6h achieved the highest relative CFRs of 12.59%, 24.31%, 32.06% and 40.90% for 20-year-old, 40-year-old,60-year-old and 90-year-old patients, respectively (Table 3).

According to previous studies, achieving a target of at least 40%–50% fT > MIC×4 during the dosing interval is generally considered sufficient to ensure adequate therapeutic efficacy in critically ill patients (Ariano et al., 2005; McKinnon et al., 2008). Based on this criterion and the MIC distribution of P. aeruginosa, the recommended empirical regimens are LD 1500 mg followed by MD 1000 mg every 6 h for 20-year-old patients, LD 1000 mg followed by MD 750 mg every 6 h for 40–60-year-old patients, and LD 750 mg followed by MD 500 mg every 6 h for 90-year-old patients, thereby safeguarding the antimicrobial efficacy of meropenem across different age groups.