Model-Based Individualized Treatment of Chemotherapeutics: Bayesian Population Modeling and Dose Optimization

6-Mercaptopurine (6-MP) is one of the key drugs in the treatment of many pediatric cancers, auto immune diseases and inflammatory bowel disease. 6-MP is a prodrug, converted to an active metabolite 6-thioguanine nucleotide (6-TGN) through enzymatic reaction involving thiopurine methyltransferase (TPMT). Pharmacogenomic variation observed in the TPMT enzyme produces a significant variation in drug response among the patient population. Despite 6-MP’s widespread use and observed variation in treatment response, efforts at quantitative optimization of dose regimens for individual patients are limited. In addition, research efforts devoted on pharmacogenomics to predict clinical responses are proving far from ideal. In this work, we present a Bayesian population modeling approach to develop a pharmacological model for 6-MP metabolism in humans. In the face of scarcity of data in clinical settings, a global sensitivity analysis based model reduction approach is used to minimize the parameter space. For accurate estimation of sensitive parameters, robust optimal experimental design based on D-optimality criteria was exploited. With the patient-specific model, a model predictive control algorithm is used to optimize the dose scheduling with the objective of maintaining the 6-TGN concentration within its therapeutic window. More importantly, for the first time, we show how the incorporation of information from different levels of biological chain-of response (i.e. gene expression-enzyme phenotype-drug phenotype) plays a critical role in determining the uncertainty in predicting therapeutic target. The model and the control approach can be utilized in the clinical setting to individualize 6-MP dosing based on the patient’s ability to metabolize the drug instead of the traditional standard-dose-for-all approach.


Introduction
Cancer has been a perennial challenge to clinicians and researchers ever since its existence has come to light. Despite the remarkable progress in healthcare in the past few decades, cancer remains the second leading cause of deaths, accounting for about 1 in 4 deaths in the US [1]. One of the primary reasons for the failure of cancer treatment can be attributed to high interpatient variability in response to such treatment. The existing treatment modalities are "effective" only in subsets of patient population due to a significant genetic and phenotypic variation among patients. It is this intrinsic variation, even within the same genotypic group, which renders the clinical decision far from straightforward. The dose regimen determined during randomized clinical trials involving a small number of patients may be appropriate at best for an "average" patient because these studies are designed to define the best dose for the whole population, and not for any specific patient. Thus, it produces severe toxicity in some patients and insufficient treatment in others. Adverse drug reactions has been estimated to cause about 2 million hospitalizations and 100,000 deaths per year in the US [2], necessitating a dire need for a rational approach to individualized treatment.

Challenges of 6-MP Treatment
6-MP is one of the important drugs in a series of purine analogues. In addition to many pediatric cancers, 6-MP is a key drug in inflammatory bowel disease (IBD) and many autoimmune diseases. Common acute side-effects during 6-MP treatment include myelosuppression, pancreatitis, gastrointestinal intolerance and hepatotoxicity. Besides acute side-effects, many clinical studies reported several chronic effects related to 6-MP treatment. For instance, in patients who have undergone 6-MP treatment for acute lymphoblastic leukemia (ALL), recurrent ALL, secondary neoplasm and other multiple chronic medical conditions are prevalent [3,4]. Clinical studies show that inadequate therapy leads to recurrent ALL whereas aggressive treatment results in acute side-effects and secondary malignancies, thus calling for optimization and individualization of 6-MP dosing [4][5][6][7][8][9].
6-MP undergoes extensive intracellular metabolism to yield 6-thioguanine nucleotide (6-TGN) (active metabolite) and other methylated metabolites of mercaptopurine (MeMP) [10]. 6-TGN and MeMP are catalyzed by enzymes hypoxanthine-guanine phosphoribosyltransferase (HGPRT) and thiopurine methyltransferase (TPMT) respectively. The relative activities of HGPRT and TPMT are genetically transcribed and regulated for a given patient. Among these, TPMT enzyme activity appears to be the rate limiting step and hence dictates the net concentration of 6-TGN [11]. Much of the treatment variability during 6-MP treatment is cascaded down from the genetic polymorphism exhibited in the TPMT gene [12]. For instance, in patients with high TPMT activity, the 6-TGN pathway is suppressed, resulting in low 6-TGN concentration and hence treatment failure. On the other hand, in patients with low TPMT activity, the 6-TGN metabolism is elevated which eventually results in life-threatening myelotoxicity. Thus, TPMT genetic polymorphism is highly correlated with treatment outcome [5,[13][14][15]. Hence, the utilization of the TPMT genotype as a pharmacogenetic marker has been suggested to guide the treatment protocol for an individual patient [16].

Genotyping vs. Phenotyping
Dosing decision for a specific patient can be made with any one of the following biological markers associated with a drug-disease combination: i) DNA sequence, ii) gene expression profile, iii) protein/enzyme concentration/activity, iv) metabolite/drug concentration and v) cell/clinical response. We refer to this as biological chain-of-response. Pharmacogenomics, at the crossroads of genetics and pharmacology, sheds light on the causative genetic variants influencing treatment response under drug intervention. When introduced, pharmacogenomics was believed to alleviate many of the issues arising in current medical sciences, especially personalized treatment [17,18]. However, it has produced limited success in some of the drug-disease applications that were deemed to be classical cases for pharmacogenomics based personalization [19,20]. Although pharmacogenomics provides some information vital for predicting treatment outcome in the subgroups of patient population, its ability to explain within-group variation may be in question. A detailed account this subject can be found in [21].
Pharmacogenomics relies on a static snapshot of a specific DNA sequence or gene expression and assumes a deterministic evolution of biomolecular events resulting in a predictable cellular response for a given gene variant. In addition, it considers each gene as an independent causal factor for the observed response. However, human physiology is a complex dynamic system and cellular response is a manifestation of interplay between many levels of physiological processes and molecular entities. Furthermore, human physiology is complicated by homeostatic feedback loops, molecular crosstalk and bypass mechanisms that can lead to unexpected clinical outcomes. As a result, within a specific genotype, there is a distribution of phenotypes across patient population. For example, in the case of 6-MP, although there are only a few validated TPMT genotypes, there are as many enzyme activity levels as there are patients. To make things even more complex, within a specific enzyme activity range, there are as many active drug concentrations (the ultimate manifestation of the genotype) as there are patients. Hence, 6-TGN concentration, which is responsible for the cytotoxicity, should be regarded as the ultimate covariate for dose individualization.

Importance of Individualized Treatment
The foregoing perspective is more clearly evident from a case-study on a TPMT deficient patient who had undergone 6-MP treatment [22]. As soon as the treatment was initiated, the patient's vital cellular counts plummeted below the critical level. Following several treatment interruptions and, blood/platelet transfusions due to severe myelotoxicity, TPMT genotype and phenotype were assayed at week 20; 6-TGN concentration was also measured. Even with this information, it took another 45 weeks to arrive at the optimal dose through trial-and-error with additional treatment interruption and blood/platelet transfusions. It is clear that the current practice of patient titration and qualitative use of vital information is severely deficient. This trial-and-error process has a profound impact on immediate and long-term patient health and hence on allocation of limited clinical resources.
There is a growing body of literature that acknowledges this state of affairs and suggests the need for tailoring the dose regimen based on a patient's genetic and phenotypic make-up [8,[23][24][25]. Given the dynamic nature of physiological responses, it has to be an ongoing process rather than a 'study-and-adopt' approach. In other words, following a detailed analysis and accumulation of information during the study phase, a minimum of information must be obtained from each new patient to adapt the approach to the patient before making predictions and dose optimization. Given the significant limitation on continuous monitoring in the clinical settings, a robust model-based in silico approach, adaptable to individual patients, is indispensable. A recent report by the National Academy of Engineering and the Institute of Medicine highlights the potential of such engineering approaches, consummated through a partnership between healthcare professionals and engineers, in patient focused health care delivery [25]. Hence, these factors form the thrust of this manuscript and are summarized in Fig 1. In section 2, we describe the model and methodologies used and provide some important results in section 3. Finally in section 4, we conclude with discussion on the impact, constraints in clinical implementation and possible extension.

Problem Statement
A generic dynamic model for any given drug and patient can be described by a system of ordinary differential equations of the following form, where x i (t) 2 < n : state variables, θ i 2 < p : model parameters, u i (t) 2 < q : drug input, c 2 < n : initial conditions,ŷ i ðtÞ 2 < l : predicted model output. Patient population is characterized by the random parameter matrix Θ 2 < p × < N so that an individual patient, i, is represented by a vector θ i . Further, we assume that the model parameters θ i can be partitioned into highly sensitive parametersθ i 2 < s and less sensitive parametersθ i 2 < r ¼ < pÀs such that, Sensitive parameters are used to identify patient specific parameters whereas the rest are fixed at population means. Let y ij denote the j th measurement for the i th patient, at time t ij . It is not uncommon for these individual measurements y ij to be corrupted by measurement and assay error, besides model misspecification. If we assume that the modeling and experimental errors are additive in nature, then the observed concentration would be given by, For practical purposes, the errors can be assumed to be independent and normally distributed . ξ ij represents error model parameters. Unlike chemical systems and some biological systems, where the model parameters are more or less constant once the experimental conditions are fixed, human physiology is subject to high inter-subject variability. In other words, for a given disease and drug dosing, no two patients will respond in exactly the same way. As such, although the model structure can be assumed to be identical for individual patients, the model parameters vary significantly among patients. For the i th patient, the model parameters are given as, where y is the typical value of model parameters (population mean) and η i 2 < p are independent vectors representing the deviation of the i th patient's parameters from the population mean values. From Eq (5), it is clear that an individual patient is a part of the population described by a multivariate probability distribution, In the above equation, M p ð:; :Þ represents p-dimensional multivariate distribution and S 2 < p × < p is an inter-patient variance-covariance matrix. Much of the challenge in population modeling and dose individualization resides in robust estimation of Pðθ i Þ for a new patient within physiological and clinical constraints. A brief description of various methods available to determine Pðθ i Þ may be found in [24]. In this work, we are exploiting the Bayesian approach which is elaborated in [24,26,27]. Under the Bayesian framework, several steps are included: M p ð θ; SÞ is evaluated off-line with the availability of drug dose and drug concentration data from a large number of patients. Based on the variability in response propagated by M p ð θ; SÞ, optimal sampling times that minimize the uncertainty in parameter estimation are estimated. For each new patient, given a measurement at an optimal time, Pðθ i Þ and optimum dosing profiles are determined on-line.

Experimental Data
The data for characterizing 6-MP metabolism and cellular 6-TGN concentration were obtained from three separate clinical studies reported in the literature. In the first study (D1), the 6-TGN concentration was measured in a group of 23 patients undergoing 6-MP treatment [28]. Up to four data points were collected per patients over a period of eight weeks. The second set of data (D2) for 6-TGN concentration was collected from patients undergoing 6-MP treatment under "per protocol" group in [29,30]. 6-TGN concentration was measured from eight patients over 20 weeks. A third set of data (D3) was collected from a clinical study on ALL patients in which 102 measurements of TPMT enzyme activity and corresponding 6-TGN concentrations were collected [31]. This data set is primarily used to determine the dosing profile at the beginning of the treatment when 6-TGN measurement is not available and hence the dosing decision has to be made with TPMT enzyme activity measured at t = 0. D3 is also extensively used for evaluating uncertainty in phenotype prediction with various data along the chain-of-response. The data sets that showed a clear indication of treatment discontinuation or significant dose reduction were eliminated.

Modeling 6-MP Metabolism
Although the TPMT pharmacogenomics and the metabolism of 6-MP is one of the extensively studied systems in clinical pharmacology literature, utilization of quantitative dose optimization strategies are limited. Hawwa et al. [32] performed a population pharmacokinetics study of 6-MP in pediatric patients considering TPMT genotype as one of the main covariates to characterize inter-patient variability. Phenotypic variations and/or dose optimization strategies were not considered in their work. A simplified schematic of the 6-MP intracellular metabolism accounting for 6-TGN production is shown in Fig 2. Following oral intake to the gut, 6-MP is absorbed at the rate of k ab into the plasma where it undergoes extensive hepatic clearance at the rate of k el . From the plasma, 6-MP is transported into the intracellular space where it undergoes metabolic conversion. Since negligible intracellular concentration of 6-MP has been reported [33], we assume that 6-MP is metabolized as soon as it enters the intracellular space. The desired pathway leading to 6-TGN is initially catalyzed by HGPRT, followed by a series of other metabolic conversions at a lumped rate of k cm . 6-TGN is eliminated from the cells at a constant rate of k me .
The state variables are defined as follows: x g , amount of 6-MP in the gut (picomole (pmol)); x c , amount of 6-MP in the plasma (pmol); x m , concentration of metabolite 6-TGN in peripheral RBCs (pmol/8x10 8 RBCs). In the above model, x m is the observed variable. Conversion of 6-MP into 6-TGN follows Michaelis-Menten kinetics with reaction rate k cm and Michaelis-Menten constant K. Patient specific TPMT activity is represented as a quantity relative to its maximum level. Thus, where e denotes TPMT enzyme activity. It is observed in clinical studies that the production of 6-TGN is negatively correlated with TPMT activity [31]. Hence, it is assumed that a fraction of 6-MP, proportional to (1 − e rel ), is converted to 6-TGN. Consequently, the reaction rate is modeled as, The description and units of all the state variables and parameters are listed in Table 1.

Bayesian Parameter Estimation: Off-line
Estimation of Individual Patients' Parameter Distribution. The key step in the use of Bayesian population modeling for patient treatment individualization is the estimation of Bayesian posteriors for each patient using patient-specific data and a suitable prior. To accomplish this, first, parameter statistics are estimated to formulate a meaningful prior distribution for Bayesian calculations. The incorporation of prior knowledge is a key and unique aspect of the Bayesian framework. With large amount of data and well-defined parameters, the prior distribution may have a little impact on the posterior inference. However, when the data is far from adequate, such as in clinical settings, the prior distribution can play an important role, and hence it is expedient to devote considerable effort to obtain an informative prior. The model parameters for each of the individual patients who are part of the clinical study can be estimated using the maximum likelihood approach.
where LðD i jθ; xÞ is the likelihood function. With the assumption of experimental errors ε being independent and normally distributed, the likelihood function is given by, where D i is the data set obtained from the i th patient. Next, with the prior distribution formulated using this statistics, the posterior distribution can be calculated using Bayes theorem. According to Bayes theorem, the posterior distribution of the model parameters p i ðθ; ξjD i Þ is given by Eq (12).
where p(θ, ξ) is the prior distribution. The denominator, pðD i Þ, is the normalization factor equal to the expected value of the data irrespective of the parameters. The most common method for Bayesian inference is the Markov chain Monte Carlo (MCMC) sampling. However, MCMC can be very computationally intensive for large and complex models. An alternative methodology to MCMC sampling is the variational Bayes approximation. Variational Bayes translates the Bayesian inference into an optimization problem by approximating the posterior distribution using a known distribution form (e.g., a family of Gaussian distributions). The idea behind this methodology is that the logarithm of pðD i Þ can be separated into two elements as shown in Eq (13) to (15).
qðθ; ξjϕÞ ln qðθ; ξjϕÞ LðD i jθ; ξÞpðθ; ξÞ dθdξ ð14Þ Here, q(θ, ξ|ϕ) is the parametric probability distribution that approximates the posterior distribution pðy; xjD i Þ. ϕ is the set of parameters characterizing q (e.g., mean and covariance for a Gaussian distribution). If q(θ, ξ|ϕ) is free to be any probability function, then the maximum lower bound is obtained when the Kullback-Leibler (KL) divergence is zero. This occurs when q(θ, ξ|ϕ) exactly matches the posterior distribution [34]. Then, minimizing the KL divergence is equivalent to maximizing the lower bound of L. Therefore, the set of parameters Θ is determined as the one that maximizes L. Laínez et al. [35] developed a decomposition strategy to deal with the variational inference of models that are described by a set of DAE which will be followed for the variational Bayes approximation of the model for 6MP metabolism. Following this approach, the solution of the variational Bayes problem is decomposed into three steps: a maximum a posteriori optimization which is facilitated by using an orthogonal collocation approach, a preprocessing step which is based on the estimation of the eigenvectors of the posterior covariance matrix, and an expected propagation optimization problem [34]. The decomposition strategy has been implemented using the R [36] and GAMS [37] software packages. GAMS has been used for the optimization problems in the first and last steps.
Population Prior Distribution. One of the advantages of using the Bayesian approach is the systematic and prospective accumulation of information from each participating subject. Initially, patient-specific information from the retrospective clinical studies is accumulated to capture the population characteristics. This information is not only used for characterizing the new incoming patients but also constantly updated through the incorporation of new patients as part of the information. The individual posterior joint distributions obtained in the previous section are utilized to formulate an informative population prior distribution. A specific number of parameters (based on some weighting factor) from the converged posterior distribution are sampled to formulate the population prior distribution [24].
The weighting factor w i is usually chosen based on the quality and quantity of data used to obtain these posterior distributions. Since the original data contained varied number of measurements, w i is taken as a function of the number of data points for each patient.

Global Sensitivity Analysis
The model parameters for a new patient cannot be estimated accurately unless measurements are made on that patient. Although the availability of an informative prior helps to reduce the number of measurements required, its impact can be considerably reduced when there are a number of model parameters to be estimated. However, if we can systematically identify and reduce the number of parameters which must be estimated a priori, it will aid in the efficient estimation of parameters online, when measurements are limited. In addition, it will significantly reduce the computation time which is also an essential part of timely decision making for a new patient in the clinical setting. The basis for reducing the number of parameters to be estimated stems from the observation of uncertainty in the dynamical systems. Uncertainty in the model output is primarily propagated from the uncertainty in the model input (i.e. parameters). Although it is true that the uncertainties in model parameters will have some effect on the model output, not all parameters have the same level of influence [38]. Consequently, it can be expected that the uncertainty in the estimation of the highly sensitive parameters will have the most significant impact on the model prediction. Thus, it becomes critical to estimate the most sensitive parameters as accurately as possible with the limited data available. Although less sensitive parameters have little effect on the measured variable, they impart indirect effect through other auxiliary variables and hence are fixed at a nominal value, instead of being eliminated or neglected. The error involved in such an approximation can be estimated as follows [39].
where S tỹ v is the total sensitivity of the model output corresponding to parameterỹ v . For an arbitrarily small value of ε = 0.05, the probability of getting dðỹ v;0 Þ < 21S tỹ v is more than 0.95.
Since physiological model parameters vary over a wide range, we used global sensitivity analysis (GSA) for estimating S tỹ v . The Soboĺ method was used in this work to estimate the total sensitivity indices [39].

Optimal Experimental Design
Clinical data, especially for new patients, are constrained due to economical and logistical reasons. Hence, we are interested in characterizing patients by measuring drug concentration with a minimum number of samples. By collecting samples at optimal time points, dictated by the design of experiment principles (DoE), one can estimate the model parameters accurately. Obviously, these optimal sampling points should be determined for the whole population with a detailed study on a class of patients with full set of data. Here, the population characteristics are enriched in Mðθ; ξÞ and will be used for sampling time determination. The average drug profile among patients studied is given by Suppose we have a way to characterize the best choice of time at which these measurements are to be made, given the parameter vector θ. Denote this instant by T ðθÞ. The average of this time is given by Similarly, one could also consider T ðEθÞ as a potential choice for the optimal time. To evaluate these measures, more generally, we may propose θ Ã such that the variance below is minimized.
The solution of the foregoing minimization problem will yield θ Ã and hence the best time for the measurement by evaluating T ðθ Ã Þ. Eq (20) may be written further as Clearly, Eq (21) is minimum when T ðθ Ã Þ ¼ ET . With this conclusion at hand, T ðθÞ is found using D-optimal criteria which attempts to minimize the volume of the hyper-ellipsoid spanned by the joint parameter space. This is achieved by maximizing the determinant of the Fisher information matrix as in Eq (22) max In Eq 22, φ represents design parameters. The robust formulation is ensured by integrating the determinant over a representative population parameter space and weighted by the where S i is the sensitivity matrix defined as

Dose Optimization
With the availability of a patient-specific model, the final step in model-based individualized dosing is dose optimization to fulfill certain physiological objectives. The nature of the objective function selected is entirely dependent on the drug-disease-side effect combination. For 6-MP, the therapeutically effective range is defined in terms of 6-TGN concentration in peripheral RBCs. Moreover, the dosing strategy employed can depend on the degree to which the patient response can change over time. In some cases, the optimal dosing once determined can be repeated for a number of dosing cycles, unchanged [40]. In other situations, model mismatch and changes in patient response over time require repeated dose adjustments as patient response is tracked. This is the case in this work, thus, we utilize robust model predictive control (MPC) to optimize the dose due to its intrinsic capability to implement prediction and optimization under uncertainty [41][42][43]. In general, the MPC problem is formulated as solving the on-line finite-horizon, closed-loop optimal control problem subject to an underlying model and constraints involving state, input, and output. Based on the measurement obtained at time τ, the controller predicts the future moves of the system over a prediction horizon T and estimates the input that optimizes the predetermined open-loop performance objective function. To compensate for disturbances, model-patient mismatch and the finite nature of the optimization problem, only the first control action is implemented. The remaining samples are discarded and a new optimization problem is solved based on Y τ+1 at the next sampling step (τ + 1). The foregoing concept of MPC is analogous to the decision making process of a physician. When the patient arrives at the clinic for treatment, the physician diagnoses the patient, assesses the existing state of the disease, considers available treatment options, predicts the prognosis and administers the best available treatment to the patient but only until the next visit. When the patient returns for the next clinical visit, the physician repeats the same steps. Technically, MPC performs a similar task to what the physician does, perhaps in an optimal way, provided the model predictions are reasonable. While the traditional MPC algorithm is built with nominal dynamics of the model, we exploited the whole posterior distribution (with appropriate number of samples) of the individual patient parameters to account for the uncertainty.
Consider the equality constraint to maintain the drug concentration in the therapeutic level TLŷ i ðu i ðtÞ; x i ðtÞ; θ i ; ξ i Þ 2 fTLg; t ¼ 0; . . . ; T ð25Þ subject to the inequality constraints where T is the prediction horizon. g represent physiological constraints other than the clinical objective. For the case of regulating the system to the target concentration C, the quadratic cost function is defined as follows.
JðU;ŷÞ ¼ where u min and u max are minimum and maximum 6-MP doses allowed respectively. Δu min and Δu max signifies minimum and maximum allowed slew rate of 6-MP dose.

Off-Line Population Model Building
With the data sets D1 and D2, model parameters were estimated using the maximum likelihood approach as explained in Section 2.4.1. Feasible ranges of parameters were chosen based on experimental and clinical studies from literature. The statistics of the estimated parameters are given in Table 2. We assumed that the priors, p(θ, ξ), are log-normally distributed with this statistics. Using this prior distribution and individual patients' 6-TGN data D i , the posterior joint distribution of model and error parameters were estimated for each patient through the variational-Bayes approach outlined earlier. The posterior distributions for selected patients were verified with MCMC approach as this is a standard approach for performing Bayesian estimations. MCMC was done using Metropolis-Hastings algorithm implemented through R package 'mcmcpack' [44,45]. The posterior parameter distribution is essentially an updated form of the prior distribution in light of new information i.e. individual patient's 6-TGN concentration. Fig 3 shows marginal distribution of model and error parameters for a representative patient. Correlation coefficients for all parameters are acceptable except between k cm and k me . To examine the adequacy of the model in representing drug concentration data for various patients, we employed the global lack-of-fit test described by Blau et al. [26]. This test compares the occurrence of the experimental points within the highest probability density (HPD) regions for concentration predicted by the model. By definition, HPD is the region in which there is a 100(1 − α)% probability that the true value falls within the area under pðŷðtÞÞ which satisfies pðŷðtÞÞ ! pðyðtÞÞ. This confidence region (CR) is the smallest interval region among all credible intervals and hence is termed as the highest probability density region [46]. In our simulation, 259 out of 263 experimental points remained within the 95% HPD concentration confidence region, resulting in a confidence level for the lack-of-fit of 0.019. Since this value is less than 0.05, the selected model is adequate as measured by the global lack-of-fit test. Once the Bayesian parameter distribution is determined for all patients in the study, the population parameter distribution is formulated as explained in section 2.4.2. From each patient's posterior parameter distribution, a specified number (10000 Ã number of data points available for the corresponding patient) of samples are drawn. The final distribution for each parameter is approximated by a multivariate normal distribution. This approximate distribution is used as an informative population prior for all incoming new patients. In addition, every time a patient visits the clinic, this population prior is updated to include the new information. Fig 4 shows the population prior for model and error parameters. It is evident from the figure that most of the variability is explained by the reaction kinetic parameter k cm and, to some extent by the elimination rate of 6-TGN k me . Fig 5 shows the comparison of 95% concentration confidence region estimated with a non-informative prior (i.e. formulated with MLE) and informative prior. The red region shows the 95% HPD when parameters are estimated with non-informative prior. As expected the confidence region obtained from informative prior (gray region) is narrower. This exemplifies the advantage of the Bayesian approach where accumulation of information results in improved accuracy of prediction.

Model Reduction and Sampling Time Determination
Global sensitivity analysis was performed using the Sobol technique with 1000 set of parameters. Parameters were sampled through sparse grid sampling using statistics from the population prior distribution. Reaction and elimination rates of 6-TGN are the two parameters that emerged as highly sensitive ones. These two parameters have also displayed multimodal distributions in population priors (Fig 4). Similar observations were made in other population studies in literature [32]. Since the sensitivity indices varied with time, a cumulative error was calculated by assuming six representative time points across the treatment period. Accordingly, Eq (17) was modified as shown in Eq (30) with ε = 0.05. Table 3 shows the error associated with all the parameters. Any parameter with less than 2% of the error associated with the most sensitive parameter will be regarded as less sensitive and hence fixed at the population mean for all individual patients. The most sensitive parameter and hence the highest error involved was found to be k cm . Taking this as the reference error, the error involved was less than 2% for all parameters, except k me . However, as mentioned before, the correlation between k cm and k me is 0.96. Bayesian estimation with only k cm and k me as estimable parameters (other parameters were fixed at population mean) confirmed this trend and is shown in Fig 6. As a result, it suffices to estimate only k cm and fix all other parameters at the population mean for new patients. Consequently, the experimental design was formulated with the objective of improving the precision of parameter estimation for k cm . Fig 7 shows the evolution of Fisher's information as a function of time and parameter sets. From the figure, the maximum information is made available towards the steady state of the model. Concentration densities simulated with population prior also pointed that the maximum variation in the concentration distribution resulted when the drug concentrations are higher. However, it is not prudent to wait until the steady state to gather the data and identify the new patient. Hence, by compromising about 5% of the maximum information, 35 th day was determined as the optimal time to collect blood sample to measure the 6-TGN concentration.

Online Implementation
To this stage, we accumulated the information available in the literature by formulating a population prior distribution and showed the strategy for estimating the parameter for a new patient with a minimum number of samples. Since the optimal sampling time falls on the 35 th day from the beginning of treatment, it puts the patient in the dark for this period as there is no patient-specific model available for dose optimization. To overcome this drawback, we employed TPMT activity as an additional piece of information available up in the biological chain-of-response. For this purpose, data D3, consisting of TPMT measurement and corresponding 6-TGN concentration, are utilized to predict the optimum dose during the first 35 Table 3. List of parameters identified for deriving patient-specific model (bold face) together with other fixed parameters (regular face). Cumulative error, calculated at 1, 10, 20, 50, 75 and 100 days according to Eq 30, is given in column 2. % error is given in column 3. Individualized Dosing of Chemotherapeutics days. TPMT activity in these patients ranged from 7-to 30 U/ml/hr. With an activity window of 2 U/ml/hr. (i.e. Group 1: 7.01-9.00 U/ml/hr., Group 2: 9.01-11.00 U/ml/hr. etc.), 11 patient groups are formed based on their TPMT enzyme activity. Using the population prior distribution generated in section 3.1, a new set of parameter distributions were obtained for k cm for all patients in D3. With these distributions, group prior distributions were formulated as follows,

6-MP
where M G ðθ; ξÞ is the prior distribution for the group defined by the TPMT enzyme activity and N g is the number of patients in each activity group that ranged from 3-27. Eq. 31 is equivalent to Eq 16 in that Eq 31 considers the subset of patients having similar gene expression pattern as the population instead of the whole population. In summary, when a patient arrives at the clinic (at time t = 0), TPMT enzyme activity is measured and the patient will be placed in one of the 11 groups. Dosing decisions until t = 35 days are made with the assumption that the

Evaluation of Uncertainty
Another important objective of this work is to evaluate the impact of measuring various entities in biological chain-of-response on the uncertainty in predicting variables used for therapeutic drug monitoring and optimization. Despite the attractiveness of the DNA sequencing and gene expression profiling in characterizing patients, the more downstream variables, such as enzyme activity and active metabolite concentration, provide a much more robust indication of underlying response and thus helping to minimize the uncertainty. For this purpose, the individual parameter distributions are sampled to represent various measurements viz. whole population, groups based on TPMT genotype, groups based on TPMT enzyme activity and individual patient 6-TGN concentration measured at a single time point. Two TPMT genotype groups were formed based on enzyme activity. A TPMT enzyme activity of less than 10 U/ml/hr. is deemed as TPMT H /TPMT L (heterozygous) and an activity of more than 10 U/ml/hr. is deemed as TPMT H /TPMT H (homozygous-High) [31]. There are no TPMT L /TPMT L (homozygous-Low) patients found in this study. As Fig 8 shows, the black region is the 95% CR of concentration for the whole population included in this study. Without any attempt to collect genetic/ phenotypic measurements from a patient, the most that can be aspired for is that the concentration will fall in the range of 50-650 pmol/8x10 8 RBCs. This is indeed a reflection of the current status of clinical practice. Obviously, such uncertainties would preclude the possibility of dose individualization. The confidence regions for both the TPMT genotype groups were no different from this population CR. The gray region shows the 95% CR for the group having a TPMT activity of 15-17 U/ml/hr. With the measurement of a slightly downstream marker, the CR is narrower compared to the CR predicted for the population. However, when the 6-TGN concentration is measured and the model is adapted subsequently, the 95% CR is much narrower with notably lower uncertainty. With this patient-specific model, making accurate Comparison of 95% CR predicted using different information on biological chain-of-response for two representative patients. The black region represents the population. The gray region shows the prediction when only TPMT enzyme activity is measured. The red region shows the 95% CR predicted when just one measurement of 6-TGN is available (the solid dot). doi:10.1371/journal.pone.0133244.g008 Individualized Dosing of Chemotherapeutics prediction of drug concentration at hand, one can venture into dose optimization strategies to direct the drug concentration into a desired region that maximizes the efficacy and minimizes the side-effects.

Dose Optimization
Dose optimization is performed using the robust MPC strategy to achieve a therapeutically effective 6-TGN concentration. All dosing calculations are based on 15 days sampling horizon with 75 days treatment window as patients visit the clinic every two weeks. 50 parameter sets are sampled either from the group prior or patient-specific distribution. The combined error for all the parameters, weighted by their corresponding likelihood, was minimized within the MPC optimization. Clinical studies recommend a therapeutic 6-TGN concentration of 235-400 pmol/8x10 8 RBCs for an effective management of both efficacy and toxicity [47,48]. Hence, we optimized 6-MP input with a target 6-TGN concentration of 300 pmol/8x10 8 RBCs. Fig 9 shows the optimal 6-MP input together with resultant 6-TGN concentration for patients who had different response in relation to their respective groups. As mentioned earlier, without the patient-specific model until day 35, the dose is optimized based on the enzyme activity. The red region shows the 95% CR of concentration optimized for an enzyme activity group. When the patient-specific model is obtained with a measurement on the 35 th day, the optimized region shifts based on the drug concentration measured and reaches the target with narrow CR (shown in black). The CR in green shows the back calculation with the same dose as that of group optimum but with patient-specific model. The patient in subplot A had a lower reaction rate in relation to the group to which he/she belonged and hence when the actual measurement becomes available the dose had to be increased to push the 6-TGN concentration higher. The opposite is true for the patient in subplot B. Dose inputs for different patients suggest that the dosage varied as much as 200% and as low as 25% of the standard dose, which is not uncommon in the clinical practice. Although there are obvious and significant differences between the standard and optimized 6-MP usage, the merit of dose optimization should be viewed from the maximization of therapeutic benefits rather than the reduction of drug input as the cost of drug is only a fraction of the overall healthcare spending.

Discussion
Individualized treatment strategies are increasingly favored and extensively recommended for many drugs, especially for deadly diseases like cancer. Availability of modern analytical facilities and computing power for in-silico approaches support recommendations to further such efforts. Implementation of these advances in clinics faces several challenges. Among these, the most important issues are the technical ones which are mainly concerned with robust determination of optimal dosing for an individual patient. In this work, we have engineered a strategy for a model-based individualized treatment for an important chemotherapeutic drug. We have devised a Bayesian approach to describe the population characteristics and subsequently utilized it for robust estimation of patient-specific parameters in light of sparse clinical data. The Bayesian approach is inherently suitable to consolidate the information from various sources in various formats into a prior distribution. Highly informative priors play a significant role in dealing with situations where real-time decisions have to be made with partial information i.e. sparse data. The use of distributions of parameters arising from the Bayesian approach, instead of parameter point estimates, enabled us to perform a robust estimation of sampling times, to optimize the dose and to quantify the associated uncertainty. More importantly, we have shown the advantage of measuring the downstream biomolecular phenotypic indicators, rather than an upstream genotypic marker, in minimizing the uncertainty in predicting clinical variables of our interest. While the latest DNA microarray technologies have made genotyping an inexpensive way to characterize the patients, its utility in individualized dosing comes at the price of significant uncertainty. In section 1.3, we described how important it is to measure the downstream clinical variables for accurate prediction of dose response. Although 6-TGN concentration itself was shown to be a valuable indicator of clinical efficacy and toxicity and proven useful in certain clinical conditions where efficacy/toxicity measures are categorical and highly subjective, extension to include cellular response will certainly add value to the approach. Additional variations while the active drug imparts cytotoxicity on various cell populations may play a role and impose another level of uncertainty. Hence, an interesting and important extension to this work would be to connect this model to relevant pharmacodynamic variables (cellular response). With regard to 6-MP treatment, these cellular responses primarily involve the bone marrow cell population. Works are in progress to consider these extensions.
Besides the issue of 'technical know-how', clinical implementation is also constrained by physiological, logistical, economic and social factors. Physiological issues are concerned with whether the technical solutions are realizable within the physiological constraints. For example, if the approach demands several additional blood samples, it will be prohibitive to implement in pediatric patients, no matter how beneficial it is. Certain types of measurements, such as those from the bone marrow are highly restricted. Logistical issues surface mainly during the translational phase. They arise due to incompatible resources at the healthcare facilities. For example, if the algorithm requires a measurement 6 hours post-dose, managing patients arriving at different times in the clinics (inpatient or outpatient) would be an issue. In addition, timely analysis of samples and reporting the results would be important consideration. The economic issues are obviously concerned with the additional cost involved for procedures necessitated by the individualized approach. Individualized treatment means clinically identifying or characterizing an individual patient through genotyping and/or phenotyping that invariably requires additional testing/ lab assays. Basically, one has to show the potential benefits, both short-and long-term, weighed over the cost, to convince third party payers to cover the tests. The final and important social issues encompass communication, information sharing, education to healthcare providers and patients on the new approach, patient compliance, policy, ethical and privacy issues. It has been our aim in this work to address some of the technical issues within the physiological constraints. We believe that by comprehensively addressing the foregoing two issues and corroborating the potential incentives will pave way for overcoming the rest. This will ensure that the patients are treated with the suitable dose, dictated by their own genetic/phenotypic background, ultimately resulting in improved qualityof-life.