A Limited Sampling Schedule to Estimate Individual Pharmacokinetic Parameters of Fludarabine in Hematopoietic Cell Transplant Patients

Allogeneic hematopoietic cell transplants (HCT) have curative potential for patients with hematologic diseases and malignancies (1, 2). Myeloablative conditioning regimens, however, have high rates of nonrelapse mortality, thus limiting their availability to older or medically infirm patients (1–3). These patients can now receive a HCT with a reduced-intensity conditioning regimen (4, 5). Administered i.v. as a more soluble prodrug fludarabine monophosphate (Fludara, abbreviated with the term “fludarabine” for the rest of the article), the purine nucleoside analogue fludarabine is administered at doses ranging from 90 to 250 mg/m² in reduced-intensity HCT conditioning regimens (2). Fludarabine administration depletes lymphocytes and improves engraftment rates [versus total body irradiation (TBI) alone], but it rarely causes fatal nonhematologic toxicity (2, 5, 6). Identification of the optimal fludarabine dose, one that achieves the delicate balance between engraftment rates and toxicity, would be of significant benefit to patients receiving nonmyeloablative HCT.

After administration, fludarabine is quickly dephosphorylated to 2-fluoroadenine (F-ara-A) by nucleotidases (7–9). F-ara-A subsequently undergoes uptake and sequential phosphorylation to its active metabolite, fludarabine triphosphate (F-ara-ATP), which inhibits ribonucleotide reductase and DNA polymerase (8). The toxicity to fludarabine monotherapy is not predicted by body surface area (BSA)–based dosing (10); however, an inverse relationship has been observed between plasma F-ara-A area under the concentration-time curve (AUC) and neutropenia (11). Preliminary data from patients receiving fludarabine as part of myeloablative HCTs suggest conditioning regimen-dependent pharmacodynamic relationships with F-ara-A AUC (12, 13). Understanding the sources of variability in nonmyeloablative HCT patients is of the utmost importance because the success of HCT is based on achieving a delicate balance between the donor's and recipient's immune systems. The pharmacokinetic variability in F-ara-A may contribute to this narrow balance between the donor and recipient; however, this hypothesis has yet to be tested. Because fludarabine is administered in an ambulatory care clinic, it is critical to develop a limited sampling schedule (LSS) to facilitate patient compliance with the pharmacokinetic sampling necessary to evaluate the pharmacodynamics of F-ara-A in nonmyeloablative HCT recipients. To this end, we sought to develop a population pharmacokinetic model and LSS to minimize the patient burden for participation in pharmacodynamic studies to evaluate if F-ara-A exposure predicts outcomes. Future pharmacodynamic studies are essential for determining the relevance of the F-ara-A pharmacokinetic variability on clinical outcomes in outpatient nonmyeloablative HCT recipients.

Forty-two patients received myeloablative conditioning with busulfan and fludarabine. The busulfan doses were personalized to achieve a target (^TBU) concentration at steady state of 800 to 1,000 ng/mL. Fludarabine doses were not changed based on F-ara-A pharmacokinetic data. F-ara-A pharmacokinetic data were available on 27 patients receiving targeted oral busulfan for 4 d after fludarabine (30 mg/m²/d for 4 d; protocol 1519; ref. 13). F-ara-A pharmacokinetic data were available in 15 additional patients enrolled on a subsequent ^TBU/fludarabine protocol (protocol 2041). This regimen consisted of 50 mg/m²/d fludarabine for 5 d (days −6 to −2), targeted daily i.v. busulfan on days −5 to −2, and thymoglobulin on days −3 to −1. In both regimens, fludarabine was infused through a central venous access catheter over 30 min. Institutional Review Board approval was obtained before study conduct, and all patients or patient guardians gave informed consent before their participation in the study began.

All patients had been diagnosed with hematologic disorders and had adequate renal (i.e., serum creatinine <1.5 mg/dL and creatinine clearance or radioisotope glomerular filtration rate >60 mL/min/1.73 m²) and liver (i.e., total bilirubin <1.5 mg/dL and alanine aminotransferase <300 units/L) function. The body weight and BSA used to calculate doses were based on institutional guidelines.

We evaluated demographic data, concomitant medications, and relevant laboratory tests for 3 d before and during fludarabine administration to identify potential covariates. Potential covariates included age at the time of fludarabine administration, gender, BSA calculated using actual body weight, adjusted ideal body weight, height, serum creatinine, creatinine clearance (estimated using Cockcroft-Gault formula; ref. 14), blood urea nitrogen, and albumin. Patients received similar antiemetics, antibiotics, and antifungals per institutional Standard Practice Guidelines; corticosteroids were not used as antiemetics with fludarabine. Genetic covariates were not evaluated because fludarabine is mainly renally eliminated (11, 15, 16).

F-ara-A pharmacokinetic blood samples were drawn from a catheter not used to infuse the fludarabine dose. In the 27 protocol 1519 patients who received 30 mg/m²/d fludarabine, pharmacokinetic blood samples were scheduled to be drawn immediately before and at the end of the 30-min infusion (0.5 h) and 1, 4, 8, 12, and 24 h after the start of the infusion. The 24-h pharmacokinetic sample was drawn before the next dose. These samples were collected after all four fludarabine doses. The F-ara-A pharmacokinetic data from this study were previously reported (13). Outpatient fludarabine administration necessitated changing the pharmacokinetic sampling times in the subsequent study (i.e., protocol 2041). In the 15 protocol 2041 patients who received 50 mg/m²/d fludarabine, pharmacokinetic samples were drawn at 0.5 h and 5 min after the end of the infusion (0.583 h) and 1.5, 4.5, 6.5, and 24 h after the start of the infusion. The 24-h pharmacokinetic sample was drawn before the next dose. Pharmacokinetic samples were collected with the first fludarabine dose only (i.e., day −6) in protocol 2041 patients.

F-ara-A plasma concentrations were quantitated using liquid chromatography-mass spectrometry as previously described (13); the dynamic range was 0.23 to 9.04 μmol/L and the interday coefficient of variation (CV) was <10%.

Seven hundred twenty-eight concentration-time points were available for population pharmacokinetic modeling. Concomitant busulfan does not interact with F-ara-A (17), and therefore, the data sets were combined for the population pharmacokinetic modeling and construction of the LSS. The pharmacokinetic model of F-ara-A was fit using the first-order method in the nonlinear mixed-effects modeling software NONMEM (version VI; ref. 18). A two-compartment model of F-ara-A was initially tested to estimate the following parameters: CL, the apparent elimination clearance; V₁, the apparent volume of distribution in plasma; Q, the intercompartmental clearance; and V₂, the apparent volume of the peripheral compartment. Between-subject variability (BSV) of parameters was initially modeled using a full covariance matrix; other alternative covariance matrix structures were tested as well. We assumed a log-normal distribution of BSV for all parameters. Residual unknown variability was estimated with a combination of additive and proportional error models. Individual parameters were calculated as empirical Bayes estimates using the POSTHOC option in NONMEM. Criteria for model selection included visual inspection of goodness-of-fit plots (data versus population and individual predictions), residuals, and weighted residuals. Plots of data versus population and individual predictions, data versus residuals, and data versus weighted residuals were visually inspected for goodness of fit. Models with different numbers of parameters were selected based on expected differences in χ² distributions with appropriate degrees of freedom (19).

Subsequently, we evaluated if the model was improved by inclusion of the following covariates in the base pharmacokinetic model: age, gender, BSA, weight (actual or adjusted ideal body weight), height, creatinine, creatinine clearance, blood urea nitrogen, and albumin. Linear regression relating individual parameter values and covariates as well as visual inspection of plots were used to explore these relationships. Inclusion/exclusion of covariates was determined by parsimony criteria.

As a validation of the model, a visual predictive check was done (20). We simulated data for 1,000 subjects using the final population pharmacokinetic model with BSA-normalized dosing (see Results for a rationale for BSA-normalized dosing). The median and 90% prediction intervals were plotted with the original data.

We sought to choose an optimal LSS that balanced the practical challenges of pharmacokinetic sampling in patients with the need for sufficient data to predict an individual's F-ara-A AUC. Specific challenges were that the LSS had to have few enough samples to encourage patient compliance (chosen as a maximum of four pharmacokinetic samples) and a short enough time span that patients can receive fludarabine in an ambulatory clinic (i.e., within 6.5 h from the start of the infusion). Therefore, three of the four LSS times were chosen within the first 6.5 h after the start of the fludarabine infusion. The fourth LSS was set at 24 h, which is immediately before the next fludarabine dose. Individual pharmacokinetic parameters were determined via maximum a posteriori (MAP) Bayesian estimation, where the population prior (from the population pharmacokinetic analysis) may potentially offset the loss of individual data due to the use of the LSS.

We hierarchically simulated 2,000 subjects whose pharmacokinetic parameters were based on our population pharmacokinetic model of F-ara-A after i.v. fludarabine administration. Each subject arose from the model BSV and was replicated 20 times using different simulated residual unknown variability (data error) instantiations. For each potential LSS, we did individual MAP parameter estimation for each subject's simulated data for each error replicate and computed both the 24-h and 120-h AUCs after one simulated fludarabine infusion. The modeled AUCs were compared with the true (simulated) AUCs to compute scaled mean squared error (sMSE). The mean (across subjects) sMSE was the figure of merit for each LSS. The LSS with minimum sMSE was considered optimal. For each simulated subject,

were computed for the AUC for each potential time grid, where $θ̂i$ is the AUC for error replicate i, and $θ̂$ is the “true” AUC value.

To determine how the LSS would perform in the calculation of the AUC compared with the intensive, full pharmacokinetic sampling (FULL, at 0.5, 0.583, 1.5, 4.5, 6.5, and 24 h after dose), where all six concentration-time data points were fit to the model via individual maximum likelihood estimation (MLE; i.e., without parameter prior), the same subjects were also considered and the root mean sMSE was computed for FULL.

The framework to accomplish this LSS optimization was constructed using the R statistical software, with the needed simulation and parameter estimation done via system calls to NONMEM.

Of the 42 patients studied, 21 were male and 21 were female. At the time of fludarabine administration, the median age was 51.9 years (range, 12.6-65.5), the median BSA (based on actual body weight) was 1.97 m² (range, 1.56-2.60), and the actual body weight was 79.0 kg (range, 55.4-139.7). Creatinine clearance (n = 40) was 89 mL/min (range, 55-148).

We determined that the best population pharmacokinetic model for these data was a first-order pharmacokinetic, two-compartment model with a 30-minute fludarabine infusion as an input to the central compartment. Comparison of a full covariance matrix BSV to various other configurations indicated that the correlations between CL (apparent elimination clearance) and Q (the intercompartmental clearance) and CL and V₂ (apparent volume of the peripheral compartment) were not significantly (P < 0.05) different from zero. All covariance configurations confirmed a very strong correlation (near 1) between Q and V₂. Estimation of a correlation near 1 can lead to numerical instabilities and may disrupt determination of convergence. Thus, for the final base model, we concluded that the best BSV covariance structure was a banded covariance matrix [implying that corr(CL,Q) = corr(CL,V₂) = 0], with the additional assumption that the correlation of Q and V₂ was fixed to 1. In addition, by this construction, corr(V₁,V₂) = corr(V₁,Q). More detail about this covariance structure and its formulation in NONMEM is provided in  Appendix 1. A combined additive and proportional residual unknown variability error model performed better (P < 0.05) than either an additive or a proportional error model alone.

Visual inspection of plots and linear regression relating individual parameter and covariate values revealed the strongest association between all parameters (i.e., CL, V₁, Q, and V₂) and the various measures of body size (i.e., BSA, height, actual body weight, and adjusted ideal body weight). Because the body size measures (except height) were highly correlated, we chose BSA because fludarabine is dosed based on BSA (e.g., 30 mg/m²/d). As a first step, we evaluated BSA-normalized dosing, essentially testing BSA as a covariate (scaling) of all parameters. This inclusion was justified by a 13.3-point improvement in objective function value with no additional parameters and a 4% point decrease in BSV for each parameter. The structural parameters fixed effects for this BSA-normalized dosing are provided in Table 1. The right half of Table 1 presents the BSV of the model parameters as the lower half of a matrix, with variability (%CV) on the diagonal and parameter correlation (as Pearson r) in the off-diagonal elements.

Using the BSA-normalized dosing model, visual inspection of plots and linear regression relating individual parameter and covariate values revealed weak relationships between V₁ and gender and both CL and V₁ and height. These were added separately as linear parameter-covariate relationships to the BSA-normalized model. Inclusion of V₁ as a linear model in terms of height and (separately) of gender did not reach statistical significance (P < 0.05; objective function improvement of <3.84; ref. 19).

The inclusion of CL as a linear model of height did reach the required P < 0.05 level of statistically significant objective function improvement (the decrease was 5.5), although it was not strongly significant (i.e., not at the P < 0.01 level, which would require an objective function decrease of 6.6). The typical CL value of 5.70 L/h/m² (Table 1) was replaced by 5.74 + 0.0334 * (height − 170.5), where height is the subject's height and the median height was 170.5 cm. This inclusion also resulted in the explanation of 4.2% points of BSV in CL. Other parameter and variability values changed only negligibly.

Figures 1 and 2 provide model validation in the form of a visual predictive check. The median and 90% prediction intervals were plotted with the original data for subjects receiving i.v. 30 mg/m²/d fludarabine with pharmacokinetic sampling after all four doses (Fig. 1) and i.v. 50 mg/m²/d fludarabine with pharmacokinetic sampling after the first dose only (Fig. 2). The visual predictive check shows generally good agreement between model and data, although somewhat fewer than 10% of data points lie outside the 90% prediction interval.

LSS results are presented in Table 2 ranked by root mean sMSE for each potential LSS (so that the units are in percent SD). A comparison of schedules based on root mean sMSE for the 24- and 120-hour AUC showed several schedules performing similarly well (Tables 2 and 3). We chose the 0.583-, 1.5-, 6.5-, and 24-hour LSS as the optimal schedule. This decision was based on the low sMSE of the schedule and on our experience with patient compliance with pharmacokinetic sampling. Other potential schedules performed nearly as well. The performance of the FULL schedule (six observations; model fit by MLE parameter estimation, i.e., no population prior) was nearly as good as the chosen LSS for determining the 24-hour AUC but performed worse than the LSS for determining the 120-hour AUC.

This work represents a critical step for future F-ara-A pharmacodynamic studies in patients receiving fludarabine in the ambulatory clinic, in which the pharmacokinetic scheduling must be convenient enough to ensure patient compliance while still accurately characterizing an individual's pharmacokinetic disposition of F-ara-A. To this goal, we developed a novel population pharmacokinetic model and determined a LSS. This minimally intrusive LSS, when used in conjunction with population prior values, performs as well as (Table 2) or better than (Table 3) more intense pharmacokinetic sampling with the individual's pharmacokinetic data alone using compartmental analysis (MLE). The LSS chosen to characterize an individual's F-ara-A pharmacokinetics, when used in combination with population pharmacokinetic priors, was 0.583, 1.5, 6.5, and 24 hours after the start of the 0.5-hour fludarabine infusion.

In a subsequent investigation, we are using this LSS to test the hypothesis that plasma F-ara-A AUC is associated with day 28 T-cell donor chimerism in patients receiving fludarabine/TBI conditioning in an outpatient setting. Day 28 T-cell chimerism is associated with engraftment, relapse, and acute graft-versus-host disease rates (2, 5, 6), and therefore, we seek to determine if F-ara-A AUC is associated with day 28 chimerism. As fludarabine/TBI is one of the most common reduced-intensity regimens, it is imperative to identify biomarkers associated with donor chimerism in these nonmyeloablative HCT recipients (21). The availability of HCT has expanded due to the recent development of reduced-intensity conditioning regimens (2, 22). With the current method dosing of fludarabine (i.e., based on BSA), there is substantial variability in donor chimerism, efficacy, and toxicity in nonmyeloablative HCT recipients. We hypothesize that this interpatient variability in donor chimerism may be related to interpatient variability in F-ara-A pharmacokinetics. We are the first to develop a population pharmacokinetic model of F-ara-A in HCT recipients and determine a LSS that requires only four blood samples to estimate an individual's F-ara-A AUC when used in combination with population parameter priors. With the development of these critical tools, we are testing the hypothesis that F-ara-A AUC is associated with clinical outcomes in patients receiving nonmyeloablative HCT conditioning in the ambulatory clinic.

The interpatient variability in F-ara-A clearance is comparable with that of busulfan, the dosing of which is personalized based on pharmacokinetics at many HCT centers (23). It is critical to establish if the pharmacokinetic variability in F-ara-A AUC is associated with clinical outcomes and, if so, subsequently conduct a clinical trial to show the clinical benefit of personalized dosing. Our group has taken that approach with both busulfan (24) and cyclophosphamide (25, 26). With both these alkylating agents, the pharmacodynamics differs based on the conditioning regimens. Preliminary data suggest conditioning regimen-dependent pharmacodynamic relationships with F-ara-A AUC. In patients receiving busulfan/fludarabine myeloablative conditioning, F-ara-A AUC was not related to T-cell chimerism or engraftment (13). However, in patients receiving fludarabine in combination with cyclophosphamide and TBI, those with a low F-ara-A AUC were less likely to experience neutrophil engraftment, whereas the patient with the highest F-ara-A AUC died due to neurotoxicity (12). Therefore, future studies testing the hypothesis that variable clinical outcomes are associated with F-ara-A AUC should be conducted in a homogenous group of nonmyeloablative recipients, all of which are receiving fludarabine/TBI. With fludarabine-containing conditioning regimens, the development of this population pharmacokinetic model and LSS provides critical tools to evaluate if F-ara-A AUC is associated with clinical outcomes (2, 21, 27–29). The LSS requiring four blood samples is minimally intrusive to HCT patients treated in the ambulatory clinic while still allowing for accurate characterization of F-ara-A AUC to be used in our ongoing pharmacodynamic study in nonmyeloablative conditioned patients receiving fludarabine/TBI.

In conjunction with characterizing F-ara-A AUC, we seek to also evaluate the role of F-ara-ATP in donor chimerism. After administration, fludarabine is quickly dephosphorylated to F-ara-A by nucleotidases (7–9). F-ara-A subsequently undergoes uptake and is then sequentially phosphorylated to its monophosphate, diphosphate, and triphosphate metabolites (8). The inhibition of ribonucleotide reductase and DNA polymerase by F-ara-ATP, the active moiety, ultimately leads to cellular apoptosis (8). We have recently created a novel method to phenotype F-ara-ATP accumulation in CD4⁺ and CD8⁺ cells (30). F-ara-ATP accumulation, used in conjunction with F-ara-AUC, can be used to estimate F-ara-ATP concentrations in CD4⁺ and CD8⁺ cells.

This population pharmacokinetic model is the first to be constructed for F-ara-A plasma pharmacokinetics. The model estimates of clearance agree with previously reported data analyzed with noncompartmental or compartmental pharmacokinetic analysis in cancer patients not undergoing an allogeneic HCT (11, 15–17). Therefore, with further validation, this population pharmacokinetic model and LSS may be used in nontransplant patients receiving i.v. fludarabine in the ambulatory clinic. Outside of the setting of allogeneic HCT, fludarabine is used to treat patients with chronic lymphocytic leukemia, follicular lymphoma, or acute myeloid leukemia (9, 31, 32). In patients receiving fludarabine monotherapy, an inverse relationship has been observed between plasma F-ara-A AUC and neutropenia (11). It is also of interest to determine if elevated F-ara-A AUC is associated with neurotoxicity outside the setting of allogeneic HCT (12), which is particularly important because neurotoxicity was dose limiting in phase I studies of fludarabine monotherapy (7). Therefore, the development of this population pharmacokinetic model and LSS provides tools for future studies to test the hypothesis that elevated F-ara-A AUC is related to toxicity outside the setting of allogeneic HCT (12).

Population pharmacokinetic analysis provides a rigorously quantitative approach to incorporating population information (i.e., typical values and variability of the disposition of drugs through the body) in the determination of parameter priors, which can assist in estimation of an individual's exposure (via MAP Bayesian estimation). The benefits of population pharmacokinetic modeling are just starting to be realized in HCT recipients, as evidenced by our results with cyclophosphamide (33, 34). This population pharmacokinetic model is the first, to our knowledge, to be built for F-ara-A. It was built using F-ara-A concentration-time data on days that only fludarabine, not busulfan or thymoglobulin, was administered to patients conditioned with myeloablative doses of ^TBU/FLU. Notably, concomitant busulfan does not interact with F-ara-A pharmacokinetics (17). A limitation of this F-ara-A data set was that all patients had adequate renal and liver function to receive myeloablative ^TBU/FLU. However, some nonmyeloablative HCT recipients may have impaired renal or liver function. Notably, none of the indices of renal or liver function was covariates in our analysis (see “Population pharmacokinetic modeling” of Results). An improved understanding of the covariates (e.g., age) of BSV can lead to improved experimental design and determination of a LSS (35). BSA was the only covariate included in the model. The inclusion of subject height as a covariate for apparent clearance (CL) was found to barely reach statistical significance (P < 0.05); our opinion is that inclusion of this relationship in the model will likely not have clinical effect. We observed a very strong correlation of parameters Q and V₂ in our model. To account for this correlation, we tried alternative covariance matrix structures, which did allow us to fix a correlation to <1, but these required a block-diagonal covariance matrix that assumed a zero correlation between parameters V₁ and both Q and V₂. These alternative covariance structures were determined (by parsimony criteria) to be inferior to our final model. In addition, fixing corr(Q,V₂) to, say, 0.98 seemed as arbitrary as fixing the correlation to 1. Therefore, although it is unlikely that two model parameters will have an actual correlation of 1 among subjects, we felt this assumption was justified. Pertinent to the creation of the LSS, although D-optimality is the customary optimality criterion in population experimental design, we opted to use MSE in this individual experimental design setting. D-optimality, the determinant of the Fisher information matrix, accounts for parameter precision only and ignores estimator bias. We expect our individual MAP estimates to be biased due to the presence of the population prior (36). Using MSE accounts for both estimator precision and bias (37).

Our successful development of the first population pharmacokinetic model and LSS to characterize F-ara-A AUC represents an essential step in executing future pharmacodynamic studies in patients receiving fludarabine in an ambulatory clinic. This minimally intrusive LSS will allow us to test the hypothesis that F-ara-A AUC, estimated through the use of this four-sample LSS and population pharmacokinetic priors, is associated with donor chimerism and clinical outcomes in a homogenous population of patients receiving nonmyeloablative conditioning of i.v. fludarabine/TBI.

The base pharmacokinetic model was implemented in NONMEM. Implementation of the assumed corr(Q,V₂) = 1 is reported below, as well as some notes of explanation.

The model parameterizations and log-normal BSV are formulated as follows:

Notes:

• Cov(CL,V₁) = θ₁θ₁Cov(η₁,η₂)

• Corr(CL,V₁) = Corr(η₁,η₂)

By construction, the correlation corr(Q,V₂) = 1:

θ₅ can be interpreted as the (square root of the) scale parameter between Var(Q) and Var(V₂). The value of θ₅ corresponding to the parameter values in Table 1 is 1.04 (11.4% SE).

No potential conflicts of interest were disclosed.

We thank the study participants, Terry Furlong for his assistance with collection of the F-ara-A pharmacokinetic data, Meagan Bemer for her assistance with database management and manuscript preparation, and Rainer Storb, M.D., for his insightful comments to an earlier version of this manuscript.

The base pharmacokinetic model was implemented in NONMEM. Implementation of the assumed corr(Q,V₂) = 1 is reported below, as well as some notes of explanation.

The model parameterizations and log-normal BSV are formulated as follows:

Notes:

• Cov(CL,V₁) = θ₁θ₁Cov(η₁,η₂)

• Corr(CL,V₁) = Corr(η₁,η₂)

By construction, the correlation corr(Q,V₂) = 1:

θ₅ can be interpreted as the (square root of the) scale parameter between Var(Q) and Var(V₂). The value of θ₅ corresponding to the parameter values in Table 1 is 1.04 (11.4% SE).