This was a prospective study conducted at the Department of Pulmonology and Clinical Care Medicine, Haeundae Paik Hospital, Busan, Republic of Korea, from November 2018 to November 2020. Thirty patients (aged ≥19) who underwent ECMO for respiratory and/or cardiac dysfunction and who received meropenem for treatment or prophylaxis were included in the clinical study. A written informed consent form was obtained from and signed by the legally authorized representative of each subject before enrollment. This study was approved by the Institutional Review Board of Inje University Haeundae Paik Hospital (IRB No. 2018-06-017) and conducted in accordance with the Declaration of Helsinki and Good Clinical Practice. The baseline demographic factors between continuous renal replacement therapy (CRRT) group and non-CRRT groups were compared. If the parameters for both the groups are normally distributed, the t-test was used; if either of the two groups did not satisfy the normality, the Wilcoxon rank-sum test was used.

The ECMO system was the Permanent Life Support System (MAQUET, Rastatt, Germany) consisting of a PLS-i oxygenator and a ROTAFLOW Centrifugal Pump. The circuit was primed with 1 L of normal saline or plasma solution. The total circuit volume was 500–600 ml.

A single 500 or 1,000 mg dose of meropenem diluted in 200 ml of 5% dextrose in water was infused intravenously over 3 h to patients on ECMO. After the first dose. 7 blood samples were drawn from each patient’s arterial catheter into heparinized tubes. The predetermined three sampling schemes were as follows: Scheme 1 = 0 (predose), 3.33, 3.67, 4, 5, 6, and 8 h; Scheme 2 = 0, 3.33, 3.67, 4, 5, 8, and 10 h; or Scheme 3 = 0, 3.33, 3.67, 4, 6, 11, and 14 h after the start of meropenem administration. The blood sampling times were determined by considering the blood sampling times and distribution half-life (0.498 and 0.504 h) and elimination half-life (1.67 and 2.09 h) of two previous studies (Doh et al., 2010; Lu et al., 2016). The plasma samples were separated via centrifugation (2000g at 4°C for 10 min) within 30 min of sampling. They were transferred to polypropylene tubes by 1 ml and stored at −70°C for 1–6 months until assayed.

Plasma meropenem concentrations were analyzed using a validated high-performance liquid chromatography (HPLC)–tandem mass spectrometry assay. The HPLC was performed on an Agilent 1,200 series with an Atlantis C₁₈ column (Company, Waters, Milford, MA, USA) (2.0 mm × 150 mm, 3.0 μ). Mass spectrometric detection was performed using a triple-quadrupole mass spectrometer (SCIEX API4000, Applied Biosystems, Foster City, CA) with an electrospray ionization interface. Data acquisition and processing were accomplished using Analyst software (version 1.4.1; Applied Biosystems, Foster City, CA). The lower limit of quantitation was 1 mg/L. The assay results were linear over a range of 1–50 mg/L (R² = 0.9974). Inter-day precision and accuracy of the validation concentration range (1, 2, 5, 10, 25, and 50 mg/L) analyzed with standard samples for 3 days were 0.5–2.7% and 89.9–100.0%, respectively.

Population PK modeling was implemented using NONMEM 7.5 (Icon Development Solutions). A first-order conditional estimation with interaction method was used during analysis to account for potential interactions involving between-subject variability (BSV) for PK parameters and residual variability (RV), caused by assay error, model misspecification, errors in independent variables, and intra-individual variability, etc. One-, two-, and three-compartment models were tested using ADVAN1 TRANS2, ADVAN3 TRANS4, and ADVAN11 TRANS4 from the pharmacokinetic model library in NONMEM. First-order kinetics was assumed for all PK processes other than the zero-order infusion.

The PK parameters were assumed to follow a log-normal distribution. The parameter model was defined as θ_i = θ × exp(η_i), where θ is the median value of the PK parameter, θ_i is an individual parameter, and η_i is a random effect that is assumed to be normally distributed with a mean of 0 and variance of ω². Additive, proportional, and combined additive and proportional error models were tested for RV, which was assumed to be normally distributed with a mean of 0 and variance of σ² (Dosne et al., 2016).

Model evaluation and selection were based on objective function values (OFVs) by NONMEM, relative standard errors for parameter estimates, shrinkage of BSV, and diagnostic goodness-of-fit plots. In a log-likelihood ratio test, an OFV reduction (ΔOFV) greater than 3.84 between two nested models with one degree of freedom or greater than 5.99 with two degrees of freedom was considered a significant model improvement. Diagnostic plots included conditional weighted residuals (CWRES) v. time, CWRES v. model-predicted population concentration (PRED), observation v. PRED, and observation v. model-predicted individual concentration (IPRED) (Hooker et al., 2007).

Perl-speaks-NONMEM software (version 5.0.0, available at https://uupharmacometrics.github.io/PsN/) was used to search for significant covariates and evaluate the final model using a visual predictive check and nonparametric bootstrap method. Stepwise forward selection and backward elimination were conducted to identify significant covariates for PK parameters. The statistical significance criteria were p < 0.01 (ΔOFV < -6.63 with 1 degree of freedom) for selection and p < 0.001 (ΔOFV >10.8 with 1 degree of freedom) for elimination. A significant covariate was required to have clinical relevance and to meet the statistical criteria. The tested covariates for all PK parameters were sex, age, weight, height, body surface area (BSA), serum protein level, serum albumin level, primary diagnosis, comorbidity, Acute Physiology and Chronic Health Evaluation II (APACHE II) score, Sequential Organ Failure Assessment Score (SOFA) score, presence of CRRT, ECMO flow rate, and ECMO type (veno-arterial or veno-venous). Serum creatinine level and renal function as estimated by applying the Cockcroft-Gault (CG), Chronic Kidney Disease Epidemiology Collaboration (CKD-EPI), modified CKD-EPI, Modification of Diet in Renal Disease (MDRD), and modified MDRD equations were tested only for total clearance. The modified CKD-EPI and MDRD values were calculated using individual body surface area (BSA) values, where BSA was calculated by applying the Du Bois formula. The final PK parameter estimates between continuous renal replacement therapy (CRRT) group and non-CRRT groups were compared. If the values for both the groups are normally distributed, the t-test was used; if either of the two groups did not satisfy the normality, the Wilcoxon rank-sum test was used.

Prediction- and variability-corrected visual predictive checks were performed by comparing the observed plasma concentrations with 80% prediction intervals from 1,000 simulated datasets applying the final PK model. Virtual observations were compared with prediction- and variability-corrected observed concentrations of meropenem to evaluate the model performance (Bergstrand et al., 2011). The nonparametric bootstrap method was used to evaluate the stability of the final model. The median and 95% confidence intervals for the parameter estimates of bootstrap samples (n = 2000) were generated to compare with the final PK parameter estimates. The R statistical software package (version 4.0.3, available at https://www.r-project.org/) was used for visualization and postprocessing of modeling output.

Two Monte Carlo simulations were conducted. The first simulation was implemented to evaluate the appropriateness of the recommended dosing regimen (for a creatinine clearance [CL_CR] > 50 ml/min, 1 g q8h; for a CL_CR of 26–50 ml/min, 1 g q12 h; for a CL_CR of 10–25 ml/min, 500 mg q12 h; for a CL_CR < 10 ml/min, 500 mg q24 h) when treating adult patients infected with intra-abdominal infection by Pseudomonas aeruginosa (P. aeruginosa). A total of 10,000 individual PK parameters of virtual patients was generated assuming a log-normal distribution for each PK parameter or each covariate. Then, ten thousand MICs were randomly assigned to the virtual patients. The clinical breakpoint distribution of MICs set by the European Committee on Antimicrobial Susceptibility Testing (EUCAST) was used to simulate the MICs. The steady-state concentration-time profiles (in minutes) for the virtual patients were generated to investigate the probability of target attainment (PTA). The target index for meropenem is the cumulative percentage of a 24 h period during which the free (f) drug concentration exceeds the minimum inhibitory concentration (MIC) at steady-state conditions (fT_>MIC) (Drusano, 2004). The tested targets were 40% fT_>MIC, 100% fT_>MIC, and 100% fT_>4xMIC for meropenem. A dosing regimen was considered optimal if the PTA is equal to or greater than 90%. The f was fixed at 0.98 (Ulldemolins et al., 2011).

The second simulation was implemented to investigate the optimal dosing regimen for the three targets. Thousand individual PK parameters were generated assuming a log-normal distribution for each PK parameter. The one of the finally selected two covariates, CL_CR, was generated by applying a uniform distribution within the range 0–170 ml/min and the virtual patients were assigned to the six renal function groups (0–10, 10–25, 25–50, 50–90, 90–130, or 130–170 ml/min). The other covariate, ECMO flow rate, was fixed to the median value of 3.7 L/min. Then, the PTA for the generated steady-state concentration-time profiles were evaluated with various combinations of the three doses (0.5, 1, and 2 g), two infusion times (0.5 and 3 h), two dosing intervals (8 and 12 h), and MICs (0.060, 0.125, 0.25, 0.5, 1, 2, 4, 8, and 16 mg/L). The thousand parameters were also used to investigate the effect of ECMO flow rate on the PTA of 40% fT_>MIC. The steady-state concentration-time profiles were generated at the flow rate of 2, 4, and 6 L/min and evaluated with various combinations of the three doses (0.5, 1, and 2 g), two dosing intervals (8 and 12 h), and MICs (0.060, 0.125, 0.25, 0.5, 1, 2, 4, 8, and 16 mg/L), while the infusion time was fixed to 0.5 h.

Prospectively, 30 patients were enrolled in this study, of whom 10 patients received CRRT (Table 1). The most common primary diagnosis was pneumonia in both the non-CRRT group (n = 11) and the CRRT group (n = 5). Of the types of ECMO employed, venovenous (VV)-ECMO was used in 75% (n = 15) and venoarterial (VA)-ECMO in 25% (n = 5) in non-CRRT group and VV-ECMO used in 80% (n = 8) and VA-ECMO used in 20% (n = 2) in the CRRT group. Serum creatinine level and eGFR by MDRD, modified MDRD, CKD-EPI, and modified CKD-EPI of the patients with CRRT were significantly different from those of the patients without CRRT, while CLCR by Cockcroft-Gault formula was not. The median (IQR) ECMO flow rates of the non-CRRT and CRRT group were 3.62 (2.62-4.08) L/min and 3.83 (3.47-4.17), respectively, and were not significantly different (p = 0.6511). All patients received 1,000 mg of meropenem with the exception of one in the non-CRRT group, who received 500 mg.

The median (IQR) durations of ECMO therapy of the non-CRRT and CRRT groups were 15.5 (7.40–20.3) days and 9.00 (6.00–13.8) days, respectively. Ten patients (50%) in the non-CRRT group and five patients (50%) in the CRRT group died in the intensive care unit (ICU). The mortality rates for the patients in the non-CRRT group receiving VV-ECMO or VA-ECMO were 40% (n = 6) and 80% (n = 4), respectively. Those for CRRT group were 50% (n = 4) and 50% (n = 1), respectively. The most common causes of death were multi-organ failure and sepsis.

A total of 210 plasma samples was used to build the population PK model. The time course of meropenem concentrations was well described by a two-compartment model. The objective function values (OFVs) for one-, two-, and three-compartment models without covariates were 547.667, 409.786, and 395.043, respectively. However, the three-compartment model failed to achieve model convergence and generated poor parameter estimates. The structural PK parameters for the two-compartment model were total clearance (CL), central volume of distribution (V_C), volume of distribution for the peripheral compartment (V_P), and intercompartmental clearance between V_C and V_P (Q), as shown in Table 2. All PK parameter estimates were not significantly different between the CRRT and non-CRRT groups (Table 3). Individual model fits are shown in Supplementary Figure S1.

In the final PK model (OFV = 369.100), CL_CR as estimated by the Cockcroft-Gault formula was identified as a significant covariate for CL, while the OFV of a reduced model without this covariate increased to 391.854. The random BSV for CL was reduced from 57.9 to 39.6% after including the covariate. The ECMO flow rate was a significant covariate for V_C, while the OFV of a reduced model without ECMO flow rate on V_C increased to 387.207. The random BSV for V_C was reduced from 67.1 to 48.5% after including the covariate. The RV was well described by a combined additive and proportional error model. Model robustness was supported by the bootstrap median values and the 95% confidence intervals for the parameter estimates (Table 2).

Diagnostic plots for the final PK model are presented in Figure 1. The conditional weighted residual values (CWRES) were evenly distributed around zero (Figure 1A,1B), indicating no major bias in the structural model. The observed concentrations were evenly distributed around the line of identity, indicating that there was no bias in the population parameters (Figure 1C) and the structural model was appropriate for most individuals (Figure 1D). A prediction- and variability-corrected visual predictive check is presented in Figure 2. This plot shows that 176 of the 210 observed concentrations (83.8%) fell within the 80% prediction intervals, and the observed 10th, 50th, and 90th percentiles fell within the 95% confidence intervals (CIs) of the simulated 10th, 50th, and 90th percentiles. These results suggest that the final PK model appropriately describes the observed data and has acceptable predictive performance.

The currently recommended dosing regimen for patients with intra-abdominal infection by P. aeruginosa was optimal when a target was 40% fT_>MIC and a MIC was equal to or less than 2 mg/L (Figure 3). When a MIC was 4 mg/L, the PTA was close to 90%, but not reached. For the target 100% fT_>MIC, the recommended regimen was optimal when a MIC was equal to or less than 0.25 mg/L. For the target 100% fT_>4xMIC, this regimen did not attain 90% when a MIC was greater than 0.125 mg/L.

The optimal dosage regimen for the three treatment targets was investigated intensively under various conditions (Figure 4; Supplementary Tables S1A–S1C). For patients with normal renal function (CL_CR of 90–130 ml/min), a dosing regimen of 1 g q8h by i. v. infusion over 0.5 h was optimal when the target was 40% fT_>MIC (Figure 4; Supplementary Table S1A) and the MIC was equal to 4 mg/L; 2 g q8h over 3 h was optimal when the target was 100% fT_>MIC (Figure 4; Supplementary Table S1B) and the MIC was equal to 1 mg/L; and 2 g q8h over 3 h was optimal when the target was 100% fT_>4xMIC (Figure 4; Supplementary Table S1C) and the MIC was less than 0.5 mg/L. For patients with CL_CR of 10–25 ml/min, a dosing regimen of 2 g q8h by i. v. infusion over 3 h was optimal when the target was 100% fT_>4xMIC and the MIC was equal to 2 mg/L. For patients with CL_CR of 25–50 ml/min, a dosing regimen of 0.5 g q8h by i. v. infusion over 0.5 h was optimal (PTA ≥90%) when the target was 40% fT_>MIC and the MIC was equal to 4 mg/L. For patients with CL_CR of 50–90 ml/min, a dosing regimen of 1 g q8h by i. v. infusion over 3 h was optimal when the target was 100% fT_>MIC and the MIC was equal to 1 mg/L.

The ECMO flow rate, the only significant covariate for V_C, had a slight effect on the PTA (Supplementary Figure S2). For patients with CL_CR of 50–90 ml/min, a dosing regimen of 1 g q12h by i. v. infusion over 0.5 h was optimal for the MIC of 4 mg/L when the flow rate was 4 L/min or 6 L/min, while that was not optimal when the flow rate was 2 L/min.