The data used in this work were collected in a study comparing the WBPK of radiolabeled glyburide ([11C]glyburide) in the absence or the presence of rifampicin, a potent OATP inhibitor [9]. The study was designed to investigate the impact of age and sex on OATP function. Therefore, 16 healthy volunteers (5 females and 11 males) were included to compare males > 50 years old with males < 30 years old (age effect) and males > 50 years old with females > 50 years old (sex effect). Table 1 presents the demographic characteristics of each group.

All subjects underwent a baseline WB4D PET acquisition, which started with the injection of a microdose of [11C]glyburide (intravenous [i.v.] bolus). In total, 10 of the 16 subjects underwent a second [11C]glyburide scan in the presence of OATP inhibition, at least 3 h after the end of the first acquisition to allow for decay of carbon-11 radioactivity [10] (half-life = 20.4 min). To this end, an infusion of rifampicin (9 mg/kg diluted in glucose 5% perfused over 30 min) was administered immediately before [11C]glyburide injection and WB4D PET acquisition. Potential carryover of drug within occasion was not taken into account, as the radioactivity from the first dose had completely disappeared and the amount of residual drug was negligible.

The study protocol was approved by an ethics committee (CPP IDF5: 17041, study reg. no. EudraCT 2017-001703-69). Experiments were conducted with respect to the 1975 Declaration of Helsinki. Written informed consent was provided after the subjects had received a medical examination and been given all information about the study. A blood sample was collected for pharmacogenetic analysis of OATP transporter genes. Polymorphism of the SLCO1B1 gene coding for OATP1B1 (c.521T>C genotype, rs4149015) is relatively common in Caucasian and Asian populations, with a frequency of about 15–20%, but rare in subjects of African origin [11]. Homozygous carriers of the C allele (CC) show lower elimination rates, higher plasma exposure, and enhanced toxicity risk to OATP1B1 substrates such as statins or some antidiabetic drugs [12, 13]. Consequently, the subject harboring the CC allele in our study showed a ~60% increase in plasma exposure to [11C]glyburide associated with a ~48% decrease in its liver-uptake transfer rate) and a lower response to rifampicin [6]. Therefore, this subject was excluded from model building to avoid skewing the results with one outlier.

WB4D PET acquisitions were performed using a Signa positron emission tomography–magnetic resonance imaging (PET–MR) scanner (GE Healthcare, USA). First, dynamic mono-bed acquisition focused on the abdomen was performed for 3 min after the injection of [11C]glyburide to capture early time-points (16 frames of 10 s) in the Liver, spleen, pancreas, kidneys, and myocardium. Then, 15 whole-body PET images (five bed positions, every 2.5 min) were acquired over 37 min (Fig. 1).

Whole-body PET images of [11C]glyburide distribution in healthy volunteers, without and with pre-infusion of rifampicin. Three subjects were randomly selected, one from each group: subject 14 (group 1), subject 9 (group 2), and subject 3 (group 3).

PET images were reconstructed using a three-dimensional (3D) iterative reconstruction algorithm, and WB4D PET images were generated to capture the kinetics of radioactivity in the different organs over 40 min. Volumes of interest (VOIs) were manually delineated for each patient for selected tissues (liver, kidneys, spleen, myocardium, pancreas, brain, testis, muscle, and eyes) using PMOD® software (version 3.9) on the last images of the acquisition and used to compute the radioactivity per mL. An additional VOI was delineated over the left ventricle and the aorta representing arterial blood pool. Kinetics in this VOI have been previously validated as an image-derived input for [11C]glyburide, which is a metabolically stable probe.

Concentrations of radioactivity expressed in kBq/mL were extracted. Measurements were then converted for pharmacokinetic modeling to concentrations in \(\upmu\) g/mL using radioactivity expressed as concentrations (\(_}\) in kBq/mL), specific radioactivity (SRA; in kBq/\(\upmu\) mol), and the molar mass of glyburide (\(_}\)) is 494 \(\upmu\) g/\(\upmu\) mol.

$$C=\frac_}}}\times _}$$

(1)

The same conversion was used to convert the individually injected radioactivity to a dose in grams.

Supplementary Fig. S1 shows the total radioactivity, expressed in percentage of injected dose corrected by radioactive decay, as a function of time in the different organs. From these data, we decided to model only the organs in which most of the radioactivity was found to be concentrated. Brain, muscle, and eye tissue were not included in the dataset for modeling because the levels of radioactivity in these organs were negligible. Bladder and gallbladder presented very different profiles, with a marked accumulation with time and very high interindividual variability. As attempts to include them in the model were unsuccessful, we present the modeling performed for only the following five organs: blood, kidney, liver, spleen, as well as pancreas as the active site of action of glyburide.

We denote \(_}\) the \(}^}\) concentration observed in subject i at time tijko in tissue k, where k denotes the tissue by its first letter, i.e., B(lood), L(iver), K(kidney), P(ancreas), and S(pleen), at occasion o (= 1 for [11C]glyburide alone and 2 when [11C]glyburide was administered after rifampicin infusion). To account for intervariability (IIV) across subjects (i = 1,..N) and interoccasion (IOV) variability over the two acquisition periods, we describe observed concentrations through the following non-linear mixed-effect model:

$$_}=_}(_},_})+_}(_},_},_})_}$$

(2)

$$_}=H(\mu ,_},_},_},\beta )$$

\(_}\) describes the structural model associated with tissue k, and gk quantifies the variance of the residual error for this tissue, with \(_}\) a normally distributed random intervariable with a variance of 1. For each tissue, the error model was either a constant (\(_}=_}\)), proportional to the predicted concentration (\(_}=_}_}\)) or a combined error model (\(_}=_}+_}_}\)). \(_}\) will denote the vector of parameters of the error model, \(_}=(_},_})\). We assumed that the vector \(_}\) of individual parameters was log-normally distributed by modeling it as the exponential of a linear combination of fixed effects \(\mu\), covariate effects \(\beta\) multiplying known individual covariates \(_}\), individual random effects \(_}\), and occasion-level random effects \(_}\) (H = exponential in equation 2). We assumed multivariate normal distributions for \(_}\) and \(_}\) with variance–covariance matrices \(\Omega\) and \(\Gamma\), respectively (\(_}\sim N(0,\Omega )\), IIV, and \(_}\sim N(0,\Gamma )\), IOV).

We introduce the covariate effect in the model. We investigated the effect of weight, age (known to affect liver and kidney function), sex (which can influence drug exposure through differences in enzyme activity and hormone levels), and rifampicin coadministration (which can affect hepatic transport). Covariates were introduced in the model as:

$$_,\text}=_}+_}\text\left(}_,\text}\right)+_}+_,\text}$$

\(\text\left(_,\text}\right) =^\left(_,\text})-^(_}\right)\) if the covariate is continuous and \(\text\left(_,\text}\right) =__,\text}= 1\}}\) if the covariate is categorical. Here, \(_}\) is the vector of individual parameters of subject i for occasion o, \(\left(_,\text}\right)\) is the corresponding vector of covariates, \(_} \sim \mathcal(0, \Omega )\) (respectively \(_,\text} \sim \mathcal(0, \Gamma )\)) is the vector of random effect, which describes the interindividual variability (respectively intraindividual variability), and \(_\) is the effect of covariate c on parameter \(\psi\).

For weight, \(_}\) was fixed according to the allometric model to 1 for weight.

We denote \(\theta =(\mu ,\Omega ,\kappa ,\sigma )\) the population parameters of the model to estimate. Thereafter, we used the Stochastic Approximation Expectation-Maximization (SAEM) algorithm for parameter estimation [14]. Parameter uncertainty was estimated using the Fisher information matrix, computed using a first-order linearization of the model.

Each tissue was represented by at least one compartment in equilibrium with the blood compartment. Some tissues could be represented by subcompartments, one compartment connected to the blood and others connected to this compartment.

In the first step, the structural model was built using the data collected after [11C]glyburide was administered alone. Each tissue was first modeled separately. We used the log-likelihood ratio test (LRT) to determine the appropriate number of compartments and to obtain initial parameter estimates. For model building purposes, we excluded one woman who had a CC mutation on gene SLCO1B1, as she presented a markedly different kinetic profile with much higher blood concentrations and much lower liver concentrations. In addition, we considered the splenectomy that one man had undergone by considering this tissue was missing for this subject. A joint model was then fitted simultaneously including the data for all five tissues, with a combined error model for all tissues. For each tissue, we selected the number of compartments, testing linear and non-linear processes for distribution and elimination, and finally chose the residual error model by testing additive and proportional error models separately for each tissue. For the pancreas and spleen, which showed rapid distribution, we also tested a simple model with a linear relationship between the tissue and blood concentrations, representing a fast equilibrium. Indeed, if the rate \(_}\) (respectively \(_}\)) is much smaller than \(_}\) (respectively \(_}\)), we can model the pancreas (respectively spleen) concentration as a ratio of the blood concentration. For the liver and kidney, we also tested a two-compartment structure with a superficial compartment exchanging with blood and a deep compartment exchanging with the first (Fig. 2). The model corresponding to the best corrected Bayesian information criterion (BICc) was selected:

$$}_}=-2}_}(\hat)+\text(_})\text(N)+\text(_})\text(_})$$

(3)

where \(\hat\) is the maximum likelihood estimate of \(\theta\) that maximizes the log-likelihood function \(}_}(\theta )\), \(_}\) represents the variances of the random effects, \(_}\) represents the fixed effects, and \(_}\) represents the total number of observations [15].

Diagram of the method. Each tissue is represented by a colored background and a letter that is used to index the volume of the corresponding compartment: blood (B, red), liver (L, yellow), kidneys (K, purple), pancreas (P, pink), and spleen (S, blue). Each block represents a modeling step. The text in gray represents the elements tested, while the text in bold black represents the elements included in the final model. B blood, L liver, K kidneys, S spleen, P pancreas, V volume, Cc concentration, k exchange rate between tissues, kₑ elimination rate, BW body weight, Rif rifampicin, COSSAC Conditional Sampling Use for Stepwise Approach Based on Correlation Tests, IIV interindividual variability, IOV interoccasion variability

We then tested the effect of body weight (BW) on the volumes of distribution using an allometric model. Allometric scaling was kept if the BICc did degrade by > 3 points.

In the second step, the model selected using [11C]glyburide alone was applied to the full dataset, including both administrations. Interoccasion variability (IOV) was introduced on all parameters (transfer rate constants between tissues, and volumes of distribution). We performed an exploratory analysis of parameter–covariate relationships using the Conditional Sampling Use for Stepwise Approach Based on Correlation Tests (COSSAC) algorithm [16], considering as potential covariates rifampicin, age, gender, and combined age-gender groups. The COSSAC algorithm uses the correlation between conditional samples of individual parameters and covariates to include the most relevant parameter–covariate relationship iteratively. The covariate model can include continuous covariates (such as body weight or age) and categorical covariates (such as sex, group, or rifampicin effect). The model with IOV on all parameters proved unstable. We therefore used a stepwise bottom-up approach for covariate building, including the effect of rifampicin first. For this, we introduced IOV only on the parameters associated with the liver and kidney (transfer rate constants between blood and liver or kidney), adding a fixed effect corresponding to a systematic change in these parameters when coadministered with rifampicin. Those relationships were tested one by one, removing both the IOV and the fixed rifampicin effect on parameters when the relative standard error (RSE) on the fixed effect exceeded 100%. The IOV for the effects remaining in the model was then fixed to 0.01 to help parameter identification. At this stage, we refined the covariance structure: We removed interindividual variabilities if their relative standard error (RSE) exceeded 100%, and we tested for correlations between parameters.

In a final step, we tested each covariate highlighted by the exploratory COSSAC proposal one by one and kept the covariate if the BICc improved. The error models were also refined by testing proportional, combined, or additive error models for each compartment.

The final model was evaluated using graphical diagnostics: prediction-corrected visual predictive check (pcVPC) [17] and normalized prediction distribution errors (NPDE) [18]. The pcVPC compares empirical percentiles (observed data) with theoretical percentiles and prediction intervals computed using data simulated under the estimated model. In our case, we chose to compare the 10th, 50th, and 90th percentiles (80% confidence interval [CI]). The corrected predictions take into account the large variability in actual dose received, a variability built in the computation of the NPDE. The NPDE was plotted against time and against predicted concentrations and can be interpreted as residuals for non-linear mixed-effect models. The goal is to align the empirical percentiles with the theoretical percentiles.

Convergence assessment was performed at key steps in the analysis to calibrate the number of iterations and chains. We ran five models starting from different initial estimates (within 20% of the final estimates) and assessed the stability of parameter and likelihood estimates by visual diagnostics of the confidence intervals (CI).

In a sensitivity analysis, we applied the final model to the full dataset including the subject with a mutation on SLCO1B1.

We used Python for data management and analysis, and Monolix-2023R1 for population parameter estimates. Algorithm settings included four chains, and the minimum number of iterations was 300. The maximum number of iterations was chosen for stability, starting with 500 iterations and increasing to 1000 if run assessments showed instability in parameter estimates.