## Visual Overview

## Abstract

Paclitaxel is a commonly used cytotoxic anticancer drug with potentially life-threatening toxicity at therapeutic doses and high interindividual pharmacokinetic variability. Thus, drug and effect monitoring is indicated to control dose-limiting neutropenia. Joerger et al. (2016) developed a dose individualization algorithm based on a pharmacokinetic (PK)/pharmacodynamic (PD) model describing paclitaxel and neutrophil concentrations. Furthermore, the algorithm was prospectively compared in a clinical trial against standard dosing (Central European Society for Anticancer Drug Research Study of Paclitaxel Therapeutic Drug Monitoring; 365 patients, 720 cycles) but did not substantially improve neutropenia. This might be caused by misspecifications in the PK/PD model underlying the algorithm, especially without consideration of the observed cumulative pattern of neutropenia or the platinum-based combination therapy, both impacting neutropenia. This work aimed to externally evaluate the original PK/PD model for potential misspecifications and to refine the PK/PD model while considering the cumulative neutropenia pattern and the combination therapy. An underprediction was observed for the PK (658 samples), the PK parameters, and these parameters were re-estimated using the original estimates as prior information. Neutrophil concentrations (3274 samples) were overpredicted by the PK/PD model, especially for later treatment cycles when the cumulative pattern aggravated neutropenia. Three different modeling approaches (two from the literature and one newly developed) were investigated. The newly developed model, which implemented the bone marrow hypothesis semiphysiologically, was superior. This model further included an additive effect for toxicity of carboplatin combination therapy. Overall, a physiologically plausible PK/PD model was developed that can be used for dose adaptation simulations and prospective studies to further improve paclitaxel/carboplatin combination therapy.

## Introduction

Cytotoxic drugs used in cancer treatment typically have a narrow therapeutic window, with severe toxicity caused by their unspecific mode of action affecting rapidly dividing cells on one side and the need for sufficiently high doses for efficacy on the other side (Gurney, 1996). Paclitaxel is a cytotoxic agent that is used to treat different cancer types. Paclitaxel enhances and stabilizes the polymerization of microtubules, leading to clinically relevant toxicity, especially dose-limiting neutropenia. This frequent and severe adverse event is caused by the cytotoxic effect on proliferating cells, particularly bone marrow progenitor cells, leading to potentially life-threatening infections (Mitchison, 2012). Neutropenia caused by paclitaxel-containing combination therapies can even be cumulative, meaning a worsening of neutropenia over the repeated treatment cycles (Huizing et al., 1997). In addition, the pharmacokinetics (PK) of paclitaxel administered in a formulation with Cremophor EL is nonlinear and schedule dependent (e.g., influence of infusion duration) and exhibits high interindividual variability. The time of paclitaxel plasma concentration above the threshold of 0.05 *µ*M (T_{C>0.05}) was found to be a relevant PK exposure surrogate for toxicity (Gianni et al., 1995; Ohtsu et al., 1995; Huizing et al., 1997a; Lichtman et al., 2006; Joerger et al., 2007b) and efficacy (Huizing et al., 1997a; Joerger et al., 2007b) (i.e., 26–31 hours; Joerger et al., 2012). The combination of a narrow therapeutic window and high interindividual variability highly favors dose individualization based on therapeutic drug monitoring.

Currently, the dose of paclitaxel is adjusted for body surface area (BSA), which influences paclitaxel PK (Smorenburg et al., 2003), but not for other relevant factors such as organ (dys)function or age (Mielke, 2007). Furthermore, paclitaxel is typically part of a combination therapy, with carboplatin or cisplatin displaying toxic effects on the hematopoietic system as well. In summary, individualized therapy of paclitaxel combination therapy is needed to balance toxicity and efficacy.

PK/pharmacodynamic (PD) modeling and simulation is suggested as a tool to improve the dosing of paclitaxel by reducing toxicity without compromising efficacy, given the complex PK of paclitaxel and the combination therapy (Joerger et al., 2016). Based on a pooled PK/PD analysis from several clinical trials including different cancer types, a PK/PD model to characterize paclitaxel plasma concentrations (PK) and resulting neutropenia (PD) was developed (Joerger et al., 2012). A three-compartment model including nonlinear distribution and elimination was used to describe the PK, while the neutropenia model structure developed by Friberg et al. (2002) (referred to hereafter as the “Friberg model”) was applied for the PD for paclitaxel-associated neutropenia. The PK/PD model was then used to develop a dose individualization algorithm, considering different covariates (BSA, age, sex) as well as target drug exposure (T_{C>0.05}) and toxicity (grade of neutropenia) from the previous cycle. Next, the dosing algorithm was applied in the Central European Society for Anticancer Drug Research Study of Paclitaxel Therapeutic Drug Monitoring (CEPAC-TDM; Joerger et al., 2016). The aim of this prospective, multicenter, randomized controlled clinical study was to investigate whether therapeutic drug monitoring and target-concentration intervention, based on sparse sampling, of paclitaxel was able to reduce toxicity while maintaining efficacy compared with standard BSA-based dosing of paclitaxel. The results showed that the exposure target of T_{C>0.05} (26–31 hours) was not met for more than 50% patients, and no significant improvement of grade 4 neutropenia was achieved for patients in the experimental arm (Joerger et al., 2016). Nevertheless, paclitaxel-related neuropathy was substantially improved. With this clinical trial, an important step was made toward an optimal and individualized dosing of paclitaxel in combination with cisplatin/carboplatin. Reasons for the lack of improvement in paclitaxel-related neutropenia may include suboptimal PK and PD models used for the paclitaxel dosing algorithm as used in the CEPAC-TDM trial. Thus, further improvement in the models and dosing algorithm is needed to reduce the toxicity burden for patients.

This study aimed to investigate the appropriateness of the original PK/PD model for the data from the CEPAC-TDM trial and to detect potential model misspecifications. Furthermore, the PK/PD model should be optimized in order to adequately capture cumulative neutropenia after repeated treatment cycles of paclitaxel in combination with cisplatin/carboplatin. The optimized model will consider physiologic plausibility in order to simulate different dose optimization scenarios of the paclitaxel combination chemotherapy.

## Materials and Methods

#### Clinical Data.

To evaluate the previously developed PK/PD model, we used data from the CEPAC-TDM study (Joerger et al., 2016), which was carried out in accordance with the Declaration of Helsinki and approval from the respective institutional review boards (Ostschweiz, St. Gallen, Eberhard-Karls-Universitaet, Tuebingen). Briefly, patients with newly diagnosed advanced non–small cell lung cancer were treated with paclitaxel (3-hour intravenous infusion) in combination with either cisplatin or carboplatin every 3 weeks for up to six cycles. In the conventional study arm, 182 patients (740 treatment cycles in total) received the standard dose of paclitaxel (200 mg/m^{2}); in the experimental study arm, 183 patients (720 treatment cycles in total) received a paclitaxel dose adapted to sex, age, and BSA for the first treatment cycle. In the following cycles, the dose for these patients was further adapted based on 1) the grade of neutropenia experienced and 2) paclitaxel exposure (expressed as T_{C.>0.05}), both of the previous cycle. PK samples were obtained from patients in the experimental arm only, since only here the dose adaptation based on T_{C.>0.05} was performed. Samples were taken once per cycle, approximately 24 hours (16–30 hours after the start of the infusion). T_{C>0.05} was determined by post hoc estimation with the original PK model in NONMEM software (version 7.3; Icon Development Solutions, Ellicott City, MD). PD samples (neutrophil measurements) were obtained from participants in both arms at baseline, on day 1 and day 15 ± 2 in each cycle, and at the end-of-treatment visit. Paclitaxel concentrations were determined by liquid chromatography and UV detection (lower limit of quantification, 0.015 mg/L, which equals 0.017 *µ*M; recovery, 90.6 ± 9.63; overall precision, <10% CV) (Zufía López et al., 2006). Neutrophil concentrations were measured via routine clinical chemistry analysis at each study site.

The low number of measurements below the lower limit of quantification (PK, 0.30%; PD, 0.40%) as well as the number of missing samples (PK, 8.61%; PD of conventional arm, 11.4%; PD of experimental arm, 9.25%) were assumed to be missing completely at random and were removed from the subsequent PK/PD analysis. For this PK/PD analysis, data from the experimental arm, data from the were primarily considered, since only these patients underwent paclitaxel PK sampling. Only for evaluation of the newly developed PK/PD model data (dosing information, neutrophil concentration measurements and covariate information) from the conventional arm was used. A summary of patient characteristics is available in Joerger et al. (2016) and a summary of the entire modeling and simulation workflow is provided in Supplemental Fig. 1.

#### Original PK Model and External Evaluation.

This analysis was based on a previously published PK/PD model (referred to hereafter as the “original model”) (Joerger et al., 2012), which was also used to develop the CEPAC-TDM dosing algorithm. The PK of paclitaxel was described by a three-compartment model (Fig. 1, upper left) with nonlinear distribution to the first peripheral compartment and nonlinear elimination. BSA, sex, age, and total bilirubin concentration were implemented as covariates on the maximum elimination capacity (VM_{EL}) using power relations. An exponential model was assumed for interindividual, interoccasion, and residual variability.

This PK model was evaluated externally using data from the experimental arm of the CEPAC-TDM study. For this purpose, post hoc estimation with the original PK model was performed to obtain individual PK parameters (empirical Bayesian estimates or EBEs) and the individually predicted paclitaxel concentration-time profiles. Model performance was then evaluated by basic goodness-of-fit plots, shrinkage (Savic and Karlsson, 2009), and comparisons of 1000 simulations with the study data in prediction-corrected visual predictive checks (pcVPCs) (Bergstrand et al., 2011). As in the original model, the first-order conditional estimation method with interaction was used for all PK analyses.

#### PK Model Optimization Using Prior Information.

To refine the PK model for the population in the CEPAC-TDM study, the final PK parameter estimates and their precision (diagonal elements of the variance-covariance matrix) were retrieved from the original PK model and were implemented in a maximum a posteriori (MAP) Bayesian approach (originally referred to as a frequentist approach) as prior information using the normal-inverse Wishart distribution (Gisleskog et al., 2002). Interoccasion variability was not re-estimated but was assumed to be the same as originally estimated, since only one PK measurement per cycle (equal to occasion) was available. We calculated the degrees of freedom for estimating interindividual variability as described in Bauer (2014); for residual variability, 1000 was chosen as the lowest number allowing stable estimation.

The optimized PK parameters were compared with the original parameters at the population and individual levels. A bootstrap analysis (Efron and Tibshirani, 1986) of 1000 bootstrap datasets was performed to evaluate parameter precision, and 95% confidence intervals and population parameter estimates were compared. Confidence intervals of the original model were calculated based on the standard error (retrieved from the original model), assuming normal distribution. To evaluate improvements in model prediction, a pcVPC was generated (*n* = 1000 simulations).

Prediction errors between individual parameter estimates, taking the optimized PK model as a reference and comparing it with the original PK model, were calculated as bias and imprecision according to the following equations (Sheiner and Beal, 1981):(1)(2)(3)in which P_{i,o,original} and P_{i,o,optimized} are the EBEs and exposure parameters (T_{C>0.05} and area under the curve) of individual “i” at occasion “o” (if interoccasion variability was applicable) based on the original and optimized models, respectively. RPE_{p,i,o} is the relative prediction error between the original and optimized parameters for each model parameter “p.” MRPE_{p} is the median relative prediction error indicating bias. MARPE_{p} is the absolute (unsigned) value of the relative prediction error expressing imprecision for each parameter taking the optimized parameters as a reference.

#### Original PK/PD Model and External Evaluation.

The neutrophil concentrations in Joerger et al. (2012) were described by the semimechanistic neutropenia PK/PD model developed by Friberg et al. (2002). For the sequential PK/PD analysis, EBEs of the optimized population PK parameters and the associated individual concentration-time profiles of paclitaxel were used to inform the drug effect (E_{drug}) of the PK/PD model by introducing a linear PK/PD relationship with the slope factor (SL). This method implies a cytostatic drug effect if E_{drug} ≤ 1 and a cytotoxic effect if E_{drug} > 1, since the proliferating rate constant *k*_{prol} was multiplied by 1 − E_{drug}.

In the Friberg model, granulopoiesis was described by a chain of five compartments, the first of which represents the proliferating cells (denoted as “Prol”). These cells replicate (with the proliferating rate constant, *k*_{prol}) and then mature and differentiate via a chain of three transit compartments with a transition rate constant (*k*_{tr}) obtained from the mean maturation time (MMT) [*k*_{tr} = MMT/(number of transitions)] until they are circulating neutrophils that can be measured in the blood. Finally, circulating neutrophils can be eliminated by a first-order process. To gain homeostasis of the system, this elimination rate constant was assumed to be equal to *k*_{prol} and *k*_{tr}. At the same time, k_{prol} was influenced by 1) E_{drug} in an inhibitory manner and 2) a feedback mechanism which was responsible for the neutrophil recovery after drug administration. An exponential model was assumed for interindividual variability in SL and MMT and for residual variability.

For the external evaluation of this PK/PD model, the same steps used for the PK model were performed. As in the original model, the first-order estimation method was used.

#### PK/PD Model Optimization.

Based on the external evaluation of the PK/PD model and the observed long-term decay in neutrophil concentrations over the cycles, a pattern of cumulative neutropenia over the repeated treatment cycles was detected. To account for this pattern, PK/PD model extensions of the Friberg PK/PD model were investigated. Two published models (Bender et al., 2012; Mangas-Sanjuan et al., 2015) were identified (models A and B, respectively) and a new PK/PD model was developed and evaluated (model C). Figure 1 provides a schematic representation of all three models, including their modifications from the Friberg PK/PD model.

Briefly, model A (Fig. 1A) (based on Bender et al., 2012) was originally developed to describe thrombocytopenia in patients treated with trastuzumab emtansine. In this model, the circulating cells (originally thrombocytes, here neutrophils) at baseline [BASE_{tot}(*t*_{0})] were assumed to belong to two different subpopulations (BASE_{1} and BASE_{2}). BASE_{2} was affected by a second time-dependent drug effect (E_{drug2}), leading to a reduction in the total baseline over time [BASE_{tot}(*t*)], which was the target value of the feedback mechanism. While originally the average concentration of trastuzumab emtansine was determining the extent of this second drug effect E_{drug2}, here T_{C.0.05} of paclitaxel of each cycle was assumed to drive this second drug effect. In addition to the PD parameters estimated in the original model (MMT, *γ*, SL), the fraction of BASE_{tot}(*t*_{0}) not affected by the second drug effect (frB) as well as the depletion rate constant of *E*_{drug2} on BASE_{2} (*k*_{depl}) were estimated.

Model B (Fig. 1B) (based on Mangas-Sanjuan et al., 2015) was originally developed to describe neutropenia after different doses of diflomotecan. Compared with the Friberg model, two additional compartments accounting for a cyclic cell pathway from Prol to Prol with two stages of nonproliferating, nonmaturating quiescent cells were implemented (Q1 and Q2). Thus, two additional parameters were estimated: *F*_{prol}, which is the fraction of proliferating cells entering the maturation chain and not transitioning to the quiescent stages; and *k*_{cycle}, which is the circulation rate constant within the quiescent cell cycle.

For physiologic plausibility, model C (Fig. 1C) was additionally newly developed and also mimicked slowly replicating pluripotent stem cells in the bone marrow as a prior additional compartment (Stem) in the chain of granulopoiesis. Their proliferation was controlled by a different proliferation rate constant (*k*_{stem}). Model C was described by the following set of ordinary differential equations:(4)(5)(6)(7)(8)(9)in which Stem, Prol, Transit1, Transit2, Transit3, and Circ represent different stages of granulopoiesis. Both proliferation rate constants (*k*_{prol} and *k*_{stem}) were affected by the same PD parameters (SL) and feedback (FB) mechanism. To ensure homeostasis of the system without therapy, *k*_{prol} was estimated as a fraction f_{tr} of the transition rate constant *k*_{tr} as follows:

Hence, f_{tr} is the fraction of *k*_{prol} (f_{tr}) and *k*_{stem} [i.e., the remaining fraction (1 − f_{tr})] of the sum of both input processes for Prol and (due to equilibration) the corresponding fractions of *k*_{tr}. Thus, f_{tr} determines the ratio between *k*_{stem} and *k*_{prol} [*k*_{stem}/*k*_{prol} = (1 − f_{tr})/f_{tr}]. As approximations, if f_{tr} is estimated to be >0.5, then *k*_{stem} becomes smaller than *k*_{prol}. By contrast, if f_{tr} = 1, then *k*_{stem} = 0, which would result in the original Friberg model. This parametrization enabled that *k*_{tr} could still be computed as in the original Friberg model, by estimating MMT.

To describe baseline neutrophil concentrations, individual pretreatment concentrations were used, allowing for residual variability (baseline method 2; Dansirikul et al., 2008). This individual baseline value (Circ_{0}) was also the initial condition for all PD compartments of the differential equations in all PK/PD models [i.e., Stem(*t*_{0}) = Prol(*t*_{0}) = Transit1(*t*_{0}) = Transit2(*t*_{0}) = Transit3(*t*_{0}) = Circ(*t*_{0}) = Circ_{0}], except in model B for Prol [Prol(*t*_{0})= Circ_{0}/*F*_{prol}] and the quiescent cell compartments [Q1(*t*_{0}) = Q2(*t*_{0}) = (1 − F_{prol}) × Prol(*t*_{0})] as described in Mangas-Sanjuan et al. (2015). If two measurements were available before the first drug administration (baseline visit and the first day of the first cycle), the mean of the two values was used as the individual neutrophil baseline. In addition, for model C, other baseline methods suggested by Dansirikul et al. (2008) were investigated. Furthermore, an *E*_{max} model was evaluated, instead of a linear relation between drug concentration and drug effect. Finally, the number of transit compartments was varied and investigated.

An exponential model was applied for interindividual variability as well as for residual variability. A first-order conditional estimation method with interaction was used for all three models.

#### Comparing Models A–C: Model Evaluation and Model Performance.

Model evaluation and selection was based on parameter precision [relative standard error (RSE)], condition number (ratio of the highest to lowest eigenvalue of the correlation matrix), goodness-of-fit plots of individual predictions, population predictions, and weighted residuals, as well as on pcVPCs. Furthermore, if models were not nested, the Akaike information criterion (AIC) was used: AIC = OFV + 2 × k, in which OFV is the objective function value and k is the number of model parameter estimates. In addition, physiologic plausibility was considered to ensure reliable extrapolations for subsequent simulations of different dosing scenarios.

Deterministic simulations were performed to explore the performance of models A–C for a typical patient [56-year-old man with a total bilirubin concentration of 7 *µ*M; BSA of 1.8 m^{2}, and baseline neutrophil concentration of 6.48 × 10^{9} cells/l (median baseline of males, experimental arm)] undergoing 3-weekly dosing of 185 mg/m^{2} for six cycles using the typical parameter estimates of each model. Cell concentrations were explored not only in the circulating cell compartment but for all compartments. For model A, the second drug effect was either assumed to be active only until the end of the last cycle (3 weeks after the last dose administration) or to continue beyond the end of the observation period.

The relative change in the highest and lowest neutrophil concentrations (nadir) from the first to sixth cycles of each model was calculated based on this deterministic simulation and used to quantify cumulative neutropenia as follows:

(12)(13)#### Differentiation between the Effect of Paclitaxel and the Platinum-Based Drugs.

Model C was eventually further optimized by considering the combination of paclitaxel with carboplatin or cisplatin. None of the patients received carboplatin and cisplatin within one cycle. However, some of the patients receiving cisplatin were switched to carboplatin during the treatment. To account for both platinum-based drugs, two published structural and covariate PK models (de Jongh et al., 2004; Lindauer et al., 2010) were integrated. Since no drug concentration measurements were available, variability was neglected. Thus, the typical concentration-time profiles of carboplatin and cisplatin, respectively, were retrieved and these population predictions were used to estimate two additional slope factors for carboplatin and cisplatin. An additive drug effect was assumed for the combination therapy as follows:(14)in which SL_{platin} and *C*_{platin}(*t*) are the slope factor and the plasma concentration at time *t* of the coadministered platinum-based drugs, while SL and CPTX(t) are the slope factor and plasma concentration at time t of paclitaxel, respectively.

This combination PK/PD model was evaluated as described for models A–C. In addition, as described for the original PK/PD model, an external model evaluation was performed using the data from the conventional study arm. For this step, a pcVPC was generated, simulating 1000 data sets with the combination PK/PD model (in a simultaneous PK/PD analysis) considering interindividual variability in the PK and PD parameters.

#### Software for Data Analysis.

All modeling and simulation activities were performed in NONMEM (version 7.3) in combination with PsN software (version 4.4.0; Uppsala Pharmacometrics, Uppsala, Sweden). Data set preparation and analysis of the results was performed in *R* (version 3.2.2; R Project for Statistical Computing, Vienna, Austria) including the vpc R package (version 0.1.1). Deterministic simulations were done with Berkeley Madonna software (version 8.3.18; University of California, Berkeley, CA).

## Results

#### Data Analysis: CEPAC-TDM Study.

The two different dosing strategies of the two study arms led to a higher median paclitaxel dose of 315 mg (range, 205–438 mg) in the conventional arm versus 270 mg (range, 111–505 mg) in the experimental arm. In addition to paclitaxel, patients received a platinum-based drug, primarily carboplatin (149 versus 33 and 153 versus 30 in the conventional and experimental arms, respectively). Of 63 patients in the cisplatin group, 9 were switched to carboplatin cotreatment during the therapy, due to cisplatin-related toxicity.

A total of 658 paclitaxel PK samples were obtained from the experimental arm (no PK measurements available from the conventional arm), whereas 1635 and 1639 PD measurements were available for the conventional and experimental arms, respectively.

#### External PK Model Evaluation and Optimization Using Prior Information.

Paclitaxel concentrations from the CEPAC-TDM clinical trial were used to evaluate the original PK model for its applicability to our data. This external PK model evaluation showed a good prediction of the individual measurements (Supplemental Fig. 2). However, high shrinkage was observed (interindividual variability: VM_{EL}, 47.4%; Km_{TR} 92.9%; VM_{TR}, 92.3%; *k*_{21}, 97.6%; V_{3}, 58.5%; and Q, 52.8%; interoccasion variability: V_{1}, 98.6%; and VM_{EL}, 57.2%; and residual variability, 97.4%). Moreover, underprediction of the paclitaxel plasma concentrations was observed at the population level (pcVPC; Fig. 2A), especially in the relevant target concentration range of 0.05 *µ*M .

To account for the observed misspecification and to improve the predictivity of the PK model, the MAP Bayesian approach (Gisleskog et al., 2002) was applied combining the original PK parameter estimates (Table 1) as prior information with the newly obtained CEPAC-TDM data. For most of the parameters, only slight differences in the newly estimated PK parameter were observed compared with the original PK parameters. However, residual variability and two fixed-effects parameters, the concentration at half of the maximum elimination capacity and the covariate effect of bilirubin, were outside the 95% confidence interval of the original PK model, although with overlapping confidence intervals. Covariate parameters in general were estimated to be less influential on VM_{EL}. The highest relative changes in parameter estimates were observed for the volume of distribution of the second peripheral compartment (V_{3}) and for the intercompartmental clearance between the central and second peripheral compartments (Q) (+9.5% and +7.7%, respectively), although they were still within the 95% confidence interval of the original parameters. The higher estimates of V_{3} and Q translated into a larger paclitaxel distribution to the second peripheral compartment and a reduced distribution back to the central compartment. These estimates, together with lower elimination (due to an increased concentration at half of the maximum elimination capacity), led to higher predicted paclitaxel concentrations, especially in the later phases of the concentration-time profile resulting in higher exposure. These changes in PK parameters resulted in an improved description of the paclitaxel concentration-time profiles as illustrated in Fig. 2B.

Aside from the population level, differences in the individual prediction resulting from the two PK parameter sets (original and optimized) were assessed by comparing their EBEs and individual exposure parameters (Fig. 2C). For EBEs, V_{3} showed the highest difference in the EBEs (MRPE_{p}, −14.7%), whereas overall higher individual exposure estimates of T_{C>0.05} and area under the curve (MRPE_{p}, −5.82% and −4.28%, respectively) were revealed for the optimized PK model.

#### External PK/PD Model Evaluation.

Subsequently, the EBEs obtained with the optimized PK model were used to predict the individual PK concentration-time profiles, which eventually informed the PD response using the original PD parameters of the Friberg model. This PK/PD model overpredicted the neutrophil concentrations (and thus underpredicted neutropenia severity) during the whole observed period (pcVPC; Fig. 3A); overprediction was already present for day 15 values in the first cycle and increased for day 1 and day 15 measurements over the cycles. Although cumulative neutropenia led to decreasing neutrophil concentrations on days 1 and 15 over the cycles, the model predicted increasing neutrophil concentrations on day 15 instead, since the dose was reduced during therapy. This led to an increased discrepancy between the observed neutrophil data and model prediction over time.

#### PK/PD Model Optimization.

To better describe cumulative neutropenia over repeated treatment cycles, the PD part of the model was optimized. We evaluated the performance and behavior of three different structural PK/PD models, two from the literature (models A and B; Fig. 1, A and B) and one newly developed for physiologic plausibility (model C; Fig. 1C).

Parameters of all models were estimated with high precision (RSE <15%; Table 2), except for two parameters for model B (SL and *k*_{cycle}; RSE >25%). Overparameterization was not detected for any of the models (condition number <130). The AIC favored model C (lowest value with >100-point differences in AIC to models A and B), followed by model A and then model B.

Slope factors between the three models were highly variable, with the highest value for model B. In addition, the interindividual variability in SL was the highest for this model. For model A, the estimated frB implied that BASE_{tot}(*t*) decreased exponentially to approximately 50% of its original BASE_{tot}(*t*_{0}) (half-life of approximately 15 days, based on *k*_{depl} and assuming T_{C>0.05} to be 30 hours), which was the trigger for homeostasis and affected all cell compartments. For model B, the estimates of F_{prol} and *k*_{cycle} indicated that approximately 30% of the cells in Prol continued the maturation process, whereas 70% underwent additional circulation as quiescent cells. In this circulation, the mean transition time was 43.3 hours (=3/k_{cycle}). The estimate of f_{tr} in model C (determining *k*_{stem} and *k*_{prol} as a fraction of *k*_{tr}) expressed that the replication in Stem was approximately 73% slower than in Prol. Interindividual variability in MMT was negligible in the optimized PK/PD models, due to small estimates (coefficient of variation <10%) and high *η* shrinkage (>50%).

The pcVPCs (Fig. 3, B–D) showed that all three PK/PD models captured some of the observed pattern of cumulative neutropenia, but to a different extent. Models A and C described the observed neutrophil concentrations similarly well compared with model B, especially in the later cycles.

Furthermore, model C was not improved by implementation of an *E*_{max} model compared with a linear drug effect model. In addition, different baseline methods and variations in the number of transit compartments were explored but did not improve the model in terms of objective function value, parameter precision goodness-of-fit plots, and pcVPCs.

Regarding the deterministic simulation with constant doses of 185 mg/m^{2} for six cycles (Fig. 4), model A predicted approximately the same maximum value for the first (baseline) and second cycles; thereafter, the relative change in maximum value over all six cycles was the highest (31.4%) comparing models A–C. The relative change in nadir value (36.5%) was between models B and C. Model B predicted the lowest degree of cumulative neutropenia (relative change between first and last cycle in nadir and maximum value of 7.65% and 17.0%, respectively), despite the most significant initial relative change in maximum value (from the first to second cycle); from the second cycle on the maximum value was not changing further. Model C predicted a relative change of the maximum value of 27.6%, which was between models A and B; at the same time, the highest cumulative neutropenia effect for the nadir value (relative change of the nadir value of 67.8%) was observed. This high degree of cumulative neutropenia resulted in grade 2 neutropenia in the first cycle, grade 3 in the second cycle, and grade 4 from the fifth cycle onward. For concentrations in the stem cell compartment of model C, recovery was not achieved within a cycle length of 3 weeks (Supplemental Fig. 3).

Predictions of the three PK/PD models using deterministic simulation for neutrophil values after the end of the sixth treatment cycle varied substantially (Fig. 4). For model A, two possibilities regarding how to handle the second drug effect were investigated: considering E_{drug2} either 1) only during the six cycles (18 weeks) or 2) continued after treatment. If continued, BASE_{tot}(*t*) did not recover and the pretreatment baseline was not reached again. If not continued after 18 weeks, the original baseline was reached approximately 5 weeks after the last dose. Model B also showed a fast recovery of the neutrophil concentration to the pretreatment baseline with almost no oscillation (within 3 weeks after end of cycle 6), even though the *γ* exponent of the feedback function was the highest. This might be caused by the reduced proportion of cells entering the maturation chain. All models (except if continuation of E_{drug2} in model A is assumed) predicted that baseline would be reached after the end of therapy, but the time required was longer for model C (approximately half a year after the last drug administration), whereas for model A (assumption of no E_{drug2} after the last cycle) and model B baseline was reached again after approximately 5 weeks after the last dose. Overall, the deterministic simulation of model C showed the most physiologic prediction, given that baseline recovery in model A was dependent on the assumption taken for E_{drug2,} and model B predicted a decrease in the neutrophil concentrations only for the first cycle.

In summary, although all models described the observed cumulative neutropenia, the deterministic simulation of model C was the most plausible and the model performance was also superior regarding the AIC.

#### Differentiation between the Effect of Paclitaxel and the Platinum-Based Drugs.

Model C was further expanded to distinguish between the drug effect of paclitaxel and platinum-based drugs. Based on the pcVPCs (Fig. 3E), the performance of the expanded model was similar to model C; nevertheless, the AIC decreased by 42.5 points. Parameters were estimated with sufficient precision (RSE <20%) for all parameters, but those of the slope factors slightly increased (RSEs of 16.9% and 26.1%; Table 2). The condition number (113) did not indicate an overparameterization of the model. Estimated system PD parameters remained within the range of those of model C. Only the drug effect of paclitaxel was reduced by 60.5% compared with model C, given that the effect was now split into two effects, for paclitaxel and the platinum-based drugs. Interindividual variability on the slope factor of paclitaxel was increased compared with model C. Although the slope factor of carboplatin was estimated in a reasonable range, the slope factor of cisplatin was implausibly high (32.1 l/mg), very likely due to the low number of patients (*n* = 33) treated with this drug. Hence, dose individualization based on this model should be recommended only for the combination of paclitaxel with carboplatin.

Finally, the optimized PKPD platinum combination model, including separate drug effects for paclitaxel and the platinum-based drugs, was used to predict the neutrophil concentration-time profiles of the conventional arm of the CEPAC-TDM study (pcVPC; Fig. 3F). A good prediction of the median and 95^{th} percentile was achieved, although an underprediction of the low concentration-time profiles (5^{th} percentile) was visible but still the same neutropenia grade was predicted.

## Discussion

Due to the narrow therapeutic window and the high interindividual PK and PD variability of paclitaxel, model-based dose individualization can help to reduce toxicity without compromising efficacy, as shown in a first step in the CEPAC-TDM study (Joerger et al., 2016) for neuropathy but not for neutropenia. To further improve therapy for patients suffering from severe neutropenia, we evaluated, optimized, and expanded the PK/PD model on which the dosing algorithm for the CEPAC-TDM study was developed. Thus, the Friberg neutropenia PK/PD model was extended to describe cumulative neutropenia after repeated chemotherapy cycles in combination with cisplatin/carboplatin.

The individual PK predictions of the original model matched the paclitaxel concentrations measured in the CEPAC-TDM study well, due to the sparse sampling and the high interindividual variability in the model. However, we observed an underprediction when we evaluated the population level, pointing toward suboptimal original PK parameter values for the investigated population. Since PK/PD modeling was undertaken sequentially, PK misspecifications may also affect PD. Therefore, PK model optimization was performed using the MAP Bayesian approach (Gisleskog et al., 2002). The re-estimated parameters described the paclitaxel concentrations in the terminal phase of the concentration-time profile better by a faster distribution to the second peripheral compartment and slower elimination. Furthermore, the optimized PK parameter set predicted higher exposure for the population and individuals. Overall, the re-estimated PK parameters were in line with previously published parameters, with a high volume of distribution (>100 liters) (Gianni et al., 1995; van Zuylen et al., 2001), fast elimination (clearance >30 L/h) (Wiernik et al., 1987; van Zuylen et al., 2001), and improved predictive performance. Thus, the re-estimated PK parameters were used for the subsequent steps. An external evaluation of the Friberg neutropenia model was performed with two major findings. First, overprediction of neutrophil concentrations at day 15 (i.e., less toxic effect on neutrophils), already in the first cycle, pointed toward the need for a higher slope factor. Second, cumulative neutropenia is often seen in clinical practice after repeated cycles with cytotoxic drugs, but only rarely investigated and also not accounted for in the Friberg model. We hypothesized that this pattern was caused by bone marrow exhaustion (BME) due to the damage of early primitive bone marrow stem cells, such as pluripotent long-term hematopoietic stem cells (Mauch et al., 1995).

Since the Friberg model for neutropenia predicted approximately the same nadir concentration and thus extent of neutropenia for each cycle if doses were kept constant, the original model was not capable of describing the observed cumulative neutropenia and signs of BME during repeated chemotherapy cycles. Three different structural models were investigated (models A–C) to improve the predictivity of the PK/PD model for long-term treatment, as in the CEPAC-TDM study. Model A (adapted from Bender et al., 2012) followed a more empirical approach: by decreasing a part of BASE_{tot}(*t*) time dependently, the feedback effect and thus the recovery of the neutrophil concentration was reduced. Initially, model B (adapted from Mangas-Sanjuan et al., 2015) was not intended for describing cumulative myelosuppression, instead thrombocytopenia after different dosing schedules was aimed to be described. Nevertheless, this model predicted lower neutrophil concentrations than the Friberg model for repeated dosing. Model B used a more mechanistic approach by implementing two stages of quiescent cells accounting for progenitor cells that can potentially replicate but rest and therefore are not harmed by the treatment (comparable to the G_{0} phase in the cell cycle). Finally, a new model (model C) was developed to implement the BME hypothesis by adding a stem cell compartment as the first compartment in the maturation chain. In accordance with physiology, proliferation of the stem cells was slower (73%) but was affected by the same drug and feedback effect as cells in Prol.

All three PK/PD models were successfully implemented and parameters were estimated with high precision. The slope factor of these models was estimated to be higher than in the original model, whereas high differences between them were observed but could not be compared across models due to differences in the implementation of the drug effect.

Cumulative neutropenia was described by all three models to a certain extent; however, the potential for extrapolation of models A and B may be limited. For model A, the parameter estimate *k*_{depl} for the time-dependent effect is specific for the treatment cycle length of the underlying data, and predictions after the end of treatment strongly depend on the assumptions for the duration of E_{drug2}. For model B, full recovery and equilibrium was reached very quickly after the last dose. Thus, an increased cycle length would reduce or let disappear the cumulative effect, which is critical, since oncologists often empirically reduce dosing (as in the conventional study arm). Model C, on the other hand, described cumulative neutropenia in a more semimechanistic way, supporting the BME hypothesis. For nadir and maximum values of the neutrophil concentrations, a decrease over the cycles was observed, but maximum values were more affected, as supported by the data. The inability of the stem cell compartment to recover within a cycle mimicked the aforementioned BME hypothesis that long-term hematopoietic stem cells are damaged. Note that although stem cell proliferation was affected by the same drug effect parameter [E_{drug} = SL × *C*_{PTX}(*t*)], the overall toxic effect for each of the two cell types in bone marrow is dependent on the proliferation rate [i.e., *k*_{prol} × (1 − E_{drug}) and *k*_{stem} × (1 − E_{drug})] (Berenbaum, 1972), which is smaller for stem cells compared with progenitor cells (Steinman, 2002). Hence, stem cells were ultimately assumed to be less susceptible, in accordance with the literature (Harrison and Lerner, 1991). Damage to these cells was especially seen for the second and following cytotoxic dose administrations, when the cells are more vulnerable by stimulated proliferation due to the hematopoietic stress caused by the first drug administration (Harrison and Lerner, 1991; Trumpp et al., 2010). The BME hypothesis is also supported by the fact that the platelet concentrations showed the same decreasing pattern over time (Supplemental Fig. 4), which is plausible when the stem cells, from which neutrophils and platelets originate, are disturbed. Model C further predicted a long span (approximately half a year after the last dose) for the stem cells to fully recover. Since no neutrophil measurements were available for such a long period, the hypothesized long recovery time for the entire bone marrow beyond the sixth cycle represent extrapolations needing further validation. Unfortunately, this period is typically not monitored even though baseline is not achieved, since patients rapidly recover to non-neutropenic grade 0. Nevertheless, a long recovery time seems plausible, since previous chemotherapies have been identified to reduce baseline (Kloft et al., 2006), supporting the assumption that the bone marrow could still be disturbed. Model C also showed good agreement of the system-related parameters (MMT and *γ*) between the original model and the literature (Friberg et al., 2002), indicating that previous knowledge gained from the Friberg model, could be used to inform model C. Overall, the best description of the neutrophil data was obtained using models A and C. However, model C integrated the most physiologic explanation for the effect and has the potential to describe long-term treatment with paclitaxel.

Thus far, only the neutropenic effect of paclitaxel had been considered, but combination therapy with cisplatin/carboplatin, known to have an impact on neutropenia (Go and Adjei, 1999), was administered. Differentiation of drug effects is necessary for future simulation, in which dose recommendation of paclitaxel and the platinum-based components shall be adapted. Thus, model C was extended and the platinum combination model assuming additive drug effects slightly improved the model prediction, while allowing for a more realistic drug effect.

The estimated paclitaxel slope in the combination PK/PD model was lower than in model C, which is reasonable since the drug effect was divided into two components (paclitaxel and platinum-based drugs). Compared with the literature, paclitaxel and carboplatin slope factors were approximately 2-fold higher, as previously reported for the respective monotherapy, with values of 2.21 l/*µ*mol (Friberg et al., 2002 and 1.85 l/*µ*mol (Joerger et al., 2007a) for paclitaxel and 0.460 l/mg (Schmitt et al., 2010) for carboplatin. Deviations point toward a synergistic pharmacodynamic interaction, which is not considered in the additive assumption underlying this model. Synergism was suggested previously (Choy et al., 1998; Engblom et al., 1999) and is mechanistically plausible: paclitaxel increases the proportion of cells in the G_{1} phase of the cell cycle, where cells become more sensitive to carboplatin (Long and Fairchild, 1994; Engblom et al., 1999). However, more complex models (e.g., general PD interaction model; Wicha et al., 2016) were not supported by the data. Hence, the simplified (additive effect model) was selected, but predictions of neutropenia grade from extrapolated combination drug concentrations should be regarded with care, since the synergistic interaction depends on both drug concentrations. Despite this limitation, the platinum combination model well described the data from the conventional study arm as well, indicating good predictive performance.

In conclusion, a comprehensive PK/PD model able to describe and predict cumulative neutropenia after paclitaxel combination therapy was developed. Describing this long-term effect semimechanistically improved the understanding of neutropenia and provided further evidence of the BME hypothesis that the effect on pluripotent stem cells might cause cumulative neutropenia. Because of the mechanistic character of the model, the framework can be applied to other myelotoxicity drugs. In addition, the model can differentiate between the effects of paclitaxel and carboplatin, which allows for better predictions of dose adaptations for patients with non–small cell lung or ovarian cancer. Further data assessment with other doses of the combination therapy as well as neutrophil concentration evolvement after the end of treatment should be performed. This model showed a high predictive performance for the conventional arm in which no drug concentration had been determined. In this work, we used a confirmatory phase IV study not only to assess the dosing algorithm but also for further learning to improve understanding, as suggested by Sheiner. Our paclitaxel-carboplatin combination PK/PD model provides the basis for further explorations of different dosing strategies to increase our knowledge and thus to contribute to further individualized anticancer treatment.

## Acknowledgments

The authors thank the High-Performance Computing Service of ZEDAT (Zentraleinrichtung für Datenverarbeitung) at Freie Universitaet Berlin (https://www.zedat.fu-berlin.de/HPC/EN/Home) for computing time, the Central European Society of Anticancer Drug Research for providing the study data, and Dr. Niklas Hartung for discussions about prior information.

## Authorship Contributions

*Participated in research design:* Henrich, Joerger, Huisinga, Kloft, Parra-Guillen.

*Performed data analysis:* Henrich.

*Wrote or contributed to the writing of the manuscript:* Henrich, Joerger, Kraff, Jaehde, Huisinga, Kloft, Parra-Guillen.

## Footnotes

- Received February 28, 2017.
- Accepted June 5, 2017.
C.K. and W.H. are supported by grants from an industry consortium comprising AbbVie Deutschland GmbH & Co. KG, Boehringer Ingelheim Pharma GmbH & Co. KG, Grünenthal GmbH, F. Hoffmann-La Roche Ltd., Merck KGaA, and Sanofi. C.K. is supported by grants from the Innovative Medicines Initiative Joint Undertaking (DDMoRe: Drug Disease Modelling Resources) and Diurnal Ltd.

C.K. and Z.P.P.-G. are co-senior authors.

Part of this work was previously presented at the following meetings:

Henrich A, Joerger M, Huisinga W, Kloft C, and Parra-Guillen ZP (2015) External evaluation of a PK/PD model describing the time course of paclitaxel and neutropenia in patients with advanced non-small cell lung cancer.

*24th Population Approach Group Europe (PAGE);*2015 Jun 02-05; Hersonissos, Greece.Henrich A, Joerger M, Huisinga W, Kloft C, and Parra-Guillen ZP (2015) Dose-individualisation of paclitaxel in patients with advanced non-small cell lung cancer: exploiting modelling and simulation to refine dosing algorithm.

*Annual Meeting of the Central European Society of Anticancer Drug Research (CESAR)*; 2015 Sep 17-19; Innsbruck, Austria.↵Henrich A, Joerger M, Kraff S, Jaehde U, Huisinga H, Kloft C, and Parra-Guillen ZP (2016) PK/PD model extension to characterise bone marrow exhaustion in cancer patient making use of a prior paclitaxel PK model.

*25*2016 Jun 07-10; Lisbon, Portugal. PAGE 25: 5879 [www.page-meeting.org/?abstract=5879], (2016).^{th}Population Approach Group Europe (PAGE);Henrich A, Joerger M, Huisinga W, Parra-Guillen ZP, and Kloft C (2016) Semi-mechanistic PK/PD modelling to characterise long-term deterioration of neutropenia in cancer patients.

*Annual Meeting of the Central European Society of Anticancer Drug Research (CESAR)*; 2016 Sep 08-10; Munich, Germany.Henrich A, Joerger M, Huisinga W, Kloft C, and Parra-Guillen ZP (2016) Characterisation of bone marrow exhaustion and long-term neutropenia by applying PK/PD modelling.

*Annual Meeting of the Deutsche Pharmazeutische Gesellschaft (DPhG);*2016 Oct 04-07; Munich, Germany.This article has supplemental material available at jpet.aspetjournals.org.

## Abbreviations

- AIC
- Akaike information criterion
- BASE
- baseline
- BME
- bone marrow exhaustion
- BSA
- body surface area
- CEPAC-TDM
- Central European Society for Anticancer Drug Research Study of Paclitaxel Therapeutic Drug Monitoring
- EBE
- empirical Bayesian estimate
- E
_{drug} - drug effect
- MAP
- maximum a posteriori
- MMT
- mean maturation time
- MRPE
_{p} - median relative prediction error for parameter p indicating bias
- pcVPC
- prediction-corrected visual predictive check
- PD
- pharmacodynamic/pharmacodynamics
- PK
- pharmacokinetic/pharmacokinetics
- RSE
- relative standard error
- SL
- slope factor
- T
_{C>0.05} - time of paclitaxel plasma concentration above the threshold of 0.05
*µ*M - VM
_{EL} - maximum elimination capacity

- Copyright © 2017 by The American Society for Pharmacology and Experimental Therapeutics