Adaptive Estimation of Personalized Maximum Tolerated Dose in Cancer Phase I Clinical Trials Based on All Toxicities and Individual Genomic Profile

Background Many biomarkers have been shown to be associated with the efficacy of cancer therapy. Estimation of personalized maximum tolerated doses (pMTDs) is a critical step toward personalized medicine, which aims to maximize the therapeutic effect of a treatment for individual patients. In this study, we have established a Bayesian adaptive Phase I design which can estimate pMTDs by utilizing patient biomarkers that can predict susceptibility to specific adverse events and response as covariates. Methods Based on a cutting-edge cancer Phase I clinical trial design called escalation with overdose control using normalized equivalent toxicity score (EWOC-NETS), which fully utilizes all toxicities, we propose new models to incorporate patient biomarker information in the estimation of pMTDs for novel cancer therapeutic agents. The methodology is fully elaborated and the design operating characteristics are evaluated with extensive simulations. Results Simulation studies demonstrate that the utilization of biomarkers in EWOC-NETS can estimate pMTDs while maintaining the original merits of this Phase I trial design, such as ethical constraint of overdose control and full utilization of all toxicity information, to improve the accuracy and efficiency of the pMTD estimation. Conclusions Our novel cancer Phase I designs with inclusion of covariate(s) in the EWOC-NETS model are useful to estimate a personalized MTD and have substantial potential to improve the therapeutic effect of drug treatment.


Introduction
It is common for a group of patients with the same cancer type to receive the same treatment. However, some patients will experience substantially better therapeutic effects than others, and some anticancer therapies may benefit only a subset of treated patients. Several reasons account for the heterogeneous therapeutic effect observed at the same dose level of the same drug. Patients have different genetic and environmental profiles, including demographic characteristics, concomitant diseases, concomitant drugs, biomarkers, SNP copy number, etc. [1,2,3]. Genetic and environmental factors interactively affect the therapeutic effect of a treatment intervention. Tumor heterogeneity is another significant reason for the heterogeneity of the toxicity and therapeutic effects of a drug. Tumors of a primary site in many cases represent a heterogeneous collection of diseases that differ with regard to the mutations that cause them and drive their invasion, thus are heterogeneous with regard to natural history and response to treatment. Personalized medicine has evolved recently as an advanced approach to achieve optimal medical effect in the context of a patient's genetic, environmental, and tumor profiles [1,2,3,4,5].
The first critical step toward personalized medicine is the estimation of personalized maximum tolerated dose (MTD) in a Phase I clinical trial, which is the first trial of a new drug in humans after animal studies, with the main purpose to determine the MTD of a new drug under safe administration. In a Phase I clinical trial, there is substantial heterogeneity in dose limiting toxicity (DLT) response at the same dose level of the same drug among different patients because of different genetic and environmental profiles and tumor heterogeneity. Some known factors include the vulnerability to an exaggerated pharmaco-dynamic effect (potentially mediated by receptor differences), differences in genetic susceptibility (e.g. biomarker, G6PD deficiency), and drug-drug interactions [6,7]. Ignoring the potential heterogeneity may lead to serious bias in MTD estimation for different groups of patients [7,8], and as a result, the therapeutic effect is substantially decreased. Hence, in order to achieve the optimal therapeutic effect of a drug for every patient, estimating a personalized MTD offers greater potential than estimating an overall MTD across different patients [8].
The main goal of this study is to develop a practical and leading Phase I design that can facilitate the estimation of personalized MTD for the implementation of personalized medicine. Currently available Phase I designs can be classified as rule-based or model-based. Rulebased (or non-parametric) Phase I designs fail to estimate MTDs while adjusting for covariates due to their simple up and down algorithms. Therefore, a parametric or semi-parametric model-based design is desired so that covariates, especially genomic profiles, can be included into the dose response curve. Among several parametric Phase I designs available in the literature, Escalation With Overdose Control (EWOC), proposed by Babb, et al. [9], can control the probability of exceeding the MTD during the dose escalation phase and has been used in trials at Emory University, Fox Chase Cancer Center, Miami University, Novartis, and other institutions. EWOC can detect the true MTD with high accuracy, compared with traditional 3+3 designs. Nevertheless, EWOC only considers the worst toxicity event that a patient experiences. A binary outcome is used to denote whether the worst toxicity event that occurs has DLT status. Therefore, EWOC design is limited by its binary outcome when investigators want to consider not only DLT events but also non-DLT events. In a real clinical trial, patients often experience multiple toxicities with event grades varying from 0 (no toxicity of that type) to 5 (death). Some patients may experience multiple DLTs. In addition, the toxicity events are not equally severe. For example, a grade 4 renal toxicity of severe and possibly irreversible may be more severe than a grade 3 reversible toxicity of another type, such as grade 3 neutropenia. Moreover, some, but not all, low-grade non-DLT events that occur at a lower dose may be a sign of the occurrence of the same toxicity of a higher grade at a higher dose. Different toxicity events often occur together, such as fatigue and nausea/vomiting, myelosuppression and fever, etc. Considering both DLT and non-DLT events will allow the more accurate evaluation of the toxicity of a new drug [10,11,12,13].
We previously proposed a novel toxicity scoring system to overcome the above limitations by incorporating multiple toxicity responses as a quasi-continuous variable and thus, fully utilize all toxicity information [14]. The novel toxicity scoring system has two major components: Normalized Equivalent Toxicity Score (NETS) and Target NETS (TNETS). It has been successfully combined with the non-parametric Isotonic Design, to create a new semi-parametric design named ID-NETS. In Chen et al. [15], we further proposed to replace the binary outcome with NETS in the EWOC design so that the extended method (EWOC-NETS) can be used to evaluate the MTD by considering both DLT and non-DLT events. In this study, in order to estimate a personalized MTD, we propose to include baseline covariates that represent personal tumor, genetic, and environmental factors into the EWOC-NETS model. The selected baseline covariates must be biomarkers or classifiers with strong biological evidence that they may have some kind of linear relationship with the dose and drug effect, or that their categorical status can successfully predict the possibility of patient benefit from the agent. The varying MTD estimated in the novel method is the personalized MTD according to a patient's tumor, genomic and environmental profile, instead of an overall MTD. In the remainder of this manuscript, we will describe in Section 2 details of the NETS system, EWOC-NETS, and the new model including the covariate. In Section 3, extensive simulations will be presented to examine the performance of EWOC-NETS after considering the covariate. In Section 4, we will discuss the limitations and applications of the new method.

Major differences between EWOC and EWOC-NETS
EWOC-NETS was developed based on the framework of EWOC with two new concepts. The first is to fully utilize all toxicity information by replacing the binary outcome with the NETS (ranging from 0 to 1) which can be interpreted as the severity of a patient's toxicity response. A quasi-likelihood function is further introduced to update the posterior marginal distribution of the MTD (γ). The second concept is that the MTD is re-defined as a dose corresponding to a pre-specified TNETS, θ, ranging from 0 to 1 instead of a Target Tolerated Level (TTL). The new definition satisfies where ANETS is the Average NETS of the dose level. TNETS, θ, is an analog to the Target Toxicity Level (TTL) when using a binary DLT in a traditional clinical trial design and is determined by the Target Toxicity Profile, which is a more detailed toxicity response of a patient treated at MTD γ. EWOC-NETS assumes that the dose ANETS relationship follows a logistic model and NETS can also be interpreted as the observed fractional events ranging from 0 to 1 (see Fig 1 of Chen et al. [12]).

Quantitative measurement of toxicity response with NETS
Chen et al. [14] provide a detailed description of the method to derive NETS and determine TNETS for a clinical trial. Suppose in a trial with K patients, the i th patient has a total of N i toxicities. Toxicity events are assigned to adjusted grade according to 1 for grade 1 toxicity, 2 for grade 2 toxicity, 3 for grade 3 non-DLT, 4 for grade 4 non-DLT, 5 for grade 3 DLT, and 6 for grade 4 DLT (see Table 1 of Chen et al. [12]). Let G i,j be the adjusted grade of the j th toxicity (1 j N i ) of the i th patient. The maximum adjusted grade, G i,max , of the i th patient can be obtained by Eq (2): Let S i be the NETS of the i th patient. If the i th patient has no toxicity or only grade 1 toxicity (adjusted grade = 1), the S i is arbitrarily assigned to 0 or 1/60, respectively. Otherwise, S i can be derived by Eq (3).
The parameter r i,j ranging from 0 to 1 is a weight for the correlation of the j th toxicity with other toxicities of the i th patient. The weight decreases as the correlation increases. According to Chen et al. [14], the parameters G i,max and c are suggested to be fixed at 6 and -2, respectively. The parameter β is the slope for the increasing rate of NETS controlling for toxicity events. The higher the value of β that is chosen, the more conservative the clinician is. β is recommended to be a value ranging from 0.1 to 0.5, such as 0.25.

Determination of TNETS
TNETS, θ, is an analog to the Target Toxicity Level (TTL) when using a binary DLT in a traditional clinical trial design. It is determined by the Target Toxicity Profile, which is a more detailed toxicity response of a patient treated at MTD. The calculation of TNETS, θ, is defined by Here, p l is the probability that the worst toxicity event is an adjusted grade l event and m l is the corresponding mid-range NETS value when the toxicity with an adjusted grade of l is the "worst" toxicity for a patient with varying additional less severe toxicities in the target toxicity profile and the maximum adjusted grade for the patient is l. The range of the parameter l is from 0 to 6. In the above equation, m 0 = 0, m 1 = 0.092, m 2 = 0.25, m 3 = 0.417, m 4 = 0.583, m 5 = 0.75, and m 6 = 0.917.

EWOC-NETS model including covariates
In an ideal situation, we need to fully consider all possible covariates and their interaction effects in order to estimate a personalized MTD. The full model of the dose and toxicity response relationship should include a constant, a vector of covariates including dosage, patients' characteristics, other biomarkers, and the overall interaction terms as below. ; However, the small sample size of a Phase I trial makes it inadvisable to use an overly complicated model with too many covariates. Therefore, in the models proposed in this study, we limit to one covariate that is either discrete or continuous and have omitted the interaction term.
2.4.1 Discrete covariate. Let X min and X max denote the minimum and maximum dose levels pre-specified by the clinician. Suppose a discrete covariate C with a value 1 means that the patient is in group A and a value 0 means that the patient is in group B. According to the preclinical trial, we know that group B generally has higher MTDs in males than in females. Let γ max denote the MTD for group B and γ 0 denote the MTD for group A. Therefore, the relationship would be satisfied. The dose assigned to the first patient is X min and we shall select only dose levels between the interval of X min and X max. Let X i denote the dose assigned to the i th patient, i = 1, 2, . . ., K, then x 1 = X min and x i 2 [X min ,X max ], 8i = 1,. . .,k. We model the relationship between dose and ANETS by where F is a pre-specified distribution function, called a tolerance distribution, and β 0 and β 1 are unknown. We assume that β 1 > 0 and δ < 0 so that the ANETS monotonically increases when the dose level increases, adjusting for the covariate, c i . The MTD is the dose level, denoted by γ, such that the TNETS is θ. It follows from Eq 6 that where F is a logistic regression model here. We can further simplify Eq (7) to obtain Eq (8).
Let ρ 1 denote the ANETS at the starting dose x 1 = X min for group A and ρ 2 the ANETS at the starting dose x 1 = X min for group B.

Continuous covariate.
The basic method to include a continuous covariate into the EWOC-NETS model is the same as the method to include a discrete covariate. However, the re-parameterization in Eq (13) is different. Let Z denote a continuous variable where the MTD increases when z increases, conditional on the same dose level. Then Eqs (8)-(10) should be changed by the new dose response relationship. Let z min denote the smallest value of Z and z max be the maximum of Z in the Phase I clinical trial. After considering age, the dose-response curve satisfies A different re-parameterization for δ and β 0 is derived by solving Eqs (14)- (16).

Quasi-Bernoulli likelihood
Frequentist quasi-likelihood methods are designed to model overdispersion observed in binomial or Poisson data. When the "quasi" distributions belong to linear exponential families such as the binomial family, Quasi Maximum Likelihood Estimates (QMLEs) obtained by maximizing the quasi-Bernoulli likelihood function are strongly consistent [10,16,17,18].
Recently, the quasi-likelihood approach has been successfully combined with Bayesian generalized linear models [18] and the continual reassessment method (CRM) [10]. NETS can be viewed as fractional events. So, we assume that the variance structure of NETS (S) is m S i jX i ð1 À m S i jX i Þ conditional on dose level X i for the i th patient where S 1,. . ., S N are assumed to be mutually independent in a clinical trial with N patients. The quasi-Bernoulli likelihood can be applied to update the posterior distribution of (γ max , ρ 1 , ρ 2 ) under these assumptions.
If we want to include a continuous variable in our model, the data after observation of k patients can be expressed as D k = {(x i , z i , s i ), i = 1, 2, . . ., k}. D k would include the dose assigned (x i ), the NETS observed (s i ), and the continuous covariate (z i ) of each previously treated patient. The quasi-Bernoulli likelihood of (γ max , ρ 1 , ρ 2 ) given D k is Obviously, the true distribution in NETS is unknown so we try to make inferences based on the known variance structure and the corresponding quasi-likelihood. According to the quasi-Bayesian theory, the quasi-likelihood function can be interpreted as the "limited information likelihood" and it is the best approximation of the true likelihood. After plugging in the corresponding dose ANETS curve, the likelihood function can be written as

Overdose control in EWOC-NETS
The traditional Bayesian decision theory makes inferences of the posterior mean, median or mode so that the corresponding expected loss function can be minimized. In order to control the overdose rate, we do not choose the posterior median of EWOC-NETS as an estimator of MTD at the beginning of the clinical trial. Instead, the α th percentile of the posterior marginal distribution is chosen. As a result, the marginal posterior overdosing probability is equal to α for the next patient. α was referred to as the feasibility bound by Babb, et al [9]. Here, we generally present how to extend overdose control after including a continuous covariate Z. The posterior distribution after the observation of k patients can be derived by where p(γ max , ρ 1 , ρ 2 ) is the joint prior of the new parameter (γ max , ρ 1 , ρ 2 ). We assume that the priors for γ max , ρ 1 , ρ 2 are mutually independent and follow uniform distributions. The posterior joint distribution of (γ max , ρ 1 , ρ 2 ) is dominated by the observed data. The MTD of the (k+1)th patient adjusting for z k+1 can be re-parameterized as The posterior distribution of g z kþ1 can be updated by the joint posterior distribution of (γ max , ρ 1 , ρ 2 ). We denote the cumulative density function (CDF) for g z kþ1 as p k;z kþ1 ðÁjD k Þ: p k;z kþ1 ðÁjD k Þ. It is not only related to the cumulative data D k, but also related to the next patient's covariate z k+1 . In order to ensure the ethical constraint of overdose control, the selected dose x k+1 for the new patient must satisfy equation: At the end of the trial, the personalized MTD is the posterior median of the corresponding posterior marginal distribution adjusting for the covariate.
During the implementation of Bayesian procedures, the posterior distributions are estimated using the Markov Chain Monte Carlo (MCMC) method. The MCMC method has been widely used in Bayesian frameworks to sample posterior distributions with high dimensional parameters [19]. The Metropolis-Hastings algorithm is used to obtain a sequence of random samples). In the EWOC-NETS design, the burn-in period is 1000 iterations with another 1000 iterations conducted to sample the posterior distribution. Tighiouart et al. [20] reported a successful example of using MCMC to study a large class of prior distributions in EWOC.

Additional advantages of EWOC-NETS
Simulation studies and their application to real clinical trial data demonstrate that EWOC-NETS maintains the advantages of EWOC. In addition, it provides new advantages: 1) treats toxicity response as a quasi-continuous variable; 2) improves the MTD accuracy by differentiating beyond the DLT; 3) increases trial efficiency by fully utilizing all toxicities [15].

Models incorporating patient characteristics and biomarkers
In this manuscript, we use EWOC-NETS as a framework to incorporate patient characteristics and biomarkers. The original EWOC-NETS without considering any covariates is used as a baseline model for comparison. Therefore, a total of 3 models is considered: 1) Model 1 is the original EWOC-NETS model not including any covariates (baseline model); 2) Model 2 is the EWOC-NETS model considering only a binary covariate C (0 or 1); 3) Model 3 is the EWOC--NETS model considering only a continuous covariate Z (0~1). The logistic model of the relationship between NETS and dose x for each of the 3 models is summarized below: The corresponding parameters are re-parameterized as (γ max , ρ 1 , ρ 2 ) as mentioned in Section 2. In order to evaluate the additional advantages in the performance of EWOC-NETS after considering covariates, model 2 and model 3 are compared to the baseline model 1 using bias, standard error (SE) and mean square error (MSE), respectively, under different scenarios.

Simulation settings
In order to compare the performance of the extended methods and the original method, models 2 and 3 are evaluated mainly under four scenarios each. The γ max is defined as the true MTD for the group with C = 1 and γ 0 is defined as the true MTD for the group with C = 0 when the covariate is discrete or continuous. The simulation set up under the 8 scenarios (S1 to S8) is summarized in Table 1. The first 4 scenarios (S1 to S4) are for discrete covariates and the other 4 scenarios (S5 to S8) are for continuous covariates. The covariates considered actually have true effects on MTDs under S1, S2, S3 for discrete covariates and under S5, S6, S7 for for continuous covariates, respectively, so that γ max and γ 0 are both set to 0.5 under both scenarios. In order to obtain comparable results, the dose X is standardized within the range from 0 to 1; the continuous covariate Z is also standardized within the range from 0 to 1; and the binary covariate C is either 0 or 1. The ANETS is generated by the true tolerated function. The NETS is generated by a truncated normal distribution with the mean, ANETS, and the assumed variance structure. Under S1, S2, S3 and S4, we enroll 15 cohorts with covariates of value 1 (C = 1) and 15 cohorts with covariate of value 0 (C = 0). The group with C = 1 is treated first because we assume that their MTD is higher and it is safer to treat them first. Under S5 to S8, the continuous covariate values are generated from a uniform distribution (0, 1). The comparisons of the performances of the different models are based on: 1) whether the estimated MTD is a personalized MTD; 2) the amount of bias of the final estimation ofĝ max , defined as P n i¼1 ðĝ max À g max Þ=n, where n is the total number of simulations; 3) the mean square error (MSE) of the final estimation ofĝ max where a smaller MSE yields better performance of the models; and 4) the standard error (SE) of the estimatorĝ max where the lower the standard error the more stable is the estimator.
In each simulated trial, the TNETS level is set to 0.476, which is calculated based on a target toxicity level of 33% DLT with a target toxicity profile consisting of equal probability for each non-DLT or DLT toxicity, respectively. The sample size is set to 30 in each simulated trial. The feasibility bound, α, is set to start at 0.25 for every model. α is increased with an increment unit of 0.05 when we assign the dose to the next cohort. The maximum of the feasibility bound is 0.5. The trial starts with the lowest dose level and the recommended dose level for the next cohort is the α th percentile of the posterior marginal distribution of the MTD adjusting for its covariate. Each scenario is simulated 1,000 times. We define the limiting NETS status (LNETS) as 1 if the observed NETS exceeds 0.476, and 0 otherwise. Meanwhile, since we choose a dose level based on a continuous scale, we assume a tolerance, 0.05, for both overdosing and limiting NETS status.
γ max , ρ 1 and ρ 2 are respectively assumed to follow mutually independent priors: uniform distribution (0, 1), uniform distribution (0, 0.476), and uniform distribution (0, 0.476), respectively. The posterior sample of (γ max , ρ 1 , ρ 2 ) is directly sampled by the Metropolis-Hastings algorithm implemented by JAGS, an efficient software for sampling the MCMC chain. The burn-in period is 1000 iterations. We use an additional 1000 iterations as the posterior sample of (γ max , ρ 1 , ρ 2 ). Trace plots and histograms are used to diagnose whether the MCMC chain has converged. The histogram and trace plot for the parameter of primary interest, γ max have shown that after 2000 iterations, the MCMC chain is stable and becomes a unimodal curve so that the posterior sample can provide more evidence to infer the parameter of interest γ max by the posterior median.

Inclusion of a discrete covariate
Comparison of the accuracy of the estimated parameters between model 1 and 2 is summarized in Tables 1 and 2. The advantage of model 2 is obvious when there is a huge difference between two groups (C = 0 and C = 1) in S1. Model 2 successfully detects a different MTD for the two groups when the true MTD for the C = 0 group increases from 0.27 to 0.5. In contrast, model 1 always estimates only one MTD for both groups although these two groups should have two different MTDs (Table 1 and 2). The bias for the MTD estimation using model 2 is the smallest in S1, respectively, 0.058 and -0.013 ( Table 2). The standard error and MSE are also the smallest in this scenario ( Table 2). The bias, standard error and MSE increase when the difference between the MTD of the two groups decreases from S1 to S4 ( Table 2). The results show that the MTD estimation for the two groups using model 2 will be more precise when the covariate effect on ANETS is larger (Table 2). Under the worst scenario S4, model 2 estimates γ max as 0.6 and γ 0 as 0.487, both of which are close to their true values (0.5), suggesting the robustness of model 2 when the covariate considered has no true effect. On the other hand, the bias, standard error and MSE for model 1, which only estimates a mean MTD for two groups, tends to decrease as the covariate effect on the toxicity outcome in terms of ANETS decreases from S1 to S4 (Table 1). Therefore, the MTD estimation for model 1 will be more precise when the covariate effect on ANETS diminishes (Table 1). More patients are required when the covariate effect on ANETS is small.
Compared with model 2, MTD estimation is a pooled estimate of the two groups of patients. As a result, the bias of the pooled estimator for the two groups is larger when the covariate effect on ANETS increases. Because of the small sample size in each group in the trial, the standard error of the MTD estimation is larger in model 2. The patient distributions among different dose levels under model 1 and model 2 are shown by the box plots in Fig 1  and Fig 2, respectively. Patients in different groups are treated at dose levels concentrated around their personalized MTD under model 2 (Fig 2). By contrast, patients are concentrated around the pooled estimation of MTD under model 1 (Fig 1). Since a dose level near the MTD is more effective and safer, the therapeutic effect for patients is better under model 2 than under model 1. The comparison of operating characteristics between model 1 and model 2 in terms of the overdosing rate, as measured by MTD+ percentage and LNETS percentage, is summarized in Table 3. The rates for patients being treated at dose levels higher than the MTDs (MTD+%) in all scenarios are smaller under model 2 than under model 1. Similarly, the rates for dosing above the limiting NETS status (LNETS+%) are lower under model 2 than model 1 among all scenarios. The difference between the overdosing rates of model 2 and model 1 increases as the difference in MTDs of the two groups increases when stratified by a discrete covariate. For example, the LNETS+ rate in scenario S1 is decreased by 22% under model 2 (39%) compared to the rate under model 1 (61.6%). The inclusion of the discrete covariate not only makes the final MTD recommendation more precise but also increases the therapeutic effect, and is also more ethical for patients.

Inclusion of a continuous covariate
Comparison of the accuracy of the estimated parameters between model 3 and model 1 under the 4 scenarios (S5 to S8) is summarized in Tables 1 and 2. Under the ideal scenarios (S5, S6, and S7) in which the continuous covariate considered has a true effect on the MTD, model 3 can estimate a personalized MTD depending on the value of the continuous covariate of a  (Table 2). However, the bias, standard error, and MSE for model 1 decrease when the true effect of the continuous covariate decreases from S5 to S8 (Table 2). Under the worst scenario S8, model 3 estimates γ max as 0.593 and γ 0 as 0.468, both of which are close to their true values (0.5) ( Table 2). This demonstrates the robustness of model 3 under all scenarios. Adaptive Estimation of pMTD with EWOC-NETS Comparisons of the overdosing rates between model 3 and model 1 are summarized in Table 3. Under scenario S5, consideration of a continuous covariate in EWOC-NETS (model 3) can reduce the MTD+ rate from 37% to 30% and the LNETS rate from 48% to 43%, resulting in better therapeutic effect for the participating patients in the Phase I trial (Table 3). When the considered continuous covariate actually has no true effect in MTD (S8), the MTD+ rate and LNETS rate of model 3 are still lower than the corresponding rates under model 1. This suggests that we will lose little when we consider a continuous covariate without a true effect on the MTD.

Sample size and number of covariates
From Sections 3.3 and 3.4, we can see that covariates should be included in the model when they have a true effect and that considering non-effective covariates has little loss under the worst scenario. The sample size of each covariate value is important when we consider the inclusion of the covariate because a small sample size of each covariate value makes the posterior estimation unstable in 1000 simulations. In order to estimate the MTD more precisely, more patients with the same covariate value are needed in a real clinical trial.

Applications
The above simulations are hypothetical/computer-generated in order to demonstrate the operating characteristics of the novel designs. The designs could be applied to available Phase I data/covariates to confirm their utility. For example, using the designs, three Phase I clinical trials have been designed to incorporate patient characteristics and biomarkers to estimate pMTDs for head and neck cancer patients at the Winship Cancer Institute of Emory University: 1) To estimate pMTDs for panitumumab according to each patient's binary status of human papilloma virus (HPV) (positive vs negative); 2) To estimate pMTDs for folate-dextran-paclitaxel (FDT) according to each patient's folate receptor level (continuous); 3) To estimate pMTDs for luteolin according to each patient's specific type of head and neck cancer (categorical).

Discussion
The concept of personalized medicine was introduced as early as the 1970s. Byar et al. originally proposed to select optimal treatment in clinical trials using covariate information [21]. Substantial advances have been made in this century due to progress in modern sequencing technologies, the development of new therapeutics, and the significant contribution of clinical trial methodologies. Personalized medicine has become a crucial component of contemporary cancer medicine. Examples of successful personalized medicines include several that target the epidermal growth factor receptor (EGFR). Imatinib (Gleevec), a tyrosine-kinase inhibitor of EGFR, was initially invented in the late 1990s by a biochemist, Nicholas Lyndon, and then approved by the FDA in 2001 for the treatment of multiple cancers, most notably Philadelphia chromosome-positive (Ph + ) chronic myelogenous leukemia (CML) [22]. Imatinib has benefitted thousands of patients with CML and gastrointestinal stromal tumors (GIST). Cetuximab, an EGFR antibody, has been approved by the FDA for the treatment of colon cancer with wild-type KRAS, but not those with a KRAS mutation [23], and for the treatment of certain stages of head and neck cancer. Panitumumab, a fully human monoclonal antibody specific to the EGFR, has been approved by the FDA for the treatment of EGFR-expressing metastatic colorectal cancer with disease progression [24]. In 2012, the FDA approved a real time PCR companion diagnostic test for KRAS, Therascreen KRAS test, which is the first genetic test to guide the treatment of cancer. Recent advances in biotechnology have resulted in a shift toward molecularly targeted anticancer agents, which are likely to benefit only a subset of the patients with a given cancer. Due to the molecular heterogeneity of most human cancers, only a subset of treated patients benefit from a given therapy. This is particularly relevant for the new generation of anticancer agents that target specific molecular pathways. When biomarkers to identify patients who are likely to benefit from a given therapy are available, targeted clinical trials that restrict eligibility to sensitive patients can be conducted. The identification of the appropriate "sensitive" population requires definitive testing of specific biomarkers.
Ideally, such diagnostic tests should be developed and validated before designing the trial. New biotechnologies such as microarrays can be used as powerful tools to measure biomarkers or a genetic signature and identify patients that are most likely to benefit from anticancer therapies. For example, Derin et al. reported that a lower level of MAPK expression is associated with anthracycline resistance and decreased survival in patients with hormone receptor negative breast cancer [25]. Fan et al. summarized the concordance among gene expression-based predictors for breast cancer [26]. Rapid advances in biotechnology have made it possible to obtain quantitative information regarding biomarkers and differentiate patient subsets according to sensitivity to agents. For example, EGFR is a target for molecularly targeted agents in lung cancer therapy. EGFR expression is first measured immunohistochemically on a continuous or graded scale, and is later used to categorize patients into several distinct categories for clinical management or dichotomize with various cut points, such as in the DAKO kit [27]. Personalized medicine has thus become a significant and valuable new approach in the field of medicine and the estimation of a personalized MTD according to a patient's environmental and genetic profile is a critical step toward personalizing their treatment. The goal of this study is to propose an extended Phase I design which can fully utilize all toxicity information and patient's characteristics to estimate personalized MTDs for precision medicine. Through simulation studies, we demonstrate the advantages and potential loss of estimating the conditional MTD given covariates using the EWOC-NETS design proposed by Chen et al. [15]. From the simulations, we can conclude that the extended EWOC-NETS incorporating patient characteristics and biomarkers can: 1) estimate a personalized MTD; 2) reduce the probability of patients being overdosed; 3) increase the accuracy of the treated dose, thus improving the therapeutic effect; and 4) have little loss when covariates being considered have no effect. Some [31]. Our method is based on EWOC-NETS which is an extension of EWOC that uses NETS instead of the probability of DLT. The comparison of performance between EWOC and EWOC-NETS has been elaborated in previous publications [14,15,32,33]. Our results based on EWOC-NETS are consistent with the findings of O'Quigley et al. [28,29], Babb et al [30], and Tighiouart et al [31]. Besides the Phase I clinical trial methodologies developed to determine personalized MTD, other novel approaches have been proposed to use biomarkers or genomic signatures for personalized medicine in Phase II/III clinical trials in oncology and other diseases [1,2,3,4,5]. For example, Freidlin and Simon proposed an adaptive signature design to generate and prospectively test a gene expression signature for sensitive patients [5]. Jiang et al. further developed a biomarkeradaptive threshold design for evaluating treatment with a possible biomarker-defined subset effect [1]. Freidlin et al. recently extended their adaptive signature design and developed a cross validated adaptive signature design which has considerable improvement in performance [2].
Our extended EWOC-NETS design can be applied to estimate the personalized MTD of a new agent according to the status of specific biomarkers of the targeted cancer. When planning a Phase I trial, we need to decide whether or not to include patient covariates and which covariates to include in the MTD estimation during the trial. Some well-known biomarkers and their related diseases are summarized in Table 4. Each of these biomarkers as well as other patient characteristics can be treated as binary or continuous covariates and incorporated into the corresponding models. The Phase I clinical trial is typically a small study with a small sample size and as such we cannot attempt to estimate too many parameters. Therefore, when multiple biomarkers need to be considered, we can estimate their combined effect as a single score by using an additive linear model and incorporating the score into the model as a continuous or binary variable. We should not consider any covariates that are known to have no effect on the MTD. From a statistical point of view, we stand to lose little in terms of the accuracy of the MTD and therapeutic effect for participants when covariates taken into account in the model are actually not predictive of the severity of patients' toxicity response. However, it is a substantial monetary cost to quantify covariates, especially when patients need to be genotyped and certain biomarker expressions need to be determined. Therefore, careful consideration must be given to the balance of these pros and cons.
Subsequent Phase II and III trials should further test the efficacy of personalized MTD among biomarker-stratified patients instead of testing an average effect among an unselected population. In an enrichment design, all patients will be homogeneous regarding a specific binary biomarker, so that the dose recommended to all participants will be the personalized MTD according to the status of the biomarker. In a biomarker-stratified design, patients of each stratum will be treated with a specific pMTD according to the level of the categorical biomarker the patient has in the stratum. In more complicated settings, where no natural cutpoint of the biomarker is known in advance, or there are multiple biomarkers and high dimensional genomic tumor characterization, patients will be treated with pMTD according to the value of a continuous biomarker or a composite score summarizing their overall effect.
There are also some limitations in the real practice of the design. For example, the required genomic or demographic data from incoming patients may need to be retrieved from different sources and be integrated together for pMTD determination. Different data are not always available immediately and sometimes may be inconsistent. Therefore, reliability, validation, and standardization of data from different sources are critical during the data processing, storage, analysis, and interpretation to make real-time personalized assessments. Most previous studies have been conducted based on an individual study level. Patient privacy during the sharing of data is another consideration. Fortunately, increasing computation capability, rapidly expanding databases, and powerful safe networks enable the immediate retrieval and integration of bioinformatics and clinical data to deliver real-time information for pMTD estimation and generate a highly personalized outcome while protecting patient privacy. Further database development is highly recommended for the estimation of pMTD and implementation of personalized medicine, but the application of our Phase I design for pMTD should not be limited by the availability of related databases. Biomarkers have played significant roles in the adaptive randomization, group stratification, and patient enrichment of Phase II and III clinical trials. The pMTD estimated from a Phase I trial should be utilized to treat different groups of patients with personalized doses according to their genomic and clinical characteristics, thus optimizing the efficacy of precision medicine and maximizing therapeutic effects. Given the potential benefits of personalized therapy, the estimation of pMTD as the first step toward this goal should be pursued whenever at least some genomic and clinical data are available.
In summary, our study has shown that the inclusion of covariate(s) in the EWOC-NETS model is useful to estimate a personalized MTD and has substantial potential to improve the therapeutic effect of drug treatment.