Operational Feasibility of A Priori Covariate-Based Meropenem Exposure Prediction Using a Southeast Asian Population Pharmacokinetic Model in Indonesian Critically Ill Patients with Sepsis

Study design and setting

This was a retrospective pharmacometrics cohort study using an in-silico approach at a tertiary academic centre in Yogyakarta with an active clinical pharmacy service and around 20-bed medical-surgical ICU. Data for January to December 2024 were drawn from the electronic hospital information system following approval from the institutional ethics committee, as further described in the Ethical Considerations section. As this was a retrospective study using routinely collected de-identified medical record data with no direct participant contact or intervention, the requirement for informed consent was waived by the institutional ethics committee. All simulations and analysis were run in RStudio 2025.09.1 on R version 4.5.1.

Study population

Adults (≥18 years) admitted with sepsis or septic shock (Sepsis-3) who received intravenous meropenem were eligible. Inclusion required documented sepsis or septic shock, intravenous meropenem for ≥24 h, a serum creatinine within ±24 h of initiation, and a culture with documented meropenem susceptibility (susceptible or resistant). Patients were excluded for carbapenem hypersensitivity, renal replacement therapy at initiation or within the first 72 h, pregnancy or breastfeeding, or death within 72 h from a non-infectious cause. Total sampling was used, and analytical feasibility was checked against the events-per-variable rule of ≥10 events per covariate.¹²

Population PK model and individual parameter estimation

Individual pharmacokinetic parameters were derived from the two-compartment linear-elimination model of Boonpeng et al.¹¹ (Songklanagarind Hospital ICU, Thailand; n = 52). The covariate equations applied to each patient i were:

CrCL was estimated with the Cockcroft–Gault equation using ideal body weight (Devine formula). Because no measured meropenem concentrations were available for Bayesian estimation, the individual random effect was fixed at η_i = 0, yielding covariate-driven typical-population (a priori) predictions. The reported between-subject variability for CL in the source model (60.5%) is therefore not propagated into the predictions, which is acknowledged as a source of overconfidence in the classification.

Covariate data and missing-value imputation

Serum albumin measured within ±48 h of meropenem initiation was available for 42 of 69 patients (60.9%); the remaining 27 (39.1%) were imputed by multiple imputation by chained equations (MICE; mice package, m = 20 imputed datasets, predictive mean matching) using age, sex, BMI, serum creatinine, septic shock status and baseline SOFA score as predictors. Estimates were pooled across the 20 datasets. Convergence and plausibility were verified by trace plots over 20 iterations and by density overlap of observed versus imputed values; the missingness pattern was most strongly associated with baseline SOFA score and serum creatinine, supporting the plausibility of a Missing at Random assumption conditional on observed covariates.

MIC assignment

Individual pathogen MICs were not recorded, so MICs were assigned from subpopulation MIC₅₀ values of the ATLAS surveillance database (Asia–Pacific, 2020–2024; adult ICU isolates). The MIC₅₀ of the relevant subpopulation was used for the primary analysis and the MIC₉₀ for the sensitivity analysis, computed separately within the susceptible and resistant subpopulations. For patients with multiple isolates, the highest MIC was used.

In silico simulation and target attainment

Steady-state concentration–time profiles were simulated with the deSolve package in R version 4.5.1 (RStudio 2025.09.1) using each patient’s actual prescribed regimen (dose, infusion duration and dosing interval). Profiles were integrated to steady state (≈ 5 × terminal half-life) at 0.01-h resolution. The free fraction was taken as f_u = 0.98 (free concentration = total × 0.98), consistent with the minimal protein binding of meropenem.¹³ The target of 100% fT > MIC required the free concentration to remain above the assigned MIC for the entire dosing interval; patients were classified as adequate (≥100% fT > MIC) or inadequate (<100%).

Clinical outcomes

The primary outcome was a composite 72-h treatment failure: any of antibiotic escalation (a broader agent added or substituted, or a documented meropenem dose increase), death within 72 h, fever (≥38.0°C) at hour 72, or worsening organ dysfunction (ΔSOFA ≥2 from baseline at hour 72). For patients who died of infection before hour 72, the mortality component was recorded, and the others were not assessed.

Statistical analysis

Continuous variables were tested for normality (Shapiro–Wilk) and summarized as median [IQR]; categorical variables as n (%). Patient characteristics were compared with the source population descriptively only (medians/proportions), with no hypothesis testing. The 95% CI for target attainment used the exact

Clopper–Pearson method. Between-subgroup comparisons (adequate versus inadequate; imputed versus observed albumin) used the Mann–Whitney U or t-test for continuous variables and Fisher’s exact or Pearson χ² for categorical variables. Risk ratios with 95% CIs were obtained with epitools and validated by non-parametric bootstrap resampling (B = 2000) with bias-corrected and accelerated (BCa) intervals. Multivariable analysis used Firth penalized logistic regression (logistf ) adjusted for baseline SOFA score and septic shock status. A two-sided P < 0.05 was considered significant. Analyses used the mice, epitools, boot, logistf and car packages in R.

Patient characteristics

Of 292 adults admitted with sepsis or septic shock in 2024, 221 (75.7%) received meropenem, and 69 met all criteria ( Figure 1). The commonest reason for exclusion was the lack of a culture with susceptibility testing (n = 61), which reflects routine ICU practice in Indonesia, where broad-spectrum empirical therapy usually precedes microbiological confirmation.

Figure 1. Patient selection flowchart.

Patients were screened from admission between January and December 2024. ICU, intensive care unit.

The cohort was mostly male (59.4%), with a median age of 67 years [IQR 60.0–73.0]. The findings reflect the older ICU study population. Illness was severe (median APACHE II 24 [20–28]; median baseline SOFA 9^8–12); 52.2% had septic shock and 79.7% needed mechanical ventilation. Most had reduced renal function (median CrCL 45.5 mL/min [29.1–75.5]), and only 2.9% met ARC criteria. Hypoalbuminaemia was common (median albumin 2.5 g/dL [2.0–2.9]).

Against the Boonpeng et al.¹¹ source population, the covariates that drive the model, such as CrCL (45.5 versus 44.6 mL/min), serum albumin (2.5 versus 2.4 g/dL) and proportion male (59.4% versus 59.6%) matched closely, which supports using the external model. Larger differences appeared in body weight (50.6 versus 61.5 kg), BMI (19.9 versus 22.9 kg/m²) and severity (APACHE II 24 versus 20; septic shock 52.2% versus 34.6%; Table 1). Body weight is not in the model and does not affect predictions. However, the greater severity, implies wider between-patient variability through the shock term on Vp.

Individual PK parameters

Individual estimates varied widely, matching the clinical mix ( Table 2). Median clearance was 4.79 L/h [3.44–7.27], spanning more than six-fold (2.21–13.58 L/h) with renal function, which ranged across CrCL 14.3–152.0 mL/min. This covariate was consistent with the linear covariate structure of the source model. The terminal half-life ran from 1.67 to 11.58 h, capturing both rapid elimination near normal renal function and slow clearance in severe AKI. Central volume followed albumin (median Vc 11.10 L [10.49–11.67]; range 8.27–14.50 L), and peripheral volume rose with septic shock, typical Vp 10.3 L without shock versus 21.6 L with it, more than doubling the distribution volume.

Target attainment

Under their actual regimens, 54 patients (78.3%; 95% CI 66.7–87.3%) were predicted to reach 100% fT>MIC and 15 (21.7%; 95% CI 12.7–33.3%) to fall short ( Table 3). The %fT>MIC distribution was sharply bimodal, with patients

clustered near 100% or near 0% and only four in between (24.6%, 26.9%, 54.4% and 71.9%; Figure 2). The split came from MIC imputation: all 15 inadequate patients had been assigned an MIC of 32 mg/L, the resistant-subpopulation MIC₅₀, mostly carbapenem-resistant A. baumannii (CRAB), with carbapenem-resistant P. aeruginosa and K. pneumoniae also present. Susceptible isolates (MIC ≤2 mg/L) all reached 100% fT>MIC, helped by the extended-infusion regimens used in over 70% of patients.

Figure 2. Distribution of steady-state %fT > MIC across all 69 patients.

The histogram is bimodal, with most patients clustering near 100% (adequate; susceptible MIC ≤2 mg/L) or near 0% (inadequate; resistant MIC 32 mg/L). The dashed line marks the 100% fT > MIC threshold.

This sets up everything that follows in this dataset the adequate or inadequate split is effectively a susceptible–resistant split, not a gradient of exposure against the same pathogen.

Association with clinical outcomes

In-hospital mortality was high overall (47/69; 68.1%), in keeping with national reports and well above the global average; 72-h treatment failure occurred in 38/69 (55.1%). For the primary outcome, failure affected 9/15 (60.0%) of the inadequate group and 29/54 (53.7%) of the adequate group—a relative risk of 1.16 (95% CI 0.58–2.29; P = 0.665). Bootstrap resampling agreed (RR 1.12; BCa 95% CI 0.60–1.73; P = 0.659). After Firth regression for baseline SOFA and septic shock (EPV = 12.7; VIF 1.03–2.02), the aOR stayed non-significant (1.24; 95% CI 0.40–4.04; P = 0.717), and neither covariate predicted the outcome ( Table 4).

Breaking the composite into its four parts: antibiotic escalation, 72-h mortality, fever at 72 h and ΔSOFA ≥2, showed no single component carrying an association ( Figure 3a), so a real signal was not being diluted within the composite. In-hospital mortality was likewise non-significant, and its point estimate pointed the other way (RR 0.74; 95% CI 0.35–1.56; P = 0.535; Figure 3b) for reasons set out in the Discussion, is not a protective effect of underdosing.

Figure 3. Forest plots of the association between predicted 100% fT > MIC attainment and (a) the composite 72-h treatment-failure outcome and its four components, and (b) in-hospital mortality.

Estimates are relative risks with 95% CIs; the composite primary outcome is shown in bold. The dashed line at 1.0 indicates no association.

Sensitivity analysis

The three pre-specified sensitivity analysis are in Table 5. Replacing MIC₅₀ with MIC₉₀ shifted the inadequate proportion by only 1.5 points and left both associations non-significant and in the same direction (failure RR 1.26; mortality RR 0.65). Restricting to the 42 patients with observed albumin widened the intervals, as expected from losing 39.1% of the sample, but kept the direction and the non-significance (failure RR 1.30, 95% CI 0.44–3.81; mortality RR 0.83, 95% CI 0.26–2.65).

The combination-therapy analysis was the most telling. Combination therapy went to 6 of 15 inadequate patients (40.0%) but only 7 of 54 adequate patients (13.0%) (OR 4.35; 95% CI 1.15–18.97; P = 0.028), and the inadequate patients on combinations fared worst (5/6 failed; 83.3%). Among the 56 patients on meropenem monotherapy, the failure RR moved from 1.16 to 0.84 and the mortality RR from 0.74 to 0.54, both still non-significant. That reversal fits confounding by indication which is clinicians steered combinations and closer monitoring to the sickest patients with resistant organisms rather than any true effect of meropenem exposure, and it strengthens rather than weakens the main finding.

Limitations and recommendations

The study has clear limits. The most substantial limitation is the lack of measured meropenem concentrations, which rules out Bayesian updating. The predictions are covariate-driven typical-population estimates rather than true individual forecasts, and fixing η = 0 ignores the between-patient variability the covariates do not explain (60.5% for CL in the source model), so the classification is overconfident. Second, drawing MIC from surveillance data rather than measuring it per isolate turned the key analytical variable into a susceptibility surrogate. Third, although common and defensible, imputing albumin adds uncertainty to the volume estimates and a residual Missing Not at Random mechanism cannot be formally excluded given that albumin measurement was not part of the routine laboratory panel at this institution. Fourth, with 15 inadequate patients and a 60% event rate, the multivariable model was underpowered. Consequently, the aOR interval (0.40–4.04) runs from real protection to real harm, meaning the analysis cannot exclude either direction. This uncertainty is an inherent consequence of the sample size rather than a methodological flaw, and it reinforces the need for larger prospective studies before any inference about the direction of the exposure–outcome relationship can be drawn. Fifth, the single-centre design and the exclusion of RRT patients limit the generalizability of these findings. Additionally, while body weight is not a covariate in the model, the observed differences in weight and BMI between our cohort and the source population may reflect broader phenotypic differences that could affect model performance, underscoring the need for local model development or formal external validation.

These gaps define priorities for future research. Prospective studies should incorporate TDM with sparse sampling to enable Bayesian estimation and robust external validation; determine isolate-specific MIC to evaluate exposure against actual pharmacodynamic targets; recruit larger, preferably multicentre cohorts stratified by pathogen susceptibility to ensure that comparisons reflect drug exposure rather than differences in causative organisms; and evaluate the optimal loading dose for patients with septic shock and impaired renal function. If existing models prove unsuitable for this population, developing an Indonesian population pharmacokinetic (popPK) model should be a priority.

Implications of this study

In settings where therapeutic drug monitoring (TDM) is unavailable, the covariate-based simulation approach may provide a pragmatic strategy for evaluating dosing adequacy, although it cannot replace TDM or Bayesian forecasting. With independent validation, this framework may also be extended to other drugs with narrow therapeutic windows. Importantly, this study highlights that successful implementation of model-informed precision dosing (MIPD) depends not only on pharmacokinetic models but also on the availability of reliable pathogen minimum inhibitory concentration (MIC) data, underscoring the need to strengthen microbiology support within antimicrobial stewardship programs. Finally, these findings provide preliminary evidence to inform the implementation of TDM services in Indonesian tertiary hospitals, in alignment with the National Action Plan for Antimicrobial Resistance Containment.