Marie-Karelle Riviere

Ying Yuan

Frédéric Dubois

Sarah Zohar

A Bayesian dose-finding design for drug combination clinical trials based on the logistic model

Keywords: Bayesian inference, Dose-finding, Drug combination, Oncology, Phase I trial. 1

In early phase dose-finding cancer studies, the objective is to determine the maximum tolerated dose, defined as the highest dose with an acceptable dose-limiting toxicity rate. Finding this dose for drug-combination trials is complicated due to drug-drug interactions, and many trial designs have been proposed to address this issue. These designs rely on complicated statistical models that typically are not familiar to clinicians, and are rarely used in practice. The aim of this paper is to propose a Bayesian dose-finding design for drug combination trials based on standard logistic regression. Under the proposed design, we continuously update the posterior estimates of the model parameters to make the decisions of dose assignment and early stopping. Simulation studies show that the proposed design is competitive and outperforms some existing designs. We also extend our design to handle delayed toxicities.

Introduction

For oncologists, the objective of phase I dose-finding studies is to determine the maximum tolerated dose (MTD), defined as the highest dose with a relatively acceptable dose-limiting toxicity (DLT) [START_REF] Zhou | Choice of designs and doses for early phase trials[END_REF][START_REF] Chevret | Statistical Methods for dose-Finding Experiments[END_REF]. DLT is usually defined as a toxicity of grade 3 or higher according to the U.S. National Cancer Institute toxicity criteria [START_REF]Common Terminology Criteria for Adverse Events (CTCAE) and Common Toxicity Criteria (CTC)[END_REF]. In practice, patients included in phase I clinical cancer trials have already been heavily pre-treated and in many cases no alternative therapeutic options are available to them. For cytotoxic anti-cancer drugs, a dose-toxicity effect is assumed whereby the higher the dose, the greater the risk of DLT and the greater the efficacy.

Dose finding for drug combination trials is more difficult than that for conventional singleagent trials due to complicated drug-drug interactions. Moreover, when combining several agents, the order of the toxicity probabilities is not fully known. Should investigators wish to gradually increase the acceptable level of toxicity during the trial, the appropriate order in which the doses for the various drugs in the combination should be increased would be of great interest. For instance, when combining two cytotoxic agents, it remains difficult to decide how to escalate or de-escalate the dose combination, even when a partial ordering is known [START_REF] Conaway | Designs for single-or multiple-agent phase I trials[END_REF][START_REF] Wages | Continual Reassessment Method for Partial Ordering[END_REF].

Recently, many phase I dose-finding designs have been proposed for drug combination trials. Thall et al. proposed a Bayesian dose-finding method based on a six-parameter model [START_REF] Thall | Dose-finding with two agents in Phase I oncology trials[END_REF]. Wang and Ivanova developed a 3-parameter model-based method in which the parameters are estimated using Bayesian inference [START_REF] Wang | Two-dimensional dose finding in discrete dose space[END_REF]. Mandrekar et al. proposed an approach incorporating the toxicity and efficacy of each agent into the identification of an optimal dosing region for the combination by using a continuation ratio model to separate each agent's toxicity and efficacy curves [START_REF] Mandrekar | An adaptive phase I design for identifying a biologically optimal dose for dual agent drug combinations[END_REF][START_REF] Mandrekar | Model-based phase I designs incorporating toxicity and efficacy for single and dual agent drug combinations: methods and challenges[END_REF]. Yin and Yuan developed a Bayesian adaptive design based on latent 2x2 tables in which the combination's toxicity probabilities in the two-dimensional space are estimated using a Gumbel-type model [START_REF] Yin | A latent contingency table approach to dose finding for combinations of two agents[END_REF]. Yin and Yuan extended their method by changing to a copula-type model to simulate the effect of two or more drugs in combination [START_REF] Yin | Bayesian dose finding in oncology for drug combinations by copula regression[END_REF]. Bailey et al. introduced a second agent as a covariate in a logistic model [START_REF] Bailey | A Bayesian case study in oncology Phase I combination dose-finding using logistic regression with covariates[END_REF]. Wages et al. considered an approach based on the continual reassessment method and taking into account different orderings with partial order between combinations. In this case, the MTD is estimated for the order associated with the highest model-selection criterion [START_REF] Wages | Continual Reassessment Method for Partial Ordering[END_REF]. Most of these existing designs rely on complicated statistical models that typically are not familiar to clinicians, which hinder their acceptance and application in practice. In addition, the performance of these designs seems comparable and there is no consensus which design should be used [START_REF] Riviere | Competing designs for drug combination in phase I dosefinding clinical trials[END_REF]. As a result, many of the dose-finding clinical trials conducted to evaluate drug combinations still use the conventional "3+3" approach, which was developed for single agents and was shown to be inefficient in terms of dose identification [START_REF] Faries | Practical modifications of the continual reassessment method for phase I cancer clinical trials[END_REF][START_REF] Goodman | Some practical improvements in the continual reassessment method for phase I studies[END_REF][START_REF] Korn | A comparison of two phase I trial designs[END_REF][START_REF] O'quigley | Experimental designs for phase I and phase I/II dose-finding studies[END_REF][START_REF] Iasonos | A comprehensive comparison of the continual reassessment method to the standard 3+3 dose escalation scheme in Phase I dose-finding studies[END_REF].

In this paper, we propose an Bayesian dose-finding design for drug combination trials based on standard logistic regression. Under the Bayesian paradigm, data monitoring, early stopping and dose assignment occur continuously throughout the trial by updating the posterior estimates of the model parameters. Simulation studies show that in general the proposed design provide better performance than some existing designs. To facilitate the use of the proposed design by clinicians, R package will be developed for implementing the new design.

This manuscript is laid out as follows: in Section 2, we propose the simple statistical method for modeling the toxicity probabilities of the drug combination under evaluation. Moreover, we present the likelihood function and the prior specifications for the unknown parameters. In this section, we also describe our allocation and dose-finding method, as well as propose stopping rules.

We conduct extensive simulation studies to examine the operating characteristics of our design in Sections 3 and 4, and conclude with a discussion in Section 5.

Methods

Statistical method for combination evaluation (LOGISTIC)

Dose-combination model

Let there be a two-drug combination used in a phase I dose-finding clinical trial for which the dose-toxicity relationship is monotonic and increases with the dose levels. Let (j, k) denote the dose level of a combination in which j refers to Agent 1 (j = 1, . . . , J), and k refers to Agent 2 (k = 1, . . . , K). Y i is a Bernoulli random variable, denoting the toxicity that is equal to 1 if DLT occurs in patient i and 0 otherwise (i = 1, . . . , N ). We assume that n j,k patients are allocated at combination (j, k) and that a total of t j,k DLTs are observed for that combination. We define π j,k as the toxicity probabilities of combination (j, k), θ as the target probability of toxicity, and p 1 , ..., p J and q 1 , ..., q K as the respective prior toxicity probabilities of Agent 1 and Agent 2 taken alone. For simplicity purposes, we refer to the selected combination as the MTD in order to maintain the same designation as in single-agent trials. In this manuscript, we focus on finding one MTD at the end of the trial. We note that in some cases, it is of interest to find multiple MTDs that can be further tested in phase II trials, see Yuan and Yin [START_REF] Yuan | Sequential continual reassessment method for two-dimensional dose finding[END_REF] and Ivanova and Wang [START_REF] Ivanova | A non-parametric approach to the design and analysis of twodimensional dose-finding trials[END_REF] for related designs.

Let u j and v k denote the "effective" or standardized doses ascribed to the jth level of Agent 1 and k level of Agent 2, respectively. We model the drug combination-toxicity relationship using a 4-parameter logistic model, as follows:

logit(π j,k ) = β 0 + β 1 u j + β 2 v k + β 3 u j v k , ( 1
)
where β 0 , β 1 , β 2 and β 3 are unknown parameters that represent the toxicity effect of Agent 1 (β 1 ), that of Agent 2 (β 2 ) and that of the interaction between the two agents (β 3 ). These parameters are defined such that β 1 > 0 and β 2 > 0, ensuring that the toxicity probability is increasing with the increasing dose level of each agent alone, ∀k, β 1 + β 3 v k > 0 and ∀j, β 2 + β 3 u j > 0, ensuring that the toxicity probability is increasing with the increasing dose levels of both agents together, and intercept -∞ < β 0 < ∞. The standardized dose of two agents are defined as u j = log

( p j 1-p j ) and v k = log ( q k 1-q k )
, where p j and q k are the prior estimates of the toxicity probabilities of the jth dose level of Agent 1 and the kth dose level of Agent 2, respectively, when they are administered individually as a single agent. Before two agents can be combined, each of them typically have been thoroughly investigated individually. Therefore, there is often rich prior information on p j 's and q k 's and their values can be readily elicited from physicians. Using the prior information to define standardized dose has been widely used in dose-finding trial designs, and the most well known example perhaps is the skeleton of the continuous reassessment method (CRM) [START_REF] O'quigley | Continual reassessment method: a practical design for phase 1 clinical trials in cancer[END_REF] with a logistic model. Research has shown that this approach improves the estimation stability and trial performance [START_REF] Chevret | Statistical Methods for dose-Finding Experiments[END_REF][START_REF] Zohar | Using the continual reassessment method to estimate the minimum effective dose in phase II dose-finding studies: a case study[END_REF]. Our definition of u j and v k can be viewed as an extension of the skeleton of the CRM (with a logistic model) to drug-combination trials.

Likelihood and posterior inference

Under the proposed model, the likelihood is simply a product of the Bernoulli density, given by

L(β 0 , β 1 , β 2 , β 3 |data) ∝ J ∏ j=1 K ∏ k=1 π t j,k j,k (1 -π j,k ) n j,k -t j,k .
We assume the prior distributions of the model parameters are independent. For β 0 and β 3 , we assign a vague normal prior N (0, 10) centered at 0 to indicate that a priori we do not favor either positive or negative values for these parameters and let the observed data speak for themselves through posteriors. For β 1 and β 2 , we assume an informative prior distribution Exp(1) centered at 1, as those parameters should not be too far from 1. Then, the joint posterior distribution of parameters β 0 , β 1 , β 2 and β 3 is given by

f (β 0 , β 1 , β 2 , β 3 |data) ∝ L(β 0 , β 1 , β 2 , β 3 |data)f (β 0 )f (β 1 )f (β 2 )f (β 3 ). (2)
We sample this posterior distribution using Gibbs sampler, which sequentially draws each of the parameters from their full conditional distributions (see Appendix). These full conditional distributions do not have closed forms, and we use the Adaptive Rejection Metropolis Sampling (ARMS) method [START_REF] Gilks | Adaptive Rejection Metropolis Sampling within Gibbs Sampling[END_REF] to sample them. In order to impose the constrain that ∀k,

β 1 + β 3 v k > 0 and ∀j, β 2 + β 3 u j > 0,
3 ) ℓ=1,...,L denote the L posterior samples obtained from Gibbs sampler, the posterior toxicity probabilities can be estimated using Monte Carlo by:

πj,k = 1 L L ∑ ℓ=1 exp(β (ℓ) 0 + β (ℓ) 1 u j + β (ℓ) 2 v k + β (ℓ) 3 u j v k ) 1 + exp(β (ℓ) 0 + β (ℓ) 1 u j + β (ℓ) 2 v k + β (ℓ) 3 u j v k ) .

Dose finding algorithm and determination of the MTD

During the trial conduct, we use the dose-finding algorithm proposed by Yin and Yuan [START_REF] Yin | Bayesian dose finding in oncology for drug combinations by copula regression[END_REF][START_REF] Yin | A latent contingency table approach to dose finding for combinations of two agents[END_REF] to determine dose escalation and deescaltion, and propose a different criterion for MTD selection at the end of the trial. In our design, similar to Yin and Yuan [START_REF] Yin | Bayesian dose finding in oncology for drug combinations by copula regression[END_REF], we restrict dose escalation and deescalation one level at a time (i.e, we do not allow dose escalate or de-escalate along the diagonal) based on the practical consideration that physicians are conservative and typically do not allow two agent to escalate at the same time for patient safety. Nevertheless, we note that Sweeting and Mander [START_REF] Sweeting | Escalation strategies for combination therapy Phase I trials[END_REF] showed that diagonal escalation strategy may be more efficient in reaching the target toxicity level quicker with fewer patients treated at sub-optimal doses and have a higher percentage of correct selection at the end of the trial. We take this strategy in our start-up phase described later.

Let c e be the probability threshold for dose escalation and c d the probability threshold for dose de-escalation. We require c e +c d > 1 to avoid that the decisions of dose escalation and deescalation occur at the same time. Our dose-finding algorithm can be described as follows:

• If the current combination is (j, k) and P (π j,k < θ|data) > c e , we escalate the combination dose level to an adjacent combination dose level {(j + 1, k), (j, k + 1), (j + 1, k -1), (j -1, k + 1)} that has a toxicity probability that is higher than the current value and closest to θ. If the current combination is the highest of the combination space, (J, K), we retain the same combination dose level for the next cohort.

• If the current dose combination is (j, k) and P (π j,k > θ|data) > c d , we de-escalate the combination dose level to an adjacent combination dose level {(j -1, k), (j, k -1), (j + 1, k -1), (j -1, k + 1)} that has a toxicity probability that is lower than the current value and closest to θ. If the current combination is the lowest of the combination space, (1, 1), we retain the same combination dose level for the next cohort.

• If the current combination is (j, k) and P (π j,k < θ|data) ≤ c e and P (π j,k > θ|data) ≤ c d , we treat the next cohort of patients at the current combination dose level.

Once we reach the maximum sample size, we select the MTD as the combination associated with the highest posterior probability, P (π j,k ∈ [θ -δ; θ + δ]), and have been used to treat at least one cohort of patients. If, for example, the target toxicity θ is 0.3, and the length around the targeted interval is δ = 0.1, which gives a targeted interval defined as [0.2; 0.4], then, as illustrated in Figure 5, at the end of the trial we can obtain the posterior densities of the toxicity probability for each combination. Shown in light gray in the figure is the area under the curve (AUC) for a toxicity probability lower than 0.2, which is equal to the probability of under-dosing. Shown in medium gray is the AUC for a toxicity probability between 0.2 and 0.4, which is the probability of being in the targeted interval. Shown in dark gray is the AUC for a toxicity probability greater than 0.4, which is equal to the probability of overdosing. For each combination, given the probabilities of being in the targeted interval, the combination already administered to at least one cohort at the end of the trial and corresponding to the highest probability is selected as the MTD. For example, dose combination (2, 3), as illustrated in Figure 2.

The probability thresholds c e and c d are critical for the performance of the design as they control the dose escalation and deescalation. The values of c e and c d should be carefully calibrated through simulation to ensure good operating characteristics of the design. In practice, this can be done as follows: first define a set of representative dose-toxicity scenarios that may be encountered in the trial, and then conduct simulation under different values of c e and c d to evaluate the performance of the design. This is a trial-and-error process and may involve repeatedly tuning the values of c e and c d based on the simulation results. The goal is to find the values of c e and c d that yield good overall performance across different scenarios (e.g., the percentage of correct selection of the MTD, the number of patients exposed to over-toxic combinations or under-toxic combinations). Such a calibration based approach has been widely used in clinical trial designs [START_REF] Thall | Practical Bayesian guidelines for phase IIB clinical trials[END_REF][START_REF] Thall | Dose-finding with two agents in Phase I oncology trials[END_REF][START_REF] Yin | Bayesian dose finding in oncology for drug combinations by copula regression[END_REF][START_REF] Yin | A latent contingency table approach to dose finding for combinations of two agents[END_REF].

Due to the limited availability of information at the beginning of the trial, the posterior estimates of the toxicity probabilities may not be reliable. Therefore, as suggested by other authors [START_REF] Yin | A latent contingency table approach to dose finding for combinations of two agents[END_REF][START_REF] Korn | A comparison of two phase I trial designs[END_REF][START_REF] Storer | Design and analysis of phase I clinical trials[END_REF][START_REF] Ivanova | Improved up-and-down designs for phase I trials[END_REF][START_REF] Ivanova | A non-parametric approach to the design and analysis of twodimensional dose-finding trials[END_REF], we implement an algorithm-based start-up phase in order to gather enough information to estimate the π j,k . Our start-up phase shares the spirit of accelerated titration design [START_REF] Simon | Accelerated titration designs for phase I clinical trials in oncology[END_REF] and can be described as follows: treat the first cohort of patients at the lowest dose combination (1, 1); if no toxicity is observed, escalate the dose along the diagonal and treat the second cohort of patients at (2, 2); we continue this dose escalation along the diagonal until the observation of a toxicity or at least one agent is at its maximum. If one of the agents is already at its maximum dose level and still no toxicity is observed, we increase the dose of the other agent until both agents reach their maximum dose level. At any time of the start-phase, once the first toxicity event is observed, the start-up phase is completed and the design switches to the model-based dose-finding algorithm as described above.

Stopping rules

Our aim is to propose a design that, should the investigators wish, will stop the trial when all combinations are estimated to be unacceptable in terms of toxicity. If the current combination is (1, 1), at least two cohorts have been included and P (π j,k < θ) ≥ 0.975, then the trial is stopped.

Some existing designs for combination studies BCOPULA and BGUMBEL methods

Yin and Yuan proposed two Bayesian methods that use copula regression for combinations. The first method [START_REF] Yin | Bayesian dose finding in oncology for drug combinations by copula regression[END_REF] uses a Clayton-copula regression in order to express the joint probability of combination (j, k) with the marginal true probabilities of toxicity (p α j , q β k ) :

π j,k = 1 - ( (1 -p α j ) -γ + (1 -q β k ) -γ -1 ) -1 γ .
The second design [START_REF] Yin | A latent contingency table approach to dose finding for combinations of two agents[END_REF] uses a Gumbel model:

π j,k = 1 -(1 -p α j )(1 -q β k ) [ 1 + p α j q β k e γ -1 e γ + 1
] ,

where γ, α, β > 0 are unknown parameters. The parameter γ characterizes the drug interactive effect, and α and β the uncertainty of the initial guesses. The combination allocation algorithm is the same as that presented in Section 2.1. The final MTD is the combination with a toxicity probability closest to the target among the combinations already administered in the trial.

I2D method

Wang and Ivanova proposed the I2D method [START_REF] Wang | Two-dimensional dose finding in discrete dose space[END_REF], which is a two-dimensional model-based design, defined as follows [START_REF] Wang | Two-dimensional dose finding in discrete dose space[END_REF]:

π j,k = 1 -(1 -p j ) α (1 -q k ) β exp (-γ log(1 -p j ) log(1 -q k )) ,
where α > 0, β > 0, the interaction term γ > 0, and p j , q k are the working models for Agent 1 and Agent 2, respectively.

The combination that will be allocated to the next cohort is the combination closest to the target that belongs to {(j + 1, k), (j, k + 1), (j -1, k + 1), (j -1, k), (j, k -1), (j + 1, k -1), (j, k)}. The final MTD is defined as the combination with a toxicity probability closest to the target among the combinations already administered in the trial.

LOGODDS

Gasparini, Bailey and Neuenschwander [START_REF] Gasparini | Bayesian dose finding in oncology for drug combinations by copula regression[END_REF] have proposed an alternative model, π j,k , to the copula proposed by Yin and Yuan [START_REF] Yin | Bayesian dose finding in oncology for drug combinations by copula regression[END_REF] containing an interaction parameter γ which explicitly quantifies departure from no interaction on the log-odds scale:

γ = logodds(π j,k ) -logodds(π ⊥ j,k
) where π ⊥ j,k is the no-interaction model defined as π ⊥ j,k = p j + q k -p j q k . By backsolving this equation, an explicit expression for the probability of toxicity is obtained. A normal prior centered in 0 and with variance of 100 was chosen for the interaction parameter. The rest of the dose allocation process, estimation and MTD determination was the same as our proposed design in order to compare our method involving a simple interaction logistic model with other logistic models.

TMML

Thall et al. [START_REF] Thall | Dose-finding with two agents in Phase I oncology trials[END_REF] proposed a six-parameter model defined as follows:

π j,k = α 1 x β 1 1 + α 2 x β 2 2 + α 3 (x β 1 1 x β 2 2 ) β 3 1 + α 1 x β 1 1 + α 2 x β 2 2 + α 3 (x β 1 1 x β 2
2 ) β 3 where α 1 , β 1 , α 2 , β 2 , α 3 , β 3 are unknown positive parameters. Following Thall [START_REF] Thall | Dose-finding with two agents in Phase I oncology trials[END_REF], a prior gamma distribution with parameters (0.3, 0.3) was chosen for α 1 , α 2 and α 3 , and a gamma distribution with parameters (0.003, 0.03) was chosen for β 1 , β 2 and β 3 . The rest of the dose allocation process, estimation and MTD determination was the same as our proposed design in order to compare our method involving a simple interaction logistic model with other logistic models.

One-dimensional CRM

In practice, one-dimensional CRM sometimes is used to conduct dose-combination trials [START_REF] Riviere | Competing designs for drug combination in phase I dosefinding clinical trials[END_REF]. Under this method, we first preselect a subset of combinations, for which the toxicity probability order is known, and then apply the standard CRM to find the MTD. The drawback of such an approach is that we only investigate a subset of the whole two-dimensional dose space and may miss the true MTD. Following [START_REF] Riviere | Competing designs for drug combination in phase I dosefinding clinical trials[END_REF], we chose the subset of combinations as the combinations located at anti-diagonal path:

D 1,1 -→ D 1,2 -→ D 2,2 -→ D 3,2 -→ D 4,2 -→ D 5,2 -→ D 5,3 .

Simulations

We simulated 2000 independent replications of phase I trials that evaluate two agents in drug combinations, with five dose levels for Agent 1 and three for Agent 2, giving 15 possible combinations. We simulated 14 scenarios to represent the possible true underlying combination toxicities (Table 1). We studied several locations of the MTD in the combination space, as well as a number of correct combinations (or MTDs) that could shift from 0 to 3. We fixed the toxicity target at 0.3 and used an overall sample size of 60. In this paper, we compare the performance of our method with those of five other designs, including the BCOPULA [START_REF] Yin | Bayesian dose finding in oncology for drug combinations by copula regression[END_REF] and BGUMBEL [START_REF] Yin | A latent contingency table approach to dose finding for combinations of two agents[END_REF] methods proposed by Yin and Yuan, and the I2D method proposed by Wang and Ivanova [START_REF] Wang | Two-dimensional dose finding in discrete dose space[END_REF], and with two other logistic models, LOGODDS proposed by Gasparini et al. [START_REF] Gasparini | Bayesian dose finding in oncology for drug combinations by copula regression[END_REF] and TMML proposed by Thall [START_REF] Thall | Dose-finding with two agents in Phase I oncology trials[END_REF].We also compared the performance of this multidimensional designs with a one-dimensional CRM, CRM anti-diag, as described previously. The design parameters of all the designs (e.g., working model) have been calibrated via simulation before used for the comparison.

To ensure comparability, we chose a cohort size of 3 for all methods and did not use stopping rules when comparing all the methods. The start-up phase was implemented in all simulations. Due to ethical concerns, we started each trial at the lowest dose combination (1, 1). We selected these features in order to employ typical trial set-ups in our simulation study. For Agent 1, we specified the marginal initial guesses of toxicities, p j , as (0.12, 0.2, 0.3, 0.4, 0.5), and for Agent 2, we specified q k as (0.2, 0.3, 0.4). The same working models were used for BCOPULA, BGUMBEL, LOGODDS and I2D. We set the length around the targeted interval, δ, at 0.1. Based on a sensitivity analysis (data not shown), we set the probability thresholds c e and c d at 0.85 and 0.45. We recorded 5000 posterior samples of the model parameters after 2000 burn-in iterations to make inference. After comparing the proposed design to the other designs, we also investigated the performance of our design with the stopping rule that we introduced above. For the one-dimensional CRM, we used the "dfcrm" R package with restrictions to avoid skipping doses. The working model was generated using the "getprior" function with an indifference interval δ = 0.05, a initial guessed MTD at dose level 4 for a trial with 7 doses, and the same settings as for multidimensional designs were used.

Results

For each scenario, we present the correct MTD selection rate, or percentages of correct selection (PCS), in Table 2. In general, the proposed method for combination evaluation (LOGISTIC) performed better than the other existing model-based designs. Indeed, except in scenarios 4 and 5, in which the PCS were already very high (86.7% and 80.4% respectively), those in other scenarios were either the best or not less than 2% of those of the other designs. Moreover, the PCS were greater than 55% in all scenarios and higher than 60% in 12 scenarios out of 14.

Scenarios 1, 3 and 9 included three possible MTDs that were located on one diagonal of the combination space; combinations (2, 3), (3, 2), (4, 1) in scenario 1, combinations (3, 3), (4, 2), (5, 1) in scenario 3 and combinations (1, 3), (2, 2), [START_REF]Common Terminology Criteria for Adverse Events (CTCAE) and Common Toxicity Criteria (CTC)[END_REF][START_REF] Zhou | Choice of designs and doses for early phase trials[END_REF] in scenario 9 (Table 1). For these scenarios, all model-based designs gave high PCS, whereas LOGISTIC seemed to perform better than the other methods. In addition, for the LOGISTIC, the distribution of allocated dose levels around the MTDs (Table 3) seemed to be spread unevenly between possible MTDs.

In scenarios 2 and 10, in which two possible MTDs were located on the same diagonal at the lower end or at the higher end of the combination space, high PCS were observed for the LOGISTIC (> 75%). When the correct combination was in the lower (1, 1) extremity of the combination space, as in scenario 4, the TMML model and then I2D method performed best (92.1% and 89.7% respectively), but the PCS for the LOGISTIC was close to this value (86.7%). By contrast, when the correct combination was in the higher (5, 3) extremity of the combination space, as in scenario 5, the LOGODDS model gave the highest PCS with 90.0%.

In scenarios 6, 7, 11, 12, 13 and 14, when the true MTDs were not located on the same diagonal but randomly located in the combination space, the proposed method performed best with the PCS that was consistently higher than 56%. In contrast, the PCS of other methods can be as low as 30%. For example, under scenario 6, the PCS of BCOPULA, BGUMBEL, I2D, and LOGODDS were all lower than 47%; and under scenario 12, the PCS of TMML was only 47.9%.

Finally, when the true MTD was unique and located in the middle of the combination space, as in scenario 8, the PCS decreased. The TMML gave the best PCS (69.7%), nevertheless the LOGISTIC was the second most efficient design and gave good results with PCS of 56.9% for that scenario (Table 2). The PCS of the LOGODDS was merely 12.1%.

Under our simulation setting, the one-dimensional CRM (based on an anti-diagonal subset of combinations) yielded higher PCS than the multidimensional designs. However, this result should not be generalized. The reason that the one-dimensional CRM performed very well is simply because the pre-selected subset of combinations happened to include the true MTD. As the CRM focused on only a subset of doses, it had more patients per dose to find the target doses than the multidimensional designs, given the same total number of patients. It is easy to see that if the true MTD(s) is (are) not contained in the pre-selected subset of ordered combinations, then the one-dimensional CRM will perform badly and can never select the true MTD.

The mean number of DLTs observed throughout the trial in each scenario when using the LO-GISTIC was similar to that observed when using the BCOPULA, BGUMBEL, and I2D methods (Table 2). Over all the scenarios, the mean number of DLTs observed during the trial was 15.4, which is lower than 18, the expected number of DLTs for 60 patients with a 0.30 toxicity probability.

According to Table 2, the LOGISTIC tended to treat as many patients at the MTD combination as the other methods. The lowest percentage of patients treated at the MTD combination was 24.3% in scenario 8, and the highest was 78.0% in scenario 4. In general, the true MTD combination was assigned to at least 1/3 of the patients, and the mean overall percentage of patients treated at the MTD combination was 42.6%.

After employing the stopping rule for all scenarios when using the LOGISTIC, we added a scenario in which all the dose combinations were unacceptably toxic (Table 4). The addition of the stopping rule for unacceptable toxicity resulted in PCS that were similar to those presented above (Table 5), except in scenario 4 in which the first dose combination, (1, 1), was the MTD. For this scenario, the PCS decreased by 16.9%, although it was nevertheless high (69.8%). In practice, the trial would not be completely stopped under these circumstances, but lower dose combinations would be added to the trial. In the additional scenario 15, where all dose combinations were too toxic, the trial was stopped in 83.7% of the cases. Therefore, the stopping rule seemed efficient.

Sensitivity analysis

We conducted a sensitivity analysis in order to study the performance of our design using different prior distributions and parameters values. We used normal distributions centered at 0 for β 0 , β 3 with high or low variance (from 10 to 50), and exponential or gamma distributions centered at 1 for β 1 , β 2 with high or low variance (from 1 to 10). According to Table 5, we can see that the PCS for all scenarios were very similar under these different prior distributions.

Time-to-event outcome

In practice, a longer follow-up time may be required to assess the toxicity outcome. Therefore, the toxicity outcomes of some patients already treated in the trial may be unobserved (or censored) when a new patient is enrolled in the trial and is ready for dose assignment. Waiting to assess toxicity outcomes for each cohort before including a new one in the trial can greatly increase the duration of the trial. To overcome this issue, we extended our method by modeling toxicity as a time-to-event outcome. Following the TITE-CRM [START_REF] Cheung | Sequential designs for phase I clinical trials with late-onset toxicities[END_REF], we considered a weighted dose-toxicity relationship ŵF (D, β 0 , β 1 , β 2 , β 3 ) where ŵ is monotone and increasing with patient follow-up time such that 0 ≤ ŵ ≤ 1, and the toxicity probability model F (D, β 0 , β 1 , β 2 , β 3 ) = π j,k is the same as that in [START_REF] Zhou | Choice of designs and doses for early phase trials[END_REF].

Let T be a (maximum) fixed time window during which patients are followed, and y i,N , C i,N , and ŵi,N , respectively, be the indicator of a DLT for patient i prior to the entry of the (N + 1) th patient, the follow-up time of patient i prior to the entry of the (N + 1) th patient, and the weight assigned to patient i prior to the entry of the (N + 1) th patient. Let X i denote the time to toxicity of the i th patient, and (j i , k i ) the combination received by patient i. Before combination assignment, the likelihood is defined as

L(β 0 , β 1 , β 2 , β 3 |data) = N ∏ i=1 ( ŵi,N π j i ,k i ) y i,N (1 -ŵi,N π j i ,k i ) . 1-y i,N
Following Cheung and Thall [START_REF] Cheung | Monitoring the rates of composite events with censored data in phase II clinical trials[END_REF], we chose the weights ŵi,N as

ŵi,N = # {m/X m ≤ C i,N , C m,N ≥ T } + ŵ0 i,N # {m/X m ≤ T, C m,N ≥ T } + 1 ,
where m refers to patient, # {m/X m ≤ T, C m,N ≥ T } is the number of completely followed patients on who a toxicity was observed, and ŵ0 i,N is the linear weight for patient i defined as ŵ0 i,N = C i,N /T , that is the proportion of time patient i was followed compared to the full follow-up time.

We simulated the time-to-toxicity outcomes using an exponential distribution such that the toxicity probabilities at the end of follow-up matched those given in Table 1. Using the same working model, we chose a full follow-up time of 3 weeks, and simulated patient accrual using a Poisson process with parameter 1, meaning that an average of one patient arrived every week.

Table 6 shows the results of the extended LOGISTIC for all 15 scenarios. We observe that the performance of the design decreases only slightly by 2%, and the PCS for all scenarios are still very high. This demonstrates that the LOGISTIC can be used in trials when toxicity cannot be assessed quickly.

Discussion

We have proposed a statistical method for clinical trial designs that evaluate various dose combinations for two agents. This method seems to be capable of identifying the MTD with a high percentage in all scenarios. Indeed, this method works for a wide variety of drug combinations. Our aim was not to find the correct order of toxicity within the drug combination space but to identify the right dose combination to be evaluated further in terms of efficacy. We compared our method with five other model-based designs for combination drug trials. All the designs seem to be efficient when the MTDs are located on the same diagonal in the combination space. One benefit of our method compared with the other proposed designs is that it is also efficient when the MTDs are not necessarily located on the same diagonal.

We have designed our method to enable the selection of only one MTD at the end of the trial. There are various reasons why investigators may wish to select one MTD per level of Agent 2. For this purpose, our method can be modified by selecting, for each level k of Agent 2, the combination (j k , k) such that j k = argmax 1≤j≤J P (π j,k ∈ [θ -δ; θ + δ]). In the case where it cannot be assumed that one MTD exists per level of Agent 2, a minimum probability of being in the targeted toxicity interval or a maximum probability of overdosing can be added in the determination of the MTD so as not to recommend a combination in the case of an MTD not existing for this dose level of Agent 2. Our design considers only toxicity responses, which can be improved upon by taking into account efficacy when determining the combination to be used in subsequent phases.

When combining several agents, designs developed for single-agent dose-finding trials cannot be applied to combination studies. For instance, the standard algorithm-based method for phase I dose-finding clinical trials in oncology is the so-called "3+3" design, which is referred to as "memory-less" since allocation to the next dose level for an incoming group of 3 patients depends only upon what has happened to the total of 3 to 6 patients previously treated at the current dose level [START_REF] Faries | Practical modifications of the continual reassessment method for phase I cancer clinical trials[END_REF][START_REF] Goodman | Some practical improvements in the continual reassessment method for phase I studies[END_REF][START_REF] Korn | A comparison of two phase I trial designs[END_REF][START_REF] O'quigley | Experimental designs for phase I and phase I/II dose-finding studies[END_REF][START_REF] Iasonos | A comprehensive comparison of the continual reassessment method to the standard 3+3 dose escalation scheme in Phase I dose-finding studies[END_REF]. This method was not designed for cases in which the full toxicity ordering is unknown. Indeed, the full ordering between toxicity probabilities of each combination is not fully known. That is why some investigators fix the dose level of one agent when they set up a combination study. They can thereby bring the problem back to one dimension and use singleagent designs. Another approach studied in this paper is to pre-select a subset of combinations whose toxicity ordering is known, and then apply a single-agent design to this reduced number of combinations. This approach performs well if the target dose combinations happen to be included in the subset. However, in practice, due to potential interaction between combined drugs, it can never be guaranteed that the target doses are always pre-selected into the subset of combinations to be inverstigated. If target doses are not in the subset, the trial would completely fail to find the target doses.

The advantage of our method is that it uses a logistic model that can be understood by a broad panel of readers and the performance seems to outperform existing designs. We also added stopping rules to the design in order to accommodate practical issues. The method can be efficient with different cohort sizes. Nevertheless, estimating parameter values and the different probabilities continues to require some advanced computational skills. To overcome this challenge, we have developed executable files that can be used to (1) determine the next combination and the MTD and estimate toxicity probabilities from data in actual trials, and (2) perform simulations before setting up a combination dose-finding trial. These files are freely available upon request. (

(5,1)

(5,2)

(5,3)

Probability of targeted interval Combination 0.0 0.2 0.4 0.6 0.8 1.0

4 Figure 2 :

 42 Figure 2: Probability of achieving the targeted interval, P (π j,k ∈ [θ -δ; θ + δ]) at the end of the trial for each combination. The dashed vertical line corresponds to the highest probability, which determines the combination selected as the MTD ((2, 3) in this example).

  at each iteration of Gibbs sampling, if sampled β 1 , β 2 and β 3 fail to satisfy the constraint, we re-sample β 1 , β 2 and β 3 . Let (β

Table 1 :

 1 Toxicity scenarios for the two-agent combinations.

	Agent 1

The true MTD combinations are shown in boldface.

Table 6 :

 6 Percentage of correct selection (PCS) and percentage of stopping for toxicity without MTD recommendation when toxicity is

	considered as a time-to-event outcome.	scenario sc1 sc2 sc3 sc4 sc5 sc6 sc7 sc8 sc9 sc10 sc11 sc12 sc13 sc14 sc15	TITE-LOGISTIC	PCS 74.0 77.5 74.1 68.0 78.2 63.3 66.5 57.6 64.7 74.6 76.2 56.5 60.4 62.7

Acknowledgements

This work was partially funded by grants from the ANRT (Association Nationale de la Recherche et de la Technologie), "Laboratoires Servier" (CIFRE number 2011/0900), and Denis Diderot University-Paris 7 ("Aide à la mobilité internationale 2012").

Appendix Gibbs sampling

Gibbs sampling was used to generate a sequence of samples from the joint probability distribution of β = (β 0 , β 1 , β 2 , β 3 ). The Gibbs sampler relies on the availability of all complete conditional distributions. We start from an arbitrary point β (0) = (β

3 ) in the distribution, we sampled in turn in each of the full conditional distribution in updating them as the process goes on.

(1) Initialize β (0) = (β

3 )

For m from 1 to M + N

(2) Sample

We chose M = 2000 number of "burn-in" iterations to be discarded before convergence.

) converges in distribution to the posterior joint distribution f (β 0 , β 1 , β 2 , β 3 |data). The following N = 5000 values retained are then considered as a sample from this distribution.

Adaptive Rejection Metropolis Sampling within Gibbs sampling

Adaptive Rejection Sampling (ARS) can not be used to sample from non log-concave distributions. When this is the case, Gilks, Best and Tan [START_REF] Gilks | Adaptive Rejection Metropolis Sampling within Gibbs Sampling[END_REF] propose to replace the rejection sampling in favor of the Hastings-Metropolis algorithm to update one parameter at a time. But to avoid high probabilities of rejection, they adapted the proposal density to the shape of the full conditional density using ARS. They added to ARS a single Hastings-Metropolis step thus creating ARMS within Gibbs sampling. The detailed method can be found in Gilks [START_REF] Gilks | Adaptive Rejection Metropolis Sampling within Gibbs Sampling[END_REF]. 0.4 11.9 12.0 0.9 0.0 0.0 0.0 0.7 3.5 0.9

The percentages for the true MTD combinations are shown in boldface.


R Yin

N Lin

Marie-Karelle Wages

Frédéric Riviere

Sarah Dubois

Zohar

Marie-Karelle Riviere
email: marie-karelle.riviere@crc.jussieu.fr

Response to comments on 'Competing designs for drug

N. Wages

As pointed out by G. Yin and R. Lin and N. Wages, designs for drug-combination trials need to be rigorously evaluated. In our paper we performed an extensive sensitivity analysis (with some of the results given in the form of additional material) and for each design we tried to find the best parameter settings for model calibration associated with a high percentage of correct selection (PCS) on average across all scenarios.

Replying to G. Yin and R. Lin concerning the BCOPULA design, we computed this design using the C code directly provided by Y. Yuan [START_REF] Yuan | Robust EM Continual Reassessment Method in Oncology Dose Finding[END_REF], (the C code used in the original manuscript describing the method). Since we only removed the stopping rule there could not be any programming issues. When using the same setting as G. Yin and R. Lin (the authors of the comments), we obtained similar results (see Table below).

In their commentary, Yin and Lin did not use the same designs as in our paper, as (1) the working model was different in that the target toxicity was at the highest dose level; [START_REF] Wages | Continual reassessment method for partial ordering[END_REF] the spacing between prior toxicity probabilities was different from ours; and (3) the cut-offs for BCOPULA and BGUMBEL were not the same as those detailed in the manuscript. However, when we ran simulations using exactly the same settings as stated in the correspondence of Yin [START_REF] Yuan | Robust EM Continual Reassessment Method in Oncology Dose Finding[END_REF][START_REF] Wages | Dose finding design for multi-drug combinations[END_REF] 53.9 69.5 54.9 87.2 88.9 41.2 52.3 29.8 12.8 30.9

We ran simulations with a cohort size of 3 throughout the trial, and then a cohort size of 1 for the start-up phase and a cohort size of 3 for the model-based phase, as it seems this is what Yin and Lin chose in the commentary. We can observe that our results are similar to those of Yin and Lin.

Marie-Karelle Riviere Yin and Lin argue that our working model for the combination of agents is wrong. Firstly, some clinical trials are performed where only one agent is studied in phase I with the other agent being a new molecule for which no prior information is available. In this case, there is no reason to limit the settings of the working model. Secondly, the spacing between the prior toxicity probabilities in our choice of working model is not the same as that of Yin and Lin, and this choice also impacts on performance. Finally, our choice was mainly based on the sensitivity analyses we performed to assess the influence of several working models and their impact on performance. Moving the targeted toxicity from the lowest dose to the highest dose increases PCS in some scenarios and decreases them in others. Choosing the working model with the target ranked in the middle was a good compromise. Of course, this choice is subject to debate (as with any other choice). A best working model could be chosen for each scenario but there was no specific working model that gave the best performance in all cases. In addition, the choice of the working model should also be made according to the information available about each agent before beginning a clinical trial. In our paper, the choices of the working models for the studied designs were performed for the first eight scenarios, before scenarios 9 and 10 were added during the manuscript's review process. It showed that, depending on the retained scenarios and the clinical setting, the working model should be recalibrated. This was discussed in detail in our paper.

Replying to N. Wages, as stated in our paper, we used the POCRM design starting at the lowest dose level, based on current practice in Phase I oncology clinical trial. We chose to study the Bayesian version of POCRM [START_REF] Wages | Continual reassessment method for partial ordering[END_REF] by analogy to the CRM that is used in most published clinical trials under the Bayesian framework. In [START_REF] Wages | Continual reassessment method for partial ordering[END_REF], a start-up phase was proposed which was different from the one mentioned by Wages in this commentary and the dose allocation did not begin at the lowest dose level. The authors proposed in [START_REF] Wages | Continual reassessment method for partial ordering[END_REF] a start-up phase where at each step the ordering selected was sampled randomly with the weights of posterior probabilities. This start-up phase was not used by the authors in [START_REF] Wages | Continual reassessment method for partial ordering[END_REF] in the simulation study and we chose to proceed as they did. Indeed, in the MLE POCRM paper [START_REF] Wages | Dose finding design for multi-drug combinations[END_REF] published soon after, a start-up phase, necessary in the frequentist context, was proposed.

As for the Likelihood CRM, the exploration begins at the lowest dose level, escalates by diagonal zones, and ends at the first observation of a DLT. We agree that this start-up phase could be used in the Bayesian POCRM version to begin safely at the lowest dose level. This point was discussed in our paper, we agree that a safe start-up phase is required for the Bayesian POCRM [START_REF] Wages | Continual reassessment method for partial ordering[END_REF], we also have proposed one.

However, our choice of starting dose for comparison purpose did not impact the PCS of the Bayesian The results presented in our paper reflect the PCS of the Bayesian POCRM [START_REF] Wages | Continual reassessment method for partial ordering[END_REF] independently of the starting dose chosen. Regarding the comment concerning the R package [START_REF] Wages | pocrm: an R package for phase I trials of combinations of agents[END_REF] mentioned by Wages, this package does not allow implementing the Bayesian POCRM. Moreover, (1) in particular cases ther e can be large difference in performance between the two versions of the POCRM; and (2) the package was not flexible enough for the choice retained in our comparative analysis, for example the size of the start-up phase, as well as other parameters, cannot be changed.

Concerning the choice and the number of orderings we have performed our own sensitivity analysis, as the guidelines proposed in [START_REF] Wages | Specifications of a continual reassessment method design for phase I trials of combined drugs[END_REF] were not available on the public domain until late in the review process of our manuscript. However, now that these guidelines are available, we recommend using them. Our Then in the "22" and "200 orderings" used in the additional material, the previous 7 orderings were included and the remaining orderings to select were sampled in the set of all possible orderings satisfying the partial order (composed of 6006 possible orderings). This sensitivity analysis on the number of orderings performed on the first eight scenarios led to our choice of three orderings (sc9 and sc10 were added lately during the review as mentioned above).

We believe that, in general, it is always possible to find some scenarios that tumble the performance of a method. Scenarios 9 and 10, that are peculiar (as they do not respect the "zoning" assumption often supposed in combination designs as well as in the start-up phase proposed by Wages [START_REF] Wages | Dose finding design for multi-drug combinations[END_REF]), were added to highlight the benefit of multi-dimensional designs vs one-dimensional designs. Therefore, we did not observe the lack of information due to the 3 orderings compare to 6 orderings. The use of 6 orderings seems to be an efficient setting for the POCRM as it increases the performance on the last two scenarios. Nevertheless, as we mentioned above, we believe that it is always possible to find some scenarios that tumble the performance of a method. For example, if we choose the following scenari o: 0.30 0.45 0.50 0.60 0.75 0.05 0.10 0.15 0.45 0.50 0.01 0.02 0.05 0.08 0.15 then, the PCS for the POCRM given the six orderings proposed by Wages is only equal to 3.4%.

Therefore, in general, the performances of all model-based methods seem to remain comparable.

Yin and Lin, and Wages have raised an important question concerning model tuning or calibration.

Calibrating such complex models requires rigorous simulation studies. As an example, the model -based continual reassessment method (CRM) using a 1-parameter dose-toxicity model has been the subject of many papers concerning its calibration. Numerous authors have tried to highlight how to calibrate the CRM before and during a dose finding clinical trial, including for ex ample Cheung and Chappell [START_REF] Cheung | A simple technique to evaluate model sensitivity in the continual reassessment method[END_REF],

Lee and Cheung [START_REF] Lee | Model calibration in the continual reassessment method[END_REF][START_REF] Lee | Calibration of prior variance in the Bayesian continual reassessment method[END_REF], Zhang et al. [START_REF] Zhang | Adaptive prior variance calibration in the Bayesian continual reassessment method[END_REF], Yuan and Yin [START_REF] Yin | Bayesian dose finding in oncology for drug combinations by copula regression[END_REF], Daimon et al. [START_REF] Daimon | Posterior maximization and averaging for Bayesian working model choice in the continual reassessment method[END_REF] among others. For the CRM the choice of the working model and the prior (when the Bayesian framework is used) considering only a one-parameter model has been the subject of much debate and there is still no consensus. In combination model-based designs, the models are more complex and there are more parameters to calibrate. The parameter settings used in a methodological manuscript of a proposed mod el are not necessarily the ones that are the most suitable for a specific indication and disease. Therefore, before the beginning of each new dose-finding clinical trial, an extensive simulation study should be performed in order to calibrate the design.

In conclusion, based on the assumptions we have made, our numerical results are valid when choosing the same setting values. The choice of the working model and calibration is important but much more complex in the case of combination studies. Given a set of fixed scenarios, we can always tune the design to make it better, but this may not work in all situations and all scenarios.

Indeed, for model-based combination methods, this aspect still needs to be studied as it was done for the CRM. In our manuscript, we did not favor any design over the others; this should be left for discussion between the members of the team involved in the planning of the clinical trial. The choice of the design should be made according to the clinical setting, the indication, the disease, the knowledge about each agent prior to trial onset and an extensive sensitivity analysis.

SIM- 12 -

 12 0373 and sim-14-0573: Response to 'Comments on Competing designs for drug combination in phase I dose-finding clinical trials' by G. Yin, R. Lin and

  sensitivity analysis included several settings: choice of orderings, number of orderings and working model (given in additional material). In the following figure, orderings (a), (b) and (c) correspond to the three orderings, and orderings (a) to (g) to the seven orderings used in the sensitivity analysis on the number of orderings given in the additional material.

  and Lin, we obtained the following

	results:
			sc1	sc2	sc3	sc4	sc5	sc6	sc7	sc8	sc9	sc10
		Yin and Lin	57.9 68.5 58.6 84.4 91.6 40.6 49.6 30.4 19.6 29.5
	Copula	Riviere et al. cohort size (3,3)	60.4 69.9 59.1 85.6 91.5 43.3 49.4 33.1 16.3 31.9
		Riviere et al. cohort size (1,3)	53.2 67.8 55.5 85.5 88.9 40.4 51.3 25.7 12.8 31.1
		Yin and Lin	55.3 70.3 58.6 85.2 91.4 42.2 50.2 27.4 19.5 29.4
	Gumbel	Riviere et al. cohort size (3,3)	58.7 71.5 58.1 87.2 90.1 45.3 51.9 33.9 22.8 28.7
		Riviere et al.
		cohort size

  1,2,* , Frédéric Dubois 2 and Sarah Zohar 1

	1 INSERM, U1138, Equipe 22 Centre de Recherche des Cordeliers, Université Paris 5, Université Paris 6, France
	2 IRIS (Institut de Recherches Internationales Servier), Suresnes, France