Searches were performed in the National Library of Medicine up to February 16th, 2023. Although we decided to primarily use only one database, we followed the 2020 Preferred Reporting Items for Systematic Reviews and Meta-Analyses recommendations (32). The search string used was “empagliflozin AND pharmacokinetics.” Original reports with no date limits were contemplated. References listed within bibliographies of the included records and relevant articles known to the authors were also considered. Identified titles and abstracts were independently screened by two reviewers (SL and ADD) in an unblinded fashion. Full-text publications of candidate reports were reviewed in detail for eligibility. Any discrepancies were resolved by discussion among authors until consensus was reached.

We included original articles reporting population pharmacokinetic models of empagliflozin in humans. Reports in languages other than English, French, German, Italian, Portuguese, or Spanish were excluded. Similarly, articles reporting non-compartmental pharmacokinetic analyses, or not reporting a pharmacokinetic model (i.e., papers focusing exclusively on pharmacodynamics) were excluded. Data were collected according to a pre-defined checklist, including demographics, studied disease, sample size, investigated dose(s), modelling software employed, used validation methods, pharmacokinetic parameters, included covariates explaining interindividual variability, and residual error.

To qualify as a reference model for the subsequent modelling steps, the following criteria had to be fulfilled: (1) sufficient details provided in order to replicate the model, (2) model validation including visual predictive checks (VPCs), (3) accurate replication of modelling results on a virtual population with similar characteristics as the one on which the model was originally developed. Demographic and other relevant baseline characteristics (i.e., age, sex, weight, body mass index (BMI), and creatinine levels and/or estimated glomerular filtration rate (eGFR)) required to build the virtual population used in the evaluation of the predictive performance of the model were retrieved from the CDC-National Health and Nutrition Examination survey (NHANES) database (33).

Following model screening (34–39), the selected model was reparameterised taking into account the role of developmental growth on drug disposition (24).

To ensure the inclusion of relevant covariates during the assessment of the dose rationale for empagliflozin in children, anonymised clinical characteristics from a small representative sample of the heart failure paediatric population followed by some of us (Table 1) were used as basis for resampling and creation of a virtual population similar to the one to be included in the prospective phase 2a trial. Ultimately, as renal function maturation is complete within 1 to 2 years after birth, only body weight was identified as a significant covariate in this model. Consequently, a large virtual cohort was simulated based on normally distributed weights, with mean and distribution (38.1 ± 16.8 kg) reflecting the actual clinical population characteristics. Random sampling was then used to optimise sampling procedures, assuming study protocols with 12, 18 and 40 patients. Doses and dosing regimens were assessed based on weight bands, as a dose in mg/kg was deemed not feasible for the available dosage form (i.e., tablets).

A prospective trial has been proposed to investigate pharmacokinetics, ease of swallow, short-term safety and explore efficacy of empagliflozin among children 6 to < 18 years of age with heart failure, and Duchenne muscular dystrophy associated cardiomyopathy (ISRCTN12497973). This will be a prospective, open-label trial, in which children will take commercially available tablets of empagliflozin once daily for 6 months, with a full-day stay (max. 8 h) upon the first drug intake (Visit 1), and one opportunistic pharmacokinetic sample 1 week later (Visit 2). Further visits at 6 weeks (Visit 3), 3 (Visit 4) and 6 (Visit 5) months will assess the safety of empaglifozin and explore efficacy markers (Figure 2).

A preliminary assessment of the study feasibility indicates that recruitment has to be limited to a maximum of 12 participants. In addition, given the maximum blood volume allowed by current regulations (40, 41) and the need for safety monitoring blood draws, it will be possible to collect a maximum of 6 blood samples for pharmacokinetics at Visit 1, one opportunistic PK sample at Visit 2, and, potentially, one additional opportunistic PK sample at Visit 3. These restrictions represent an important risk, as both the dose rationale and eventual assessment of potential correlations between exposure and treatment response may not be established. In fact, there are various examples of inconclusive studies, which have relied on standard, empirical protocols, without any concern about the operating characteristics of the study design, other than a statistical power statement, which is often of limited value in exploratory studies, as outlined by the learning-confirming paradigm (42, 43). In a rare disease setting, the implementation of a clinical trial without a comprehensive evaluation of the informational content of a study protocol and its optimisation should be ethically questionable.

The $DESIGN feature in NONMEM^® was used to explore and optimise the proposed study design (44). This D-optimality based approach seeks to minimise the covariance of the parameter estimates, by assessing the Fisher information matrix and providing the expected model parameter precision and uncertainty (45, 46). The following constraints were set: n = 12 patients, maximum 6 samples on day 1, and 7 to 8 samples per patient over the whole trial, length of stay on day 1 of maximum 8 h, sampling time at Visit 2 (1 week after empagliflozin start) between 21 and 27 h after last drug intake. In addition, scenarios were assessed in which a single sampling schedule is used for all study participants, and individualised sampling schedule, or sampling matrices are evaluated. Scenarios with a total of 4, 5, and 6 samples were also explored (Supplementary Table 1). For the sampling matrices, two (i.e., 2 groups, 6 patients each) or four (i.e., 4 groups, 3 patients each) different sampling schedules have been considered. Moreover, alternative schedules were evaluated, including an additional opportunistic sample at Visit 3 (>3 weeks after empagliflozin start) in 50% or 100% of the participants, or to collect 8 samples in patients >20 kg body weight at Visit 1. Also, for comparison, parameter estimates and their precision were derived for a non-optimised schedule with 7 samples per patient at empirically chosen time points, for a rich sampling schedule with 12 participants contributing 13 samples each at empirically chosen time points (i.e., ethically not acceptable), and for a hypothetical study, albeit not feasible, with 40 participants contributing 13 samples each.

The D-optimisation step was implemented to identify potential designs in a more rapid and efficient manner than it would be achieved with traditional simulation re-estimation procedures. It also aimed at determining, for each of the explored sampling designs, the most informative sampling times. The impact of optimised designs on the accuracy and precision of the estimates was assessed by comparing the residual standard errors (%RSE) of the parameter estimates and the decrease in the objective function value (∆OFV).

The optimised sampling times identified in the previous step through the optimisation algorithm were used as input for four selected realistic scenarios (proposed sampling times were rounded to 5–15 min intervals according to clinical and practical considerations) with 6 samples at visit 1 based on a single sampling schedule for all participants, different sampling matrices with two or four groups as well as an individualised sampling schedule. Additionally, for comparison, we simulated a non-optimised, empirical sampling schedule; a rich sampling schedule with 13 samples; and a hypothetical design with 40 participants contributing 14 samples each (Table 2).

Given the rarity of the condition, informative and non-informative priors were used to support parameter estimation in the simulation re-estimation step, which included 500 trial replicates. The ratio between re-estimated and originally simulated pharmacokinetic parameter was calculated for each primary and secondary pharmacokinetic parameter of interest. A well-performing sampling scheme is expected to produce estimates close to 1.0 and with a relatively narrow distribution (i.e., 0.7–1.3). Whilst there is no specific range for defining accuracy, the proposed range takes into account the potential implications of varying exposure for the underlying pharmacological effect.

The selected sampling schedules were further evaluated based on a smaller sample size (i.e., 4, 6, 8 or 10 patients). Moreover, additional variants of this scenario were explored (patients staying only a maximum of 6 h at Visit 1, obtaining an opportunistic sample at Visit 3 from all or none of the participants, limiting the samples at Visit 1 to a maximum of 5, but obtaining an opportunistic sample at Visit 3 from all participants).

While summarising literature findings, weighted means were calculated. If needed, means and standard deviations were transformed into median and interquartile range by simulating a large normal sample distribution and calculating the pertaining summative statistics.

For the optimisation steps, the relative precision of secondary pharmacokinetic parameters was used as metrics for optimisation. Parameters of interest included the area under the concentration vs. time curve (AUC), the peak concentration (C_max), the time at which peak concentrations occur (T_max), and average steady-state concentration (C_ss). AUC and C_max were simulated for a virtual paediatric cohort of patients across a weight range between 10 and 90 kg, following administration of two dose levels (10 and 25 mg), which reflect the currently commercially available tablets of empagliflozin. Sampling times were handled as discrete data, where appropriate. Descriptive, summary statistics for continuous data included median and interquartile range, or mean and standard deviation. A two-sided t-test was used, with p-values of <0.05 as threshold for statistical significance.

Data handling, graphical and tabular summaries were performed in R statistics (v.4.2.2) with the support of R studio (2022.07.2 or higher) (47). Modelling and simulation steps were implemented in NONMEM^® v.7.5 (ICON Development Solutions, Ellicott City, United States), in combination with PsN^®-Toolkit (v. 5.3.0) and Pirana^® v.2.9.9.

The literature search process is summarised in Supplementary Figure 1. We identified n = 6 articles reporting two population pharmacokinetic models (34–39). Each of these was then replicated and/or adapted in slightly different populations, as summarised in Supplementary Table 2. They included 5,252 patients, mean 57.6 ± 10.9 years of age with type 1 (n = 1,316) or type 2 (n = 3,936) diabetes mellitus. Investigated doses ranged from 1 to 100 mg, from single dose up to 12 weeks once daily administration. All studies were performed using nonlinear mixed effects modelling, as implemented in NONMEM^®, four of them performed a bootstrap analysis and only two articles reported visual predictive checks. Five studies reported values for the primary and secondary pharmacokinetic parameters (35–39).

The selected model (which included details on the parameter estimates and evidence of acceptable predictive performance) was built on data from 56 type 1 diabetes mellitus patients who provided 1814 blood samples. These patients had a mean (sd) age of 41 ± 11 years, weight of 79.3 ± 14.3 kg, with normal renal function (36). Patients were given empagliflozin doses of 2.5 mg (n = 19), 10 mg (n = 19) or 25 mg (n = 18) once daily for 28 days. The model describes the pharmacokinetics of empagliflozin with 2-compartments and sequential zero- and first-order absorption, an absorption lag time to the central compartment, and a first-order elimination. Random effect parameters were identified for apparent clearance (CL/F), peripheral volume of distribution (V_p/F), intercompartmental clearance (Q/F), duration of the zero-order absorption (D1), absorption rate constant (k_a) and absorption lag time (ALAG1). Residual variability was described by a proportional error model. Body mass index was identified as a covariate explaining the interindividual variability of CL/F, V_c/F, V_p/F, and Q/F.

Prior to evaluating the effect of body weight on the disposition of empagliflozin, model performance was assessed by replicating the reported results in adults based on a VPC stratified by dose level (2.5, 10 mg, and 25 mg) (Supplementary Figure 2) (36). A similar assessment was performed after model reparameterisation according to allometric principles (Supplementary Figure 3).

The model including allometrically scaled parameters was subsequently used to simulate individual concentration vs time profiles in a virtual paediatric population (Figure 3) and support the selection of doses that ensure empagliflozin exposure that corresponds to the observed range in adults, or specifically to an exposure range that is ≥70% of the AUC obtained in normal weight adults with the standard 10 mg once daily regimen (i.e., recommended for adults with heart failure). Similar criteria were used to ensure safe exposure, defined as ≤ 130% of the AUC and C_max obtained in adults receiving the maximum recommended dose of 25 mg once daily (i.e., recommended, when needed, for adults with type 2 diabetes mellitus).

Whilst these thresholds may be somewhat arbitrary, they take into account known unexplained variability in drug disposition and formulation-related variation in exposure. Simulated AUC and C_max for virtual paediatric patients across a wide range of body weights (10, 15, 20, 30, 40, 50, 60, 70, 80 or 90 kg) are shown in Figure 4. From these results, it becomes evident that children below 15 kg should not be included in the study, because they would exceed the upper limit of exposure, as determined by the reference of 25 mg once daily in adults with type 2 diabetes mellitus. Median empagliflozin exposure ratio between adults and children ≥15 kg was within the target range. This suggests that the proposed regimen is expected to be safe, with the 75th percentiles of the exposure ratio remaining at <1.3 for both AUC and C_max (Table 3).

Given that only body weight was identified as a significant covariate, we generated a large virtual cohort with normally distributed weights, with mean (sd) and distribution (38.1 ± 16.8 kg) reflecting the actual clinical population characteristics. Random sampling was then used to optimise sampling procedures, assuming study protocols with 12, 18 and 40 patients. Decrease in the objective function value (∆OFV), and precision of the estimates were similar across the different scenarios (Supplementary Table 3), with only the scenarios with 3 and 4 samples (Optima 7 and 8) showing smaller ∆OFV and higher RSE. However, the precision of the interindividual variability parameters was poorly estimated. A hypothetically large study including 40 patients was required to allow adequate estimates of interindividual variability in clearance and volume of distribution. Yet, this remained tricky for intercompartmental clearance, given the trial constraints (max. 8 h at Visit 1, limited number of samples) (Supplementary Table 3).

The optimised sampling times to carry forward in the simulation re-estimation step were rounded taking into account feasibility criteria (Table 2). The optimisation algorithm was re-run, closely mirroring the original results (data not shown).

Twelve scenarios were simulated (Table 2; Supplementary Table 4). Interestingly, the empirical, rich sampling schedule including 12 samples on Day 1, and one opportunistic sample at Visit 2 and 3, did not perform very well. The ratio of re-estimated versus simulated parameter estimates included 1.0 in its interquartile range only for the peripheral volume of distribution, and with a pretty large spread. All the remaining pharmacokinetic parameters (CL, Q and AUC), whilst still showing a median within the pre-set acceptability range of 0.7–1.3, did not include 1.0 within their interquartile range (Figures 5c,d). The precision clearly improves with a hypothetical trial including 40 patients, each contributing with 14 samples (Figures 5a,b). However, these conditions are neither realistically feasible nor ethically allowed. On the other hand, the optimised sampling scheme with 12 patients divided into 4 different sampling groups, with both informative (Figure 5e) and, even more so, non-informative priors (Figure 5f), delivered better estimates than the empirical design with rich sampling. Of note is that only the central volume of distribution did not include 1.0 within its interquartile range [0.90 (IQR 0.89–0.92)]. The empirical, non-optimised sampling schedule with 7 samples (Supplementary Figures 4a,b), the optimised sampling schedule with 1 (Supplementary Figures 4c,d) and 2 sampling groups (Supplementary Figures 4e,f), and the individually optimised sampling schedule (i.e., one for each patient) showed lower performance (Supplementary Figures 4g,h). Specifically, these sampling matrices were clearly inferior to the selected sampling matrix with respect to central and peripheral volumes of distribution, with V_c falling outside the pre-defined 0.7–1.3 acceptability range (Supplementary Figure 4).

The optimised sampling schedule with 4 different sampling groups performed well also when considering the inclusion of just 8 patients instead of 12 (Figures 5g,h), whilst the reduction in sample size had a significant effect on the precision of estimates when a study including only 6 (two groups of 3 patients) or 4 patients (1 patient per sampling group) were considered.

Lastly, we explored the impact of not collecting any opportunistic sample at Visit 3 (Supplementary Figures 5a,b), performing it in all participants (Supplementary Figures 5c,d), performing it in all participants while limiting the samples at Visit 1 to five (Supplementary Figures 5e,f), as well as the possibility of limiting the hospital stay at Visit 1 to a maximum of 6 h (Supplementary Figures 5g,h). All these scenarios performed less well than the selected scenario with 4 sampling schedules (Figures 5e,f).