A literature search was performed in PubMed in March 2023 to retrieve PK models published after 2010 that could be employed for simulations of both pharmacogenetic and TDM-guided dosing. The search strategy included relevant single-nucleotide polymorphisms (SNPs) known to exert clinically significant effects on PK, as defined by the Dutch Pharmacogenetic Working Group and the Clinical Pharmacogenetics Implementation Consortium. The specific genes included in the search were ABCG2, CYP2B6, CYP2C9, CYP2C19, CYP2D6, CYP3A4, CYP3A5, DPD, MTHFR, NUDT15, SLCO1B1 (OATP1B1), TPMT, and UGT1A1.

The initial screening of publication titles and abstracts was performed using Rayyan QCRI by author NR, and an independent review of the results was conducted by author MC. The inclusion criteria consisted of (1) identification of PK models based on the search terms, (2) determination of whether SNPs relevant to drug exposure were incorporated as covariates within the model structure, and (3) selection of models for drugs used in TDM. In this context, TDM was defined as the direct measurement of drug concentrations without any model-based estimations, such as model-informed precision dosing (MIPD). Exclusion criteria for the identified publications included non-English language publications, physiology-based PK models, drugs not commonly used in Dutch or Swedish clinics, and models based on fewer than 30 patients.

In cases where multiple models for a specific drug were available, the model with the highest number of samples and patients, along with its potential for accurate predictions and clinical relevance, as assessed by both NR and MC, was chosen. All selected models were then applied in PK simulations to evaluate their concordance with reported observational PK values. If necessary, adjustments were made to the models to align them with the research focus.

The selected models and corresponding parameter estimates were coded in mrgsolve (version 0.11.2) in RStudio (version 4.1.1.) [10,11,12]. For each drug case, a virtual population of 1000 patients was created. Pharmacogenetic phenotypes were randomly sampled according to the original patient distribution or, in case this was not available, from the Genome Aggregation Database (gnomAD) [13]. For simplicity, other covariates were fixed to a single value (e.g., weight and co-medication): continuous variables at the median value; for categorical variables, the reference value was fixed to the mode. By comparing patients to their own baseline data across simulations that only incorporated variations in IIV, IOV, and RUV, the significance of having varying covariates among patients was less critical in this analysis. The dosing regimen used in the simulations was the one specified in the original publication of the model. In the situation that the model was based on multiple dosing categories, the median dose was selected.

Simulations for each drug case were performed in two steps. During simulation step 1 (SIM1), a single model-based simulation was performed using the reported typical and variability parameters. Each PK parameter for individual i (\(_\), e.g., clearance) is here represented by a fixed effect component (\(_\), i.e., the typical value), multiplied by a random effect component representing the deviation of \(_\) from \(_\) via a log-normal distribution. However, rather than sampling the random effect component (\(_\)) from a normal distribution with a mean of 0 and the reported standard deviation of ω, the distribution was changed to a mean of 0 and standard deviation of 1. Following on, \(_\) was multiplied by θ1,IIV with the value of ω to capture the original deviation of \(_\) from the typical value \(_\) (Eq. 1).

$$\theta_ = \theta_ \times e^}}} \times \eta_ }} .$$

(1)

Each PK parameter for individual i of time j (\(_\)) is then represented by \(_\) multiplied by a random effect component, representing the deviation of \(_\) from \(_\) via a log-normal distribution. However, rather than sampling the random effect component (\(_\)) from a normal distribution with a mean of 0 and the reported standard deviation of π, the distribution was changed to a mean of 0 and standard deviation of 1. Following on, \(_\) was multiplied by \(_}\) with the value of π to capture the original deviation of \(_\) from \(_\) (Eq. 2).

$$\theta_ = \theta_ \times e^}} \times \kappa_ }} .$$

(2)

The influence of a covariate value of individual i (\(}_\)), such as the influence of pharmacogenetic CYP3A5 variations (\(}_3\text5,i}\)) on drug clearance, was added as a relative deviation from \(_\) for each affected parameter unless otherwise specified (Eq. 3).

$$\theta_ = \theta_ \times e^}}} \times \eta_ + \theta_}} \times \kappa_ }} \times }_ .$$

(3)

The difference between the observed concentration (\(}_\)) and model-predicted individual concentration for individual i at time j (\(}_\)) is represented by a random effect component, either additively (additive error) or in relation to the magnitude of the \(}_\) (proportional error). However, the distribution of these random effects (\(_)\) was changed to a mean of 0 and standard deviation of 1 (instead of the reported standard deviation of σ). Following on, \(_\) was multiplied by \(_}\) with the value of σ to capture the original difference between \(}_\) and \(}_\) (Eq. 4, additive error; Eq. 5, proportional error).

$$}_ = }_ + \theta_}1}} \times \varepsilon_$$

(4)

$$}_ = }_ + }_ \times \left( }2}} \times \varepsilon_ } \right).$$

(5)

Because all random effect distributions were fixed to a mean of 0 and a standard deviation of 1, the simulated ηi,κij, and εij values represent the individual z-scores for each PK parameter (i.e., number of standard deviations from the mean). These individual z-scores were exported for simulation step 2 (SIM2), along with individual values of \(_\) and \(}_1i}\).

In SIM2, the magnitudes of IIV, IOV, and RUV were varied (i.e., by changing the \(_}\), \(_}\), and \(_1/2}\) values from Eqs. 1, 2, 3, 4, and 5, respectively) in each simulation. All other conditions remained the same, including the patient population, their individual z-scores, and the time point at which exposure was measured. This was done to evaluate their influence on the predicted drug exposure by the pharmacogenetic and TDM methods. IIV and IOV were only changed for the parameter of interest (e.g., clearance), while the original variability values were kept for all other parameters. The values of IIV, IOV, and RUV were varied in each simulation using different combinations of magnitudes, ranging from 0 to 1 in 0.1 increments for IIV, ranging from 0 to 0.5 in 0.1 increments for IOV, and ranging from 0 to 0.3 in 0.05 increments for RUV (proportional), unless otherwise specified. For the additive error model, a range between 0 and 0.5 of the average concentration achieved during steady state (Cavg,ss) was evaluated, in 0.05 increments. This resulted in 463 simulated scenarios (i.e., 6 × 7 × 11 + the originally reported IIV, IOV, and RUV values) for each of the original 1000 individuals, per drug case.

For the parameter of interest (e.g., clearance), the influence of the covariate of interest (e.g., CYP450 genotype/phenotype) was scaled to the original value (\(}_1i}\) from SIM1). This was done to ensure that the individual parameter value remained identical, whilst only the relative contribution of the random effect components (\(_} \times __}\times _\)) and the fixed effect component (i.e., \(}_2i}\)) differed between the subsequent simulation rounds. As such, the total variability between patients remained consistent throughout the simulations and only the relative contribution of unexplained IIV (through θIIV × ηi) and explained IIV (through the COVi) changed. The covariate factor for each simulation COVSIM2 was calculated as follows (Eq. 6):

$$}_}2i}} = \frac},}1}} \times \eta i + \theta_},}1}} \times \kappa ij} \right)}} }}},}2}} \times \eta i + \theta_},}2}} \times \kappa ij} \right)}} }} x }_}1i}} .$$

(6)

For each simulated individual (n = 1000) and under each simulation scenario (n = 463), the three different values (i–iii) of the relevant summary exposures (Ctrough or AUC, depending on the drug) were exported (Fig. 1):

(i) The true value, including the impact of all covariates (i.e., explained IIV) and unexplained IIV (Eqs. 1 and 3), but without the effect of IOV and RUV (i.e.,\(_}\text_}\) = 0)

(ii) The pharmacogenetic predicted value, including the impact of all covariates (i.e., explained IIV) (Eq. 3), but without the effect of unexplained IIV, IOV, and RUV (i.e., \(_}, _},\text_}\) = 0)

(iii) The TDM predicted value, including the impact of all covariates (i.e., explained IIV), unexplained IIV, IOV, and RUV (Eqs. 1, 2, 3, 4, and 5).

Ctrough was here represented as a single trough concentration at steady state. For TDM of the AUC, the exposure was calculated after a single dose by the trapezoidal method (NCA) using three optimized time points at a single dose interval (i.e., a limited sampling approach) as earlier proposed [14, 15]. IOV was excluded from the (i) true value as it was considered a random value that cannot be predicted for an individual from one occasion to another, consistent with previous research [16].

Summary exposures, as predicted by pharmacogenetic and concentration measurements, were compared to the “true” PK value (Eq. 7).

$$} = }_}}} /}_}}} .$$

(7)

For each individual, the approach (pharmacogenetic vs. TDM) that predicted an accuracy closest to 1, meaning that the predicted exposure is closest to the true exposure, was considered superior. This study did not explore dose adjustments based on the TDM and pharmacogenetic testing results. The percentage of the 1000 patients for whom the pharmacogenetic biomarker performed better than TDM was computed for all simulated scenarios. The final results were visualized in 3D plots using R package plot3D (version 1.4).

In order to address concerns associated with the calculation of AUC using NCA [17], the potential for improved accuracy offered by MIPD was assessed for drugs in which AUC had been suggested to be the relevant exposure metric (specifically, vincristine, and 5-FU). For MIPD, the same three concentration measurements utilized in the NCA were inputted into the corresponding PK models to estimate individual AUC values in NONMEM (version 7.4). The AUC values obtained through both NCA and MIPD were compared to the true values in order to evaluate their accuracy:

$$}_},i/}, i}} = \frac}_}, i/},i}} }}}_},i}} }}.$$

(8)