Evaluation of body-surface-area adjusted dosing of high-dose methotrexate by population pharmacokinetics in a large cohort of cancer patients

The aim of this study was to identify sources of variability including patient gender and body surface area (BSA) in pharmacokinetic (PK) exposure for high-dose methotrexate (MTX) continuous infusion in a large cohort of patients with hematological and solid malignancies. We conducted a retrospective PK analysis of MTX plasma concentration data from hematological/oncological patients treated at the University Hospital of Cologne between 2005 and 2018. Nonlinear mixed effects modeling was performed. Covariate data on patient demographics and clinical chemistry parameters was incorporated to assess relationships with PK parameters. Simulations were conducted to compare exposure and probability of target attainment (PTA) under BSA adjusted, flat and stratified dosing regimens. Plasma concentration over time data (2182 measurements) from therapeutic drug monitoring from 229 patients was available. PK of MTX were best described by a three-compartment model. Values for clearance (CL) of 4.33 [2.95–5.92] L h− 1 and central volume of distribution of 4.29 [1.81–7.33] L were estimated. An inter-occasion variability of 23.1% (coefficient of variation) and an inter-individual variability of 29.7% were associated to CL, which was 16 [7–25] % lower in women. Serum creatinine, patient age, sex and BSA were significantly related to CL of MTX. Simulations suggested that differences in PTA between flat and BSA-based dosing were marginal, with stratified dosing performing best overall. A dosing scheme with doses stratified across BSA quartiles is suggested to optimize target exposure attainment. Influence of patient sex on CL of MTX is present but small in magnitude.


Introduction
Methotrexate (MTX) is considered an efficacious, costeffective and acceptably safe drug for the treatment of many hematological/oncological disorders and autoimmune diseases [1]. The folate analogue MTX acts as an antineoplastic agent via competitive inhibition of dihydrofolate dehydrogenase, resulting in depletion of purines and thymidylate leading to impairment of DNA synthesis [2,3]. The drug can be administered via multiple routes of administrations and has a wide variation in dosing regimens including low (< 50 mg/m 2 ), intermediate (50-500 mg/m 2 ) and high (> 500 mg/m 2 ) dose regimens [1,4]. The pronounced inter-individual variability (IIV) of PK and toxicity of MTX [5][6][7] renders individualization of dosing regimens difficult.
Hepatic metabolism accounts for a considerably lower fraction of its clearance (CL) compared to renal elimination, as the main fraction (80-90%) of the drug is primarily eliminated via glomerular filtration and active tubular secretion [8,9]. Nephrotoxicity associated with MTX impairs its CL, leading to further aggravation of toxicity such as myelosuppression and mucositis. In subjects with extracellular fluid accumulations, the drug has been shown to undergo delayed elimination [10]. A recent in vitro study by Euteneuer et al. [11] showed a sex-dependent regulation of renal transport proteins, which might play a role in the CL of MTX. To handle the variability associated with MTX exposure, monitoring of its plasma concentrations (therapeutic drug monitoring, TDM) and serum creatinine (SCr) is recommended to safeguard a relatively constant drug exposure with an acceptable risk/benefit ratio particularly in patients with impaired renal function [12]. Furthermore, MTX dosing is often guided by body surface area (BSA) estimates to account for body size-related differences in CL and volume of distribution (V). However, concerns regarding potential under-and over-exposure in certain patient groups, such as with obesity, have been expressed [13]. BSA is furthermore a highly variable measure that depends on the arbitrary choice of a BSA equation [14]. Thus, further clarification of the clinical implications of BSA based dosing for MTX is required.
Modeling of PK data has the potential to optimize TDM, where tailored dose adjustments can be made according to model-predicted concentrations of a drug [15]. Bayesian population PK analysis has been used to assist TDM guided dose adjustments for MTX [15]. In addition, population PK analysis provides the possibility to identify and quantify covariate effects on drug exposure [16,17]. This may provide a better understanding of drug's pharmacology and assist adjustments in dosage regimen according to patient's individual characteristics e.g., renal/hepatic function, genotype of drug metabolizing enzymes or transporters, and/or anthropometric characteristics. Models capturing covariate relationships have been found useful in oncology for individualized dose adaptations such as in case of busulfan, topotecan and docetaxel [16].
The current study was aimed to identify and evaluate covariates influencing PK of MTX, particularly patient sex and body surface area (BSA), by developing a population PK model using the TDM data collected from patients with hematological and solid malignancies. The model was further aimed to be used for the evaluation of the ongoing clinical practice of administering MTX based on individual BSA via a simulation study.

Methods
Patients, treatment and sampling MTX plasma concentration and covariate data was obtained from the Cologne Cohort of Neutropenic Patients (CoCoNut) [18]. Experimental protocols were approved by the local ethics committee (name and email address: Ethics Committee of the Faculty of Medicine, University of Cologne, Cologne, Germany, ek-med@uni-koeln.de; date of approval: 14.01.2014, approval file number: 13-108). All methods were performed in accordance with the local and international guidelines and regulations. Data from neutropenic patients (neutrophils < 500 /mm 3 ) with hematological malignancies or solid tumors and treated with high-dose MTX at the Department I of Internal Medicine, University Hospital of Cologne, between January 2005 and February 2018 were considered. The data from clinical laboratory was imported via Health Level Seven from the laboratory information system. The dosing information was imported from the integrated software for chemotherapy using a csv export. Further patient characteristics were documented manually in the CoCoNut database.
MTX was administered via 4 h or 24 h intravenous infusions depending on underlying malignancy. TDM was routinely performed at 42 h and 48 h post-dose for both the 4 h and 24 h protocols, while an additional sample was scheduled for 4 h MTX infusion at 24 h. If target plasma concentration exceeded the desired thresholds (> 1 μmol/L at 42 h and > 0.3 μmol/L at 48 h), TDM was performed at least every 6 h. These thresholds reflect the internal guidance document developed to translate the available heterogenous evidence [10,19] to an actionable recommendation also appropriate for the organisational conditions in our hospital. On the same basis, for the 24 h MTX infusion, leucovorin rescue was routinely performed with 30 mg/m 2 (after 42 h and 48 h) and 15 mg/ m 2 (after 54 h and 60 h). If the desired plasma concentration of MTX was not reached, leucovorin was administered every 6 h at a dose (mg) equivalent to the product of MTX plasma concentration (μmol/L) and body patient weight (kg).
MTX plasma concentrations were quantified using competitive immunoassays with 0.009 μmol/L as the lower limit of quantification (LLOQ). Demographic covariates included patient's age, sex, weight and height. Covariate data from clinical chemistry analysis included SCr, plasma total bilirubin (BT), γ-glutamyltransferase (GGT), uric acid concentrations, absolute leukocyte counts (WBC), and BSA.
Dosing, concentration and covariate data was subjected to screening prior to PK analysis. R (version 3.5.1) was used to prepare the dataset for model development. Dataset preparation was assisted by visual inspection of individual concentration time profiles. Patients with missing dosing information at treatment initiation were identified for exclusion from subsequent analysis. Subjects with missing dosing information during the treatment were flagged and concentration measurements at time points subsequent to the missing dosing information were excluded. Due to the significant amount of missing covariate data throughout the treatment course, the covariate evaluation was based on baseline covariate data for the start of treatment.

PK model development
Data were analyzed by the nonlinear mixed effects modeling approach using NONMEM 7.4.3 (ICON, Development Solutions, Elliot City, MD, USA). Perl speaks NONMEM (PsN), Pirana and Xpose4 were used to assist model development, evaluation and post processing [20][21][22]. Structural model development. A combination of iterative two-stage (ITS) and first order conditional estimation with interaction (FOCE-I) methods was applied for parameter estimation. Likelihood ratio tests (LRT) or the Akaike information criterion (AIC) were used for the evaluation of nested and non-nested models, respectively. A nested model with fewer parameters or a decrease in objective function value (OFV) by 3.84 (i.e., p < 0.05, one degree of freedom) was given preference. The model with a lower AIC value in case of non-nested models was preferred.
Model evaluation criteria comprised of plausibility of parameter estimates, reduction in unexplained and residual variability, shrinkage and precision in parameter estimates. Visual inspection through goodness of fit (GOF) plots included observed versus individual/population predicted concentrations (IPRED/PRED) over time. Residual error models were evaluated with the help of conditional weighted residuals (CWRES) versus observed concentrations and versus time after first dose (TAFD). Numerical predictive checks (NPCs) were used for further assessment by comparing the empirical cumulative distribution function of the observed concentrations with the theoretical cumulative distribution, computed from simulated data.
Compartmental analysis was performed in a stepwise manner. IIV was incorporated using exponential terms (η iiv ) which describes the deviation of PK parameter values of an individual from the population estimate [17]. Interoccasion variability (IOV), defined as the variability between individual cycles of MTX therapy, was incorporated in the model via random effects (η iov ) [23]. The PK parameter P in a specific subject was parametrized as shown in Eq. 1.
Where θ is a fixed effect, representing the median PK parameter in the population. Additive, proportional and combined error models were tested to estimate the residual unexplained variability (RUV).

Covariate model development
Covariate data was analyzed to identify covariateparameter relationships. Covariate preselection was performed considering scientific plausibility as an essential criterion. Graphical evaluation of covariates was performed including CWRES vs covariate, empirical Bayes estimates (EBEs) versus covariate, and covariate versus covariate plots. Significance of covariate relationship was principally guided by decrement in OFV and/or unexplained variability. A stepwise covariate evaluation was carried out as follows. At each step, the covariate providing the largest reduction in OFV was included (forward inclusion) or the covariate providing the lowest increase in OFV was eliminated (backward elimination). Selection criteria were a ΔOFV of 3.84 (p < 0.05) for forward inclusion and a ΔOFV of 6.63 (p < 0.01) for backward elimination.
Continuous covariates were included as linear relationships (Eq. 2) or power relationships (Eq. 3) centered around their median values. BSA effect was centered around the typical value of 1.73 m 2 .
Categorical relationships were given as Covariate effect = 1 + Covariate i × θ Covariate , where Covariate i is the individual covariate value in the i th subject and θ Covariate represents the effect size of the covariate relationship to a PK parameter. Covariate inclusion and evaluation criteria are presented in the supplementary material.

Patient and treatment characteristics
In total, 229 cancer patients (83 females) with 2182 plasma concentration measurements were included in the PK analysis. The majority of patients received 4 h and 24 h infusions, while 18 patients occasionally received 12 h and 48 h infusions. Only a single patient received a 72 h infusion. A median of 3 dosing cycles (range, 1-9) per patient were part of the available data. The number of plasma concentration measurements per patient ranged from 1 to 65 with a median of 7 measurements. Patient and clinical laboratory parameters are summarized in Tables 1 and 2.    Table. RUV was appropriately described by a combined (additive and exponential) error model. Mean PK parameters with 95% CI and RSE obtained from the bootstrap analysis (1000 samples) are presented in Table 3.

Covariate analysis
SCr was found to be a significant covariate on CL with an OFV reduction by 191. Inclusion of patient's sex and age on CL further improved the model fit (ΔOFVs of 32.0 and 13.0, respectively). Inclusion of BSA provided a significant reduction in OFV by 4 Where, sex was coded as 0 for males and 1 for females. Estimates for covariate relationships are summarized in Table 3. The percentage of subjects attaining both the target criteria (probability of target attainment; PTA) was calculated for dose levels of 500, 1000 and 2000 mg/m 2 (reference dose). An optimized flat dosing regimen, i.e. a regimen in which each subject received the same dose, was identified by simulating a range of doses and choosing the dose that provided the highest PTA. Subsequently, the procedure was repeated with doses stratified according to the BSA groups (lower extreme: < 25%, middle proportion: 25-75% and upper extreme > 75%) and the above-mentioned dose optimization was repeated for each of the three BSA regions separately. Thus, the stratified dosing approach resulted in three separate doses, corresponding to the three defined BSA groups. Figure 3 presents the PTA across the BSA groups for respective dose levels under BSA-based, flat and stratified dosing regimens. Stratified dosing provided marginally higher PTA for both the upper and lower BSA extremes compared to BSA-based and flat dosing, respectively. Based on simulation results, selection of doses with the highest PTA (comparable to BSA-based doing regimen) identified under flat and stratified dosing regimens are presented in Table 4.

Discussion
High-dose MTX is essential in cancer therapies despite its high toxicity. However, the management of delayed MTX elimination challenges clinicians to prevent potentially life-threatening MTX-associated toxicities. Further, the high toxicity can cause a premature termination of the MTX administration, which decreases its potential efficiency [10]. In this study, we investigated the optimization of MTX dose adjustment as a potential factor to reduce MTX toxicity. A three-compartment PK model of MTX is presented. Patient sex, age, BSA and SCr were related to CL. A 16% lower CL was estimated for females compared to males.
Simulations using the final covariate model support dosing stratified for BSA quartiles.
The identification of clinically relevant covariates has been the main objective of numerous population PK evaluations of MTX, providing inconsistent findings on covariate relationships [24][25][26][27][28][29][30][31][32][33][34] . In contrast to our study, several previously published studies did not support a sex effect [24,[30][31][32][33][34]. Apart from differences in sample  Table 4 sizes, the particular combination of covariates in the model might have contributed to this inconsistency. For example, the inclusion of sex in a model containing SCr effect provided a distinct model improvement (OFV reduced by 32.0 points). In comparison, only a marginal improvement (OFV reduced by 4.40 points) resulted in the univariate evaluation (i.e., without considering any other covariates). This finding supports that a sex effect should be considered to account for differences in creatinine generation rates in male and female subjects. To quantify the contribution of sex effects beyond renal function, data on urinary excretion might be useful. However, such data was not available for patients in our database. The effect of sex on CL needs further investigation.
Age was related to MTX CL in a few studies [30,32] while inconsistencies exist in the majority of studies [25,[34][35][36][37][38]. Some studies presented the influence of body weight and patient's age on both the CL and V of MTX [31,33]. Mei et al. showed that V of MTX increased with increasing age and supported the preference of age over body weight as a covariate influencing V. A relationship between weight and V was reported by some other studies as well [30,31,34,39,40]. Age was found to be significant on CL in our study with a ΔOFV of 13.0.
SCr was the most significant covariate with a ΔOFV of 191. This is in line with other studies, where MTX elimination was significantly correlated with SCr [34,41,42]. The observed effect is physiologically plausible as MTX is primarily eliminated by the kidney [24]. Nevertheless, the covariate relationship between SCr concentrations and MTX CL faces disagreements in other studies [30,31,33,37].
The covariate analysis in the present evaluation was based on baseline covariate information due to missing covariate data during the treatment time course for a significant number of patients. The development of covariate models incorporating time-varying covariate data is a useful approach in general, as it may provide a better explanation of IIV and IOV of PK parameters and thereby improve the predictive ability of the model. This  Preference of BSA-based dosing over flat dosing or based on other measures, such as patient genotype / phenotype, is an ongoing debate. In contrast to the BSA-based dosing, flat dosing is proposed for several anticancer drugs where BSA has been shown not to reduce the random PK variability to a clinically relevant degree [43][44][45]. Apart from the simplified clinical handling of flat dosing, BSA-based dosing introduces additional uncertainties which are difficult to assess due to the arbitrary choice of BSA equation [46]. Therefore, BSA as a body size measure should ideally be avoided if precise dose calculations are intended. Furthermore, scaling doses with BSA is likely to provide implausibly low or high doses in subjects with exceptionally low or high BSA. Owing to the simplicity needed to implement body sizebased dosing regimens in clinical practice, a direct, proportional relationship between BSA and dose is often assumed. This is contradictory to our current knowledge on physiology and PK and further adds to the uncertainties. Although MTX PK is linear, i.e. exhibits a proportional increase in exposure (in terms of AUC) with the increase in dose, it does not imply that MTX exposure is reciprocally proportional to BSA. In our study, the change of exposure attributable to BSA was smaller and we only observed relevant differences between BSA and flat dosing for patients with either very low or very high BSA. A stratified approach is a reasonable alternative to BSAbased dosing with individuals in the upper and lower BSA quartiles. A stratified approach is a reasonable alternative to BSA-based dosing with individuals in the upper and lower BSA quartiles. It is important to mention that the current findings are based on retrospective data and need to be further validated in a prospective study. It should be noted that these findings are conditional on the defined TDM target. The TDM target is a concentration threshold associated with overexposure, while no threshold for underexposure is currently available. To avoid possible underexposure, achieving 70-130% of the AUC for a subject with a typical BSA of 1.73m 2 was used as an additional criterion in the simulation analysis. No pharmacokinetic/pharmacodynamic (PK/PD) target related to efficacy is part of the TDM at the University Hospital Cologne, and, to the best of our knowledge, no validated PK/PD target related to efficacy is currently available.
Apart from the covariate and BSA evaluation, a nonlinear CL component was identified in this study. Nonlinear elimination has been reported before and might be attributable to the transporter-mediated tubular secretion of MTX [47][48][49]. Despite the significant improvement of the model after inclusion of non-linear CL, the impact of the non-linear component on estimated exposure and the excess of TDM thresholds was negligible. Thus, non-linearity seems to be of minor clinical relevance in the current cohort of patients. This might change if additional targets, such as PK/PD targets related to efficacy, become available. In this case, the model with the non-linear component as presented in the Supplement might be re-evaluated. Furthermore, the non-linear CL component might have a more pronounced impact on PK in presence of genetic polymorphisms and when MTX is co-administered with substrates, inducers, or inhibitors of the relevant membrane transporters.
A major limitation of the current investigation is that the optimal target exposure regarding the efficacy of MTX in various malignancies is unknown. Data on minimum drug exposure needed to achieve a positive therapeutic outcome with minimal toxicity is currently scarce. Dedicated efforts are needed to draw conclusions based on the efficacy profile of the drug with respect to the underlying disease. Furthermore, TDM data is generally obtained from clinical practice and therefore provides a reduced data quality compared to clinical trial data. Although the data was checked carefully for inconsistencies, it cannot be precluded that errors in TDM procedures translate into model misspecifications.

Conclusions
A three-compartment model described PK of MTX. A lower CL estimated for the female patients needs to be investigated in future studies. Plasma SCr, patient age, sex and BSA were found additionally as statistically significant covariates on CL. Stratified MTX dosing can be a reasonable alternative to BSA guided dosing.