Randomized dose-escalation designs for drug combination cancer trials with immunotherapy

ABSTRACT

This work considers Phase I cancer dual-agent dose-escalation clinical trials in which one of the compounds is an immunotherapy. The distinguishing feature of trials considered is that the dose of one agent, referred to as a standard of care, is fixed and another agent is dose-escalated. Conventionally, the goal of a Phase I trial is to find the maximum tolerated combination (MTC). However, in trials involving an immunotherapy, it is also essential to test whether a difference in toxicities associated with the MTC and the standard of care alone is present. This information can give useful insights about the interaction of the compounds and can provide a quantification of the additional toxicity burden and therapeutic index. We show that both, testing for difference between toxicity risks and selecting MTC can be achieved using a Bayesian model-based dose-escalation design with two modifications. Firstly, the standard of care administrated alone is included in the trial as a control arm and each patient is randomized between the control arm and one of the combinations selected by a model-based design. Secondly, a flexible model is used to allow for toxicities at the MTC and the control arm to be modeled directly. We compare the performance of two-parameter and four-parameter logistic models with and without randomization to a current standard of such trials: a one-parameter model. It is found that at the cost of a small reduction in the proportion of correct selections in some scenarios, randomization provides a significant improvement in the ability to test for a difference in the toxicity risks. It also allows a better fitting of the combination-toxicity curve that leads to more reliable recommendations of the combination(s) to be studied in subsequent phases.

KEYWORDS:

• Dose-escalation
• drugs combination
• immunotherapy
• nonmonotonic
• phase i clinical trial
• randomization

1. Introduction

1.1. Background

Immunotherapies have been shown to have strong antitumor activity and better safety profile than cytotoxic cancer therapies (Disis, 2010; Kyi and Postow, 2014). It was also found that a therapeutic effect of immunotherapy can be noticeably improved by combining it with a backbone chemotherapy or other targeted agents (Sharma and Allison, 2015). Therefore, there are many Phase I combination trials investigating either (i) the added value of an immune checkpoint blocker to a backbone therapy or (ii) the added value of a new agent to an immune checkpoint blocker (Pardoll, 2012). In both cases, one agent, called a standard of care (SoC), is administered at a fixed dose and another agent is dose escalated. Such setting allows reducing the combination context to the context of a single agent trial. However, the application of well-established single agent designs can hide potential difficulties and might fail to address the needs specific to the combination trial.

For decades, the main paradigm of cancer dose-escalation trials was “the more the better” – an agent has greater activity but greater toxicity as dose increases. Relying on this assumption of monotonically increasing toxicities, the identification of the maximum tolerated dose (MTD), the key objective of a single agent Phase I trial, can be efficiently achieved without characterizing the true dose-toxicity relationship (O’Quigley and Shen, 1996). While the monotonicity assumption holds for cytotoxic agents (e.g. chemotherapy and radiotherapy), it can be violated for molecularly targeted agents (MTA), for example, an immunotherapy. Specifically, the majority of Phase I trials of immune-checkpoint protein blockers, such as anti-programmed-death-receptor-1 (PD1), have never actually reached the MTD (see, e.g. Patnaik et al., 2015; Robert et al., 2014; Yamamoto et al., 2017). Consequently, when considering combination trials involving an immunotherapy, the information beyond the conventional (and only) objective of Phase I studies, to find the maximum tolerated combination (MTC), can be also crucial for a more accurate choice of the combination(s) to be studied in subsequent phases.

In this communication, we focus on the setting that covers two important types of clinical trials with the following scientific questions added to the MTC selection:

1. The SoC, a cytotoxic agent, is given at the full single agent dose and an immune-checkpoint blocker is dose escalated.

   1. A clinician would tolerate only a slight increase in toxicity, $\tau$, over the toxicity of the SoC (Paller et al., 2014). Then, the question “is the increase in toxicity acceptable?” should be tested.

2. The SoC is MTA and is given at the full single agent dose and another agent (either a cytotoxic agent or an MTA) is dose escalated.

   1. A clinician would like to test whether there is a plateau region (defined by not exceeding the difference in associated toxicities by more than $\tau$) between the estimated MTC and the SoC. This information would help to define the therapeutic index more accurately.

   2. A clinician has an expectation of additional toxicity, $\tau$, over the SoC under the assumption of compounds independence and would like to test for a possible interaction effect (defined as an additional toxicity over $\tau$).

The primary objective of these trials is to find the MTC. The secondary objectives are to quantify the expected difference between single and combination treatments and to determine the shape of the combination-toxicity relationship. We show that standard single agent dose-escalation methods currently used for such trials may fail to address the secondary questions of interest.

To achieve all objectives, we propose two modifications that are commonly used outside of oncology dose-escalation trials to a Bayesian model-based design. We propose to include the SoC given alone as a control arm in Phase I clinical trial and, to randomize each patient to the control arm or to one of the combinations selected by a Bayesian model-based design. We demonstrate that such randomization procedure allows for a reliable statistical evaluation of additional toxicity (added over the SoC) and of general interaction effects between compounds (Hellmann et al., 2017). We also show that the Emax model provides a well-established tool to detect and to evaluate different patterns in a dose-toxicity relation and its parameters match directly the information needed to address all stated objectives.

It is important to mention that the majority of recently proposed dual agent combination dose-escalation designs allow both agents to be dose-escalated (see, e.g. Ivanova and Wang, 2004; Lin and Yin, 2017; Wages et al., 2011; Yin and Yuan, 2009) However, in the considered types of trials, one can benefit from the fact that the dose of one agent is fixed and it is, in fact, a unidimensional MTC search. Therefore, the proposals of this work are built upon the single-agent model-based approaches.

The remainder of the manuscript proceeds as follows. In Section 1.2, we present two motivating examples. In Section 2, the Bayesian model-based approach and the Emax model are described together with the proposed randomization scheme. In Section 3, operating characteristics of the proposed design are compared to currently used dose-escalation designs in various scenarios. In Section 4, the problem of fitting the four-parameter model and the impact of the randomization ratio are investigated. The paper concludes with a discussion in Section 5.

1.2. Motivating trials

This work is motivated by two recent combination trials that could potentially benefit from the implementation of the proposed design.

1.2.1. Backbone chemotherapy and increasing doses of immunotherapy

Gemcitabine is the standard chemotherapy for advanced pancreatic cancer (Aglietta et al., 2014) which has a narrow therapeutic index (Crane et al., 2002). Tremelimumab is a fully humanized monoclonal antibody that may allow and potentiate effective immune responses against tumor cells. In several clinical studies, anti-CTLA4 agents have been shown to induce durable tumor responses through modulation of the immune system in patients with metastatic melanoma. The hypothesis was that the combination of the two agents “might provide synergistic antitumor activity without increasing toxicity” (Aglietta et al., 2014). Gemcitabine 1000 mg/m² was administered in all patients while escalating doses of tremelimumab (6, 10, or 15 mg/kg) were sequentially tested to identify the MTC. This trial was conducted using the 3 + 3 design, leading to a high risk of erroneous conclusions (Reiner et al., 1999) due to the assumption of increasing dose-toxicity relationship. Furthermore, as the expected risk of the dose-limiting toxicity with Gemcitabine is more than 20%, any (if the dose-response curve is flat) or no combination (if the dose-response curve is steep) can be the MTC. Finally, the highest combination was recommended for further investigations, but the question whether the toxicity of the MTC is actually increased compared to Gemcitabine given at 1000 mg/m² was never formally tested.

1.2.2. Fixed immunotherapy and increasing doses of standard treatment

Sorafenib is a standard treatment for advanced cellular cell carcinoma. However, its efficacy remains limited as the time to progression is around six months. Despite this agent being prescribed at the MTD, its therapeutic index makes it possible to reduce the dose in case of adverse reactions (Wilhelm et al., 2006). SPLASH is a dose-escalation study of the safety of Avelumab in combination with sorafenib in patients with advanced cellular cell carcinoma that is about to be initiated at Gustave Roussy. While Avelumab will be given at a fixed dose, Sorafenib will be escalated from 200 mg/kg up to 800 mg/kg that is the MTD when administered as a single agent. Dose escalation will be driven by safety using the Bayesian Continual Reassessment Method (O’Quigley et al., 1990, CRM). The MTC is defined as the highest tested combination associated closest to 25% probability of causing a dose-limiting toxicity during cycle 1. The expected rate of toxicity with Avelumab alone is 8%.

2. Methods

2.1. Setting

Consider a clinical trial in which combinations of agents $A$ and $B$ are studied. The drug $A$ is a standard of care (SoC) given at a fixed dose, $a$, and $B$ is dose escalated. We specify increasing doses of $B$: $0=b_0<b_1<b_2<\dots <b_m$. Then, ${\tilde{d}}_0=\{a,b_0\},{\tilde{d}}_1=\{a,b_1\},\dots ,{\tilde{d}}_m=\{a,b_m\}$ are $m+1$ combinations available in the trial, where ${\tilde{d}}_0$ corresponds to the agent $A$ given as a single agent and is subsequently referred to as the control arm and the combinations ${\tilde{d}}_1,\dots ,{\tilde{d}}_m$ are referred to as the investigational arm. An investigator observes binary outcomes of a random variable $Y$, $y=0$ if no toxicity is observed and $y=1$, otherwise. Let $p_j$ be the probability for a patient to experience a toxicity given the combination ${\tilde{d}}_j$. It is assumed that toxicity is a nondecreasing function of combination, $p_0\le p_1\le \dots \le p_m$ and a reliable prior information for the toxicity probability of the control arm, $p_0$, is available.

Let $\gamma$ be the maximum acceptable toxicity probability. The primary goal of the dose-escalation clinical trial is to find the MTC ${\tilde{d}}_{j^{\star }}$ such that

$j^{\star }=\underset{j=0,\dots ,m}{\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}|p_j-\gamma |$

using estimated toxicity probabilities ${\hat{p}}_0,\dots ,{\hat{p}}_m$. The secondary goal, specific to the combination trials, is to test if the difference of toxicity probabilities associated with the estimated MTC and the SoC is as expected. To unify notation, we study the extra toxicity (ET) defined in terms of the expected difference in toxicity, $\tau$, between the MTC and the SoC. One would conclude the ET when the probability that the difference in the two toxicity risks exceeds $\tau$ is larger than some level $\alpha$ given observed data. Formally, if $\mathbb{P}\left(p_{j^{\star }}-p_0\ge \tau |{\mathbb{Y}}_n\right)\ge \alpha$ where ${\mathbb{Y}}_n$ is the vector of outcomes observed in a trial and the probability is found with respect to the posterior probability distribution, then the ET is concluded. However, the observed outcomes of the trial are random, and one would like to be sure that this difference can be found regardless these outcomes. Therefore, if the ET is present one would like the probability

(1) $\mathcal{P}\equiv \mathbb{P}\left(\mathbb{P}\left(p_{j^{\star }}-p_0\ge \tau |{\mathbb{Y}}_n\right)>\alpha \right)$(1)

to be equal to 1. In contrast, if there is no ET, the probability (1) is desired to be equal to 0.

2.2. Bayesian model-based design

In the considered type of trials, one compound is escalated while another one is kept fixed. Therefore, the underlying combinations can be considered as doses in a single agent trial and corresponding dose-escalation designs can be adapted. Particularly, we use the Bayesian continual reassessment method (O’Quigley et al., 1990, CRM) in this work.

Assume that the toxicity probability has the functional form

$p_j=\psi (d_j,\theta )$

where $\theta \in {\mathbb{R}}^h$ is a $h$-dimensional vector of parameters, ${\mathbb{R}}^h$ is a real coordinate h-dimensional space and $d_j$ is a scalar unit-less combination level (also referred to as a standardized level) corresponding to the same toxicity probability as combination ${\tilde{d}}_j$. The latter adjustment is needed to adapt the two-dimensional combination to the single-agent context.

Denote prior distributions of $\theta$ by $f_0(.)$. Assume that $n$ patients have been already assigned to doses $d(1),\dots ,d(n)$ and binary responses ${\mathbb{Y}}_n=[y_1,\dots ,y_n{]}^{\mathrm{T}}$ were observed, respectively. The CRM updates the posterior distribution of $\theta$ using Bayes’s Theorem

(2) $f_n(\theta )=\frac{f_{n-1}(\theta )\varphi (d(n),y_n,\theta )}{{\int }_{{\mathbb{R}}^h}{f}_{n-1}(u)\varphi (d(n),y_n,u)\mathrm{d}u}=\frac{f_0(\theta )\prod _{i=1}^n\varphi (d(i),y_i,\theta )}{{\int }_{{\mathbb{R}}^h}{f}_0(u)\prod _{i=1}^n\varphi (d(i),y_i,u)\mathrm{d}u}$(2)

where

$\varphi (d(n),y_n,\theta )=\psi (d(n),\theta {)}^{y_n}(1-\psi (d(n),\theta ){)}^{1-y_n}.$

Then, the posterior mean of the toxicity probability for combination ${\tilde{d}}_j$ after $n$ patients is equal to

(3) ${\hat{p}}_j^{(n)}=\mathbb{E}(\psi (d_j,\theta )|{\mathbb{Y}}_n)={\int }_{{\mathbb{R}}^h}\psi (d_j,u)f_n(u)\mathrm{d}u.$(3)

The combination ${\tilde{d}}_k$ that minimizes

(4) $\left|{\hat{p}}_j^{(n)}-\gamma \right|$(4)

among all ${\tilde{d}}_0,\dots ,{\tilde{d}}_m$ is recommended for the next group of patients. The procedure is repeated until the maximum number of patients, $N$, has been treated. While it was shown that the CRM is able to select the MTC with a high probability, it is not designed to account for the specific of the combination trial. Therefore, we propose the following modifications.

2.3. Randomization between control and investigational arms

Patients in dose-escalation trials are assigned cohort-by-cohort, where a cohort is a small group of typically 1–4 patients. By minimizing the criterion (4), the CRM design tends to assign the majority of patients in the neighborhood of the MTC. This leads to a sparse allocation of patients on other combinations and on the control arm. This will make it difficult to test for the difference in toxicity risks associated with the estimated MTC and the control arm. Following applications outside of oncology, we introduce the following randomization procedure between the investigational (combinations) and control (SoC) arms.

Denote the cohort size by $c=c_1+c_2$ and let $c_1$ be the number of patients in cohort assigned to the estimated MTC (according to criterion (4) by the CRM) and $c_2$ be the number of patients in cohort assigned to the control arm, ${\tilde{d}}_0$. This results in at least $N_2=\frac{c_2}{c_1+c_2}N$ patients on the control arm and at most $N_1=\frac{c_1}{c_1+c_2}N$ on the investigational arm by the end of the trial. For instance, taking $c_1=3$ and $c_2=1$ (denoted by 3:1), one will end up with at least $25\mathrm{\%}$ of the total sample size being assigned to the control. Note that the model-based design is allowed to select the control arm as the estimated MTC if the associated toxicity is the closest to the target $\gamma$. This avoids exposing patients to highly toxic combinations if the first combination has a toxicity probability beyond $\gamma$. Therefore, the total number of patients on the control arm can be more than $\frac{c_2}{c_1+c_2}N$.

The modified allocation rule (the randomization of patients between control and investigational arms) raises an important question of the suitable choice of the parametric model $\psi (d_j,\theta )$ which is discussed below.

3. 2.4 4-parameter Emax model

In the context of the single agent trial, the one-parameter power model

(5) $\psi (d_j,{\theta }_0)=d_j^{{\theta }_0}$(5)

is argued to be an appropriate choice for the working dose-toxicity model (O’Quigley et al., 1990; O’Quigley and Shen, 1996; Iasonos et al., 2016), where ${\theta }_0$ is a scalar parameter. There are, however, at least two reasons why such (and similar one-parameter models) are not suitable for the considered type trials.

Firstly, the parametric form (5) implies the “strict” monotonicity assumption – the toxicity increases with the dose. It does not allow modeling the plateau in the combination-toxicity relations that is relevant for trials with immunotherapies. Moreover, one can expect that such a model would almost always conclude the ET regardless the scenario which is undesirable in the considered combination trials.

Secondly, the main argument behind using the one-parameter model (and, consequently, the main critique of using many parameter models) is that the CRM design tends to collect observations in the neighborhood of the MTC only (Iasonos et al., 2016). This means that the approximation of the dose-toxicity relation in the neighborhood of one point only is of interest and the one-parameter model is rich enough to achieve it. Given the modified allocation rule above, the majority of patients are now assigned in the neighborhood of two points: the control arm and the estimated MTC. Therefore, more flexible models are essential to consider.

Finally, in the considered types of trials a clinician is interested in one more question beyond the MTC and ET identification. Particularly, one would like to estimate the combination-toxicity relationship to allow for a more reliable recommendation for the Phase II clinical trial.

These specifics of the trials with immunotherapy translate in the following conditions on the combination-toxicity relationship:

1. A model should allow for the toxicity risk at the SoC to be modeled independently – the toxicity of the SoC does not depend on the interaction between compounds.

2. A model should allow for modeling a plateau at a particular level – a plateau is a plausible shape given dose finding studies of immunotherapies.

3. A model should allow for modeling a steep increase in toxicity risk in several combinations – the interaction effect of the backbone therapy can appear for a particular dosing of a complementary agent.

There is a number of the dose-response functions that are conventionally considered in dose finding studies in various disease areas and satisfy some of the properties above. For example, a two-parameter logistic model, which is used in oncology dose-escalation trials (Neuenschwander et al., 2008; Whitehead and Williamson, 1998), allows for the toxicity risk at the SoC to be modeled independently hence satisfies property (1). A three-parameter Emax model commonly used in dose-ranging studies (see, e.g. Bornkamp, 2017, for review of various dose-response shapes) allows for a plateau and modeling toxicity at SoC independently and hence satisfies properties (1) and (2). It can, however, model the steep increase in toxicities (property (3)) only at the beginning of the combination-toxicity relationship, but not for moderate combinations.

A potential candidate model satisfying all three properties that we propose to use is the 4-parameter Emax model

$\psi (d_j,E_0,E_{max},\lambda ,E{D}_{50})=E_0+\frac{d_j^{\lambda }{E}_{max}}{d_j^{\lambda }+E{D}_{50}^{\lambda }},$

where $E_0$ is the toxicity probability on the standard treatment, $E_{max}+E_0$ is the maximum toxicity probability attributable to the combination, $E{D}_{50}$ is the combination which produces $E_0+\frac{E_{max}}{2}$ toxicity and $\lambda$ is the slope factor. Following the assumption of the nondecreasing toxicity probability, we specify $\lambda \ge 0$. The model is already adjusted to accommodate dose combinations as it uses the unit-less variable $d_j$. To construct the standardized levels, $d_j$, one needs to represent them in terms of prior estimates of toxicity probabilities ${\hat{p}}_j^{(0)}$ associated with combinations ${\tilde{d}}_j$ $j=0,\dots ,m$

$d_j=\hat{E}{D}_{50}^{(0)}\times {\left(\frac{{\hat{p}}_j^{(0)}-{\hat{E}}_0^{(0)}}{{\hat{E}}_{max}^{(0)}+{\hat{E}}_0^{(0)}-{\hat{p}}_j^{(0)}}\right)}^{\frac{1}{{\hat{\lambda }}^{(0)}}}$

where ${\hat{\lambda }}^{(0)}$, ${\hat{E}}_0^{(0)}$, ${\hat{E}}_{max}^{(0)}$, and $\hat{E}{D}_{50}^{(0)}$ are prior point estimates of the model parameters. The prior point estimate of the toxicity probability on the control arm is ${\hat{p}}_0^{(0)}$. Therefore, by definition, ${\hat{p}}_0(0)\equiv {\hat{E}}_0^{(0)}$ that leads to $d_0=0$. Modeling $E_0$ directly guarantees that the sequential update of other parameters does not contribute to the toxicity probability estimation on the control arm. In this case, the model takes the trivial form

$\psi (d_0,\cdot )=E_0.$

Intuitively, it reflects that the toxicity probability of the standard therapy does not depend on the mechanism of its interaction with compound $B$, yet all data are used for the estimation of $E_0$.

The parameters of the Emax model allows to model the toxicity on the control arm and the plateau and have a clear clinical interpretation. Importantly, Emax model is widely employed in dose-finding studies outside of oncology and being extensively studied in a medical context (Bornkamp, 2017). We compare its performance to other widely employed in dose-finding studies one-parameter and two-parameter models in Section 3.

4. Numerical results

4.1. Simulation setting

We explore the performance of the Bayesian model-based dose-escalation method incorporating randomization to the control arm into the Emax model by simulations in many different scenarios. As before, we consider a trial studying combinations of two agents with the dose for the standard therapy being fixed. Motivated by recent combination trial reviews (see, e.g. Dahlberg et al., 2014; Riviere et al., 2015) we consider a setting with $N=48$ patients and $m=7$ combinations. Following the motivating trials, we set the target toxicity probability $\gamma =0.25$, the expected difference $\tau =0.05$ and the confidence level $\alpha =0.90$.

Three main characteristics, (i) the proportion of correct selections (PCS), (ii) the proportion of times ET is found and (iii) a goodness-of-fit measure, are considered. A goodness-of-fit measure is used to capture the overall shape of the dose-toxicity relationship which is of interest in the considered type of trials. We use the scenario-normalized mean squared error defined as

(6) $\mathrm{N}\mathrm{M}\mathrm{S}\mathrm{E}=\frac{1}{S}\sum _{s=1}^S\sqrt{\frac{{\sum }_{j=0}^7{(p_j-{\hat{p}}_j^{(s)})}^2}{{\sum }_{j=0}^7{(p_j-{\tilde{p}}_j^{(s)})}^2}}$(6)

where $S$ is the number of replications, ${\hat{p}}_j^{(s)}$ is a toxicity probability estimate corresponding to combination ${\tilde{d}}_j$ obtained by the design in the $s^{th}$ simulation and ${\tilde{p}}_j^{(s)}$ is a toxicity probability estimate for combination ${\tilde{d}}_i$ obtained by the non-parametric optimal benchmark approach (O’Quigley et al., 2002) in the $s^{th}$ simulation. The estimates ${\tilde{p}}_j^{(s)}$ are computed employing the complete information profiles as introduced by O’Quigley et al. (2002) and using the same sequences of the patients toxicity profiles as for the design under evaluation. The normalization allows a comparison across different scenarios as the nonparametric optimal estimate allows for the incorporation of the each scenario specificities. One would like the NMSE to be equal to 1 as it corresponds to the curve estimated as good as by the optimal benchmark.

4.1.1. Prior specification

The standardized levels, $d_j$, are constructed using the skeleton ${\hat{\mathbf{p}}}^{(0)}$.

${\hat{\mathbf{p}}}^{(0)}=[0.08,0.25,0.35,0.45,0.55,0.65,0.70,0.75{]}^{\mathrm{T}}.$

The first value, 0.08, corresponds to mean prior toxicity probability for the standard therapy following the motivating trial of Sorafenib. The prior MTC is the first combination ensures the start of the trial at the lowest combination. Other values are obtained by taking the adequate spacing between prior values as the CRM design has been shown to be robust and efficient in this case (O’Quigley and Zohar, 2010).

In contrast to the skeleton that is the same whatever model is used for the combination-toxicity relationship, the prior distributions of the model parameters are calibrated so that each design carries the same amount of information. To ensure that the proposed approach and competitive designs are evaluated under comparable set-ups, the credible intervals for the prior toxicity probabilities associated with the control arm and the prior MTC are used. Taking into account that one usually has some reliable information about the SoC, but a limited information about combinations, we specify prior distribution to satisfy the following assumptions:

1. Control arm ${\tilde{d}}_0$: the expected toxicity probability ${\hat{p}}_0^{(0)}=0.08$ and the upper bound of the $95\mathrm{\%}$ credibility interval is 0.25.

2. Prior MTC ${\tilde{d}}_1$: the expected toxicity probability ${\hat{p}}_1^{(0)}=0.25$ and the upper bound of the $95\mathrm{\%}$ credibility interval is 0.80.

Following these conditions the prior distributions of the Emax model were specified as

(7) $E_0\sim \mathbb{B}(0.8,10-0.8),E_{max}|E_0\sim \mathbb{U}[0,1-E_0],E{D}_{50}\sim \mathrm{\Gamma }(0.4,0.4),\lambda \sim \mathrm{\Gamma }(1,1)$(7)

where $\mathbb{B}(u,v)$ denotes the Beta distribution with parameters $u,v$ and $\mathrm{\Gamma }(u,v)$ denotes a Gamma distribution with mean $\frac{u}{v}$ and the variance $\frac{u}{v^2}$. As prior distribution may have quite an impact on the parameter estimates, we first rely on the elicitation given in (7) but also investigate the robustness to the priors in Section 4.1. Since an investigator usually has a reliable prior for $E_0$ and the detection of a plateau should be determined by data alone (hence an uninformative prior on $E_{max}$ is imposed), we fix the prior distributions of $E_0$ and $E_{max}$ and consider less informative prior distributions of $E{D}_{50}$ and more informative one for $\lambda$

(8) $E{D}_{50}\sim \mathrm{\Gamma }(0.25,0.25),\lambda \sim \mathrm{\Gamma }(3,3)$(8)

and more informative prior for $E{D}_{50}$ and less informative for $\lambda$

(9) $E{D}_{50}\sim \mathrm{\Gamma }(1,1),\lambda \sim \mathrm{\Gamma }(0.05,0.05).$(9)

We would refer to the design using the Emax model with randomization to the control arm as “EmaxR.”

4.1.2. Cohort size for the control arm and combination skipping

The cohort size is fixed to be $c_1=3$ for the investigational arm and $c_2=1$ for the control arm. Therefore, at least 12 patients ($25\mathrm{\%}$) will be assigned to the standard therapy and at most 36 patients are assigned to the investigational arm. Other randomization ratios, $2:1$ and $4:1$, are studied in Section 4.2. We allow combinations to be skipped and evaluate the impact of this constraint in Section 4.2.

4.1.3. Scenarios

Nine qualitatively different scenarios, depicted in Figure 1, are used to study the performance of the novel approach. Scenarios 1–5 are consistent with the prior information of the control arm toxicity probability. Scenario 1 corresponds to the commonly used logistic dose-toxicity shape while scenario 2 considers a flat dose-toxicity curve with the MTC being ${\tilde{d}}_7$ and no additional toxicity beyond $\tau$. In scenario 2, the model monotonicity assumption implies that the toxicity probability on the control has the largest difference with combination ${\tilde{d}}_7$ and the model is expected to have a tendency to claim the ET. Scenario 2 is used to test if the randomization to the control group and the Emax model can prevent false conclusions regarding the ET. Scenarios 3–4 have a steep dose-toxicity relationship that plateaus around the target toxicity risk and outside the maximum acceptable toxicity probability, respectively. These scenarios are used to investigate under what conditions it is easier for the design to find the ET. Scenario 5 also corresponds to a sharp increase in the toxicity probability, but for medium combinations which have unacceptable high toxicity such that there is no ET between the MTC and the control arm.

Figure 1. Considered combination-toxicity shapes. The MTC is marked by a triangle.

Scenarios 6–9 reflect the case of a misspecified prior toxicity probability for the control arm. In scenario 6, the true toxicity probability for ${\tilde{d}}_0$ is below the prior value. While the prior information dictates that there is no ET between the MTC and the control ($0.11-0.08<\tau =0.05$), although in reality, it is present. In scenarios 7–9, the toxicity probability for the control is underestimated by the prior information and the control arm has a toxicity probability close to the maximum acceptable toxicity. Both cases in which the ET is present (scenario 8) and is absent (scenarios 7 and 9) are considered.

4.1.4. Competing models

The performance of the proposed approach is compared to the design which is currently used in the considering type of Phase I clinical trial – the CRM using the 1-parameter power model (O’Quigley and Shen, 1996, P1)

$\psi (d_j,{\theta }_0)=d_j^{{\theta }_0}$

without randomization, where ${\theta }_0$ is an unknown scalar parameter. We also consider the two-parameter logistic model (Whitehead and Williamson, 1998)

$\psi (d_j,{\theta }_1,{\theta }_2)=\frac{exp(log({\theta }_1)+{\theta }_2{d}_j)}{1+exp(log({\theta }_1)+{\theta }_2{d}_j)}$

without randomization (L2) and with randomization (L2R), where ${\theta }_1$ and ${\theta }_2$ are unknown scalar parameters. As there are at most 36 patients on the investigational arm in the randomized setting against 48 in the nonrandomized one, we would also consider the two-parameter logistic model without randomization and $N=36$. This would allow to spot the influence of the randomization on the operating characteristics. The following prior distributions for the parameters of P1, L2, L2R models ${\theta }_0\sim \mathrm{\Gamma }\left(5.5,3\right)$ and

$\left(log({\theta }_1),log({\theta }_2)\right)\sim \mathcal{N}\left({\left(log\left(\frac{0.08}{0.92}\right)-0.75,0.675\right)}^{\mathrm{T}},\left[\begin{array}{cc}1.5& 0\\ 0& 0.75\end{array}\right]\right)$

are chosen to satisfy two conditions on ${\tilde{d}}_0$ and ${\tilde{d}}_1$ formulated above.

The characteristics of all the models compared are evaluated in R (R Core Team, 2015) using the bcrm-package by Sweeting et al. (2013). To accommodate the randomization to the control group and the Emax model corresponding modifications to the package were made. The code is available online at https://github.com/dose-finding/Randomised-bcrm.

4.2. Operating characteristics

4.2.1. Proportion of correct selections

The results of the comparison are summarized in Table 1 (scenarios 1–5) and Table 2 (scenarios 6–9), with each figure representing the proportion of each combination selections. The last columns represents the proportion of toxic responses.

Table 1. Operating characteristics of the EmaxR, L2R, L2, and P1 models in scenarios 1–5. The figures correspond to the proportion of each combination selections. The MTC selection is in boldface. ‘TR’ corresponds to the average proportion of toxic responses. Results are based on ${10}^4$ simulations.

	${\tilde{d}}_0$	${\tilde{d}}_1$	${\tilde{d}}_2$	${\tilde{d}}_3$	${\tilde{d}}_4$	${\tilde{d}}_5$	${\tilde{d}}_6$	${\tilde{d}}_7$	TR
Scenario 1
Toxicity	0.08	0.10	0.12	0.15	0.25	0.40	0.45	0.47
EmaxR	0.0	2.1	8.8	24.9	43.1	13.6	3.8	3.7	19.1
L2R	0.0	1.8	8.3	26.5	44.6	13.8	3.1	2.0	19.2
L2 ($N=36$)	0.4	0.6	4.4	29.7	45.4	15.6	2.4	1.5	25.2
L2	0.4	0.1	3.5	24.9	52.6	16.4	2.1	1.3	25.1
P1	0.0	1.0	4.2	16.5	51.4	20.4	5.5	1.0	29.2
Scenario 2
Toxicity	0.0800	0.0900	0.0950	0.1000	0.1050	0.1100	0.1125	0.1150
EmaxR	0.0	0.3	0.5	0.9	0.5	1.6	1.9	94.4	10.2
L2R	0.0	0.6	0.8	1.5	1.8	1.6	2.5	91.2	10.1
L2 ($N=36$)	0.1	0.3	0.8	0.3	0.9	1.1	2.0	94.5	11.0
L2	0.3	0.1	0.3	0.5	1.8	1.0	0.5	95.6	11.4
P1	0.0	0.0	0.0	0.0	0.0	0.6	0.8	98.7	11.3
Scenario 3
Toxicity	0.08	0.25	0.26	0.265	0.27	0.27	0.27	0.27
EmaxR	2.4	44.3	13.9	7.3	7.1	5.0	4.1	15.9	21.2
L2R	2.8	44.7	15.8	10.1	7.6	4.4	2.5	12.1	21.3
L2 ($N=36$)	5.3	25.5	13.4	10.1	10.6	7.53	5.7	21.9	25.3
L2	5.1	27.0	13.0	9.4	11.5	7.1	6.1	20.7	25.3
P1	0.0	0.0	0.5	0.5	13.4	18.9	17.9	44.7	26.6
Scenario 4
Toxicity	0.080	0.250	0.500	0.510	0.520	0.525	0.530	0.535
EmaxR	5.1	82.1	11.5	0.9	0.3	0.1	0.0	0.0	24.4
L2R	4.13	86.4	8.4	0.9	0.3	0.0	0.0	0.0	24.6
L2 ($N=36$)	20.7	72.0	6.1	0.7	0.3	0.1	0.0	0.0	28.4
L2	20.90	75.0	4.0	0.1	0.0	0.0	0.0	0.0	26.9
P1	20.2	71.7	7.4	0.7	0.0	0.0	0.0	0.0	31.0
Scenario 5
Toxicity	0.080	0.085	0.090	0.100	0.50	0.50	0.50	0.50
EmaxR	0.0	1.9	7.8	57.1	31.2	1.6	0.1	0.4	18.8
L2R	0.0	1.5	9.5	56.4	30.2	1.4	0.8	0.3	18.2
L2 ($N=36$)	0.0	0.0	4.1	59.7	31.7	1.5	0.5	0.4	26.1
L2	0.1	0.1	1.4	62.7	35.5	0.2	0.1	0.0	25.1
P1	0.0	0.0	0.0	54.9	36.4	7.0	1.2	0.5	33.2

Table 2. Operating characteristics of the EmaxR, L2R, L2, and P1 models in scenarios 6–9. The figures correspond to the proportion of each combination selections. The MTC selection is in boldface. ‘TR’ corresponds to the average proportion of toxic responses. Results are based on ${10}^4$ simulations.

	${\tilde{d}}_0$	${\tilde{d}}_1$	${\tilde{d}}_2$	${\tilde{d}}_3$	${\tilde{d}}_4$	${\tilde{d}}_5$	${\tilde{d}}_6$	${\tilde{d}}_7$	TR
Scenario 6
Toxicity	0.020	0.030	0.050	0.080	0.090	0.100	0.110	0.115
EmaxR	0.0	0.0	0.0	0.0	0.4	1.1	1.2	97.2	8.3
L2R	0.0	0.0	0.0	0.0	1.0	1.1	1.9	96.1	8.7
L2 ($N=36$)	0.0	0.0	0.0	0.0	0.6	0.9	0.8	97.5	10.9
L2	0.0	0.0	0.0	0.1	0.1	0.2	0.5	99.3	11.3
P1	0.0	0.0	0.0	0.0	0.0	0.1	1.1	98.8	11.3
Scenario 7
Toxicity	0.20	0.23	0.40	0.40	0.40	0.40	0.40	0.40
EmaxR	12.7	64.3	15.1	4.7	1.7	0.7	0.5	0.4	25.8
L2R	9.4	67.2	15.7	5.3	1.4	0.3	0.4	0.3	25.8
L2 ($N=36$)	19.0	52.0	18.3	4.5	2.0	1.7	1.4	1.1	30.1
L2	15.0	62.5	15.1	3.6	1.8	0.8	0.8	0.5	29.2
P1	0.2	28.0	29.2	23.5	14.1	3.87	0.8	0.4	35.2
Scenario 8
Toxicity	0.17	0.20	0.25	0.35	0.40	0.40	0.40	0.40
EmaxR	5.1	40.3	29.5	16.1	5.9	1.7	0.8	0.5	23.6
L2R	2.8	45.3	31.5	12.7	4.9	1.8	0.8	0.20	23.3
L2 ($N=36$)	4.5	27.5	34.3	19.1	7.3	3.1	1.5	2.7	29.0
L2	4.5	24.9	41.1	19.1	6.5	1.9	1.1	0.9	28.6
P1	0.0	3.4	30.2	37.5	18.7	7.27	1.7	1.1	32.3
Scenario 9
Toxicity	0.25	0.255	0.26	0.265	0.27	0.27	0.27	0.27
EmaxR	21.6	37.6	10.0	7.9	4.7	4.3	2.3	11.6	25.5
L2R	19.6	43.3	12.9	7.2	5.2	2.4	2.1	7.3	26.1
L2 ($N=36$)	10.1	20.7	12.5	10.5	10.1	6.3	5.8	24.0	26.0
L2	10.1	22.9	12.7	10.3	9.5	8.0	6.4	20.0	26.3
P1	0.3	6.5	15.0	26.7	28.1	12.8	6.9	3.7	26.6

Overall, both randomized models, EmaxR and L2R, result in the same PCS with the difference no more than 4% across all scenarios. The models that do not use randomization have a comparable performance in scenarios 1, 2, 4, 6. In the rest of scenarios, P1 is less conservative and selects toxic combinations more often than L2. The nonrandomized designs select the MTC more often if the MTC is located in the middle of the curve – scenarios 1, 5, 8. Comparing L2R and L2 ($N=48$) in these scenarios, L2 outperforms its randomized version by 6–10%. The decrease in the PCS is mainly explained by the reduced number of patients on the investigational arm: the differences between L2R and L2 $(N=36)$ vary in the range 1–3%. While randomization entails a minor drop in these scenarios, it might also result in a notable increase in the PCS if the MTC located at the beginning of the curve, close to the control (scenarios 3, 4, and 7). Considering L2R and L2, the randomization leads to more accurate MTC selection by 17%, 11%, and 5%, respectively. These difference increase to 19%, 14%, and 15% compared to the L2 $(N=36)$ with the same number of patients on the investigational arm. Moreover, the inclusion of the randomization allows decreasing the proportion of toxic responses by up to 6% in scenario 1 that results in nearly two fewer patients with toxic outcomes. Note that misspecified prior distribution for the control (scenarios 6–9) does not affect the MTC selection significantly

In scenarios 2 and 6, the MTC is the highest combination and all models select the MTC in more than 90% of replications with L2R having the least PCS: 91.2% in scenario 2 and 96.1% in scenario 6. In the case when the control treatment is already associated with a maximum acceptable toxicity (scenario 9), the randomized designs selects the control by nearly 10% more often than L2. Generally, the randomized designs are more conservative and select combinations in the beginning of the curve more often with nearly 70% of ${\tilde{d}}_0,{\tilde{d}}_1,{\tilde{d}}_2$ selection against 45% for L2 and 22% for P1.

4.2.2. Probability to find the ET

While both four- and two-parameters models with randomization were shown to have comparable performance in terms of the PCS, the major differences can be found in Table 3. The upper line represents the proportion of times the ET is found ($\hat{\mathcal{P}}$) while the lower line shows the NMSE.

Table 3. Operating characteristics of EmaxR, L2R, L2, and P1 models. The upper line: the proportion of times the ET is found. The lower line: NMSE. Scenarios with no ET are underlined. Results are based on ${10}^4$ simulations.

		$\mathrm{S}c1$	$\underset{\_}{\mathrm{S}c2}$	$\mathrm{S}c3$	$\mathrm{S}c4$	$\underset{\_}{\mathrm{S}c5}$	$\mathrm{S}c6$	$\underset{\_}{\mathrm{S}c7}$	$\mathrm{S}c8$	$\underset{\_}{\mathrm{S}c9}$
EmaxR	$\hat{\mathcal{P}}$	74.7	16.3	66.7	71.5	67.9	24.2	12.2	29.3	6.6
	NMSE	1.7	2.0	2.2	3.5	3.1	1.4	2.5	1.7	3.7
L2R	$\hat{\mathcal{P}}$	71.8	14.4	59.7	64.1	75.0	20.9	21.3	31.4	11.8
	NMSE	2.0	2.2	7.2	4.4	3.4	1.5	4.4	3.1	7.2
L2 ($N=36$)	$\hat{\mathcal{P}}$	47.0	13.7	35.8	35.3	43.5	13.7	25.1	21.3	32.3
	NMSE	1.8	2.0	6.3	5.1	3.0	1.4	6.3	4.0	8.9
L2	$\hat{\mathcal{P}}$	61.5	15.8	50.1	50.2	64.7	18.1	37.5	30.2	43.5
	NMSE	2.0	2.5	7.7	5.8	3.6	1.6	7.6	4.5	9.0
P1	$\hat{\mathcal{P}}$	99.9	93.6	99.7	99.7	99.4	95.3	99.8	99.9	99.7
	NMSE	2.1	2.6	7.8	6.1	3.8	2.1	7.8	4.5	9.2

Comparing randomized designs, EmaxR results in the greater proportion of correctly identified ET in scenarios 1, 3, 4, 6, and 8 by 2–7%. In scenarios with no ET (2, 5, 7, and 9), EmaxR declares the ET less often than L2R by 2–9%. Note, however, that both approaches wrongly declare the ET in the majority of trials in scenario 5. Both models struggle to captured the sharp jump in the toxicity curve and overestimate the toxicity probability for the MTC (see also Figures 2 and 3 for the illustration of fitted curves). It is also challenging for both approaches to determine the ET in scenarios 6 and 8 when the prior for the control is misspecified. EmaxR correctly identifies the ET in 24% and 30% of all trials. At the same time, the prior misspecification is not being an issue in scenarios 7 and 9 with no ET as the EmaxR finds additional toxicity in only 12% and 6% of trials, respectively.

Figure 2. Mean values (solid lines) and 90% credible intervals (shadowed area) of the toxicity probabilities for each combination obtained by fitted Emax (green), L2R (blue), and L2 (red) models in scenarios 1–4. The true toxicity curve is marked by the dashed-dotted line. Results are based on 10⁴ replications.

Figure 3. Mean values (solid lines) and 90% credible intervals (shadowed area) of the toxicity probabilities for each combination obtained by fitted Emax (green), L2R (blue), and L2 (red) models in scenarios 5–9. The true toxicity curve is marked by the dashed-dotted line. Results are based on 10⁴ replications.

Regarding the design with no randomization, P1 almost always declares ET regardless the underlying scenario. As expected the strong assumption of monotonicity imposed by the power model results in this (often) incorrect conclusion. L2 does not share the same pattern and it is able to declare the ET with low probability if it is not present (e.g. scenario 2). At the same time, the inclusion of the randomization in L2 leads to large increase in the proportion of correctly identified ET by 10–14% in scenarios 1, 3, and 4. Again, comparing approaches with the same number of patients on the investigational arm these proportions differ by 25–30%. In the rest of scenarios with the ET (scenarios 6 and 8), L2 and L2R have the same characteristics with L2R outperforming the L2 ($N=36$) by 5–10%. In the majority of scenarios with no ET, using L2R instead of L2 leads to a decrease in the proportion of incorrect conclusions of the ET, for example, by 16% and 32% in scenarios 7 and 9, respectively. There is only one example when randomization increases the probability to find the ET while it is not there – scenario 5.

4.2.3. Goodness-of-fit

Comparing the ability of models to fit the toxicity curves, being the most flexible model EmaxR results in the smallest values of NMSE in all scenarios with the largest difference in scenario 3 $-$ 2.2 and 7.8 (for P1), and the smallest difference in scenario 6 $-$ 1.4 and 1.5 (for L2R). Similarly, P1, the least flexible alternative, results in the greatest values of the NMSE in all scenarios. L2 shows a better fit than P1, but it can be further improved by the randomization. Low values of NMSE for L2 with the reduced number of patients can be explained by the decreased accuracy of the nonparametric optimal benchmark.

Further differences in the ability of the models to detect a dose effect are given in Figures 2 and 3 presenting the mean values of toxicity probabilities and corresponding 90% credible interval. The mean values for estimated probabilities are given by solid lines while the credible interval is given by shadow areas of the corresponding color. The true toxicity probabilities are given by dashed-dotted lines.

EmaxR corresponds to the best fit of the toxicity curve and to the narrowest credible interval among all alternatives in all scenarios. The fitted curves of L2R and L2 are similar for the first combinations, but the probability values obtained by L2R are more accurate due to the shift toward the true toxicity probability for further combinations.

Overall, it is found that EmaxR is able to give a more accurate answer to secondary questions of the Phase I trial about the ET and a possible dose effect. The cost of this answer is a minor drop in the PCS in scenarios with the MTC being in the middle of the combination-toxicity curve. Comparing the same parametric models, it was shown that the randomization is generally allowed to increase the proportion of correctly determining if the ET is present and also improves the fit of the curve. At the same time, the 1-parameter power model does not meet the requirement to answer secondary goals of the study as it almost always declares the ET.

5. Sensitivity analysis

5.1. Prior distributions

In this section the influence of different sets of prior distributions is studied. We compare the previously used prior distributions given in (7) with (8) denoted by EmaxR ($\lambda$) and the set given in (9) and denoted by EmaxR ($E{D}_{50}$). The summary of the simulation results are given in Table 4.

Table 4. Operating characteristics of EmaxR using different sets of prior distributions. Scenarios with no ET are underlined. and the most noticeable differences across scenarios are in bold. Results are based on ${10}^4$ simulations.

	$Scenario$
	1	2	3	4	5	6	7	8	9
Proportion of correct selections
EmaxR	45.1	94.4	44.3	82.1	57.1	97.2	64.3	29.5	21.6
EmaxR ($\lambda$)	43.8	93.3	41.3	85.5	58.7	99.7	64.1	27.3	22.7
EmaxR ($E{D}_{50}$)	35.3	96.5	52.9	85.5	54.0	99.8	63.4	16.5	24.8
Proportion of toxicity responses
EmaxR	19.1	10.2	21.2	24.4	18.8	8.3	25.8	23.6	25.5
EmaxR ($\lambda$)	19.2	10.3	21.1	24.7	19.1	8.2	25.9	23.5	25.7
EmaxR ($E{D}_{50}$)	19.1	10.5	21.1	24.5	19.3	8.4	25.4	23.5	26.1
$\hat{\mathcal{P}}$
EmaxR	74.7	16.3	66.7	71.5	67.9	24.2	12.2	29.3	6.6
EmaxR ($\lambda$)	74.0	12.7	71.1	78.4	80.0	28.5	15.6	31.2	10.0
EmaxR ($E{D}_{50}$)	68.1	9.2	66.3	70.5	61.5	23.6	16.2	20.4	11.2
NMSE
EmaxR	1.7	2.0	2.2	3.5	3.1	1.4	2.5	1.7	3.7
EmaxR ($\lambda$)	1.9	1.8	2.2	3.3	3.7	1.4	2.5	1.7	3.6
EmaxR ($E{D}_{50}$)	2.0	2.2	2.2	4.3	3.2	1.9	3.0	2.3	3.1

In the majority of the scenarios, the choice of prior did not influence the PCS by more than 4% with scenarios 1, 3, and 8 being exceptions. In cases where the MTC lay in the middle of the curve (scenarios 1 and 8), the informative prior for $E{D}_{50}$ resulted in nearly 10% loss in the PCS. This set of prior distributions also results in slightly lower PCS in scenario 5 with the MTC being the in the middle as well. Conversely, in scenario 3 where the MTC is located at the beginning of the curve, the PCS was greater. At the same time, the proportion of toxicity responses are not influenced by the prior choice in all scenarios.

The proportion of times ET is found is again not affected by the choice of prior in scenarios 1–4, 6–7, and 9 with the difference between competing models below 10%. The major difference in $\hat{\mathcal{P}}$ can be found in scenario 5 in which the informative prior for $\lambda$ leads to 20% more false conclusions than the informative prior for $E{D}_{50}$. At the same time, it is generally harder for the EmaxR ($E{D}_{50}$) to detect the ET with the largest difference of 10% in scenario 8. Also, while NMSE is not largely affected by the choice of the prior, EmaxR ($E{D}_{50}$) results in slightly greater values in scenarios 1, 2, 4, 6–8 with the largest difference of 1 in scenario 4.

Overall, the choice of prior does not seem to have a noticeable effect in the majority of the cases, but EmaxR ($\lambda$) is more robust and results in a more accurate toxicity curve fitting. Therefore, an informative prior for $\lambda$ and an uninformative prior for $E{D}_{50}$ appear to be the preferable choice.

6. Randomization ratio

The randomization was shown to be the key factor in the dose-toxicity curve fit improving and in establishing the ET. We now consider the impact of different randomization ratios on the performance of the design. The prior given in (7) is used and the summary of the results is given in Table 5.

Table 5. Operating characteristics of EmaxR using different randomization ratios and no-skipping combination constraint. Scenarios with no ET are underlined and the most noticeable differences across scenarios are in bold. Results are based on ${10}^4$ simulations.

	$Scenario$
	1	2	3	4	5	6	7	8	9
Proportion of correct selections
EmaxR (2:1)	43.5	91.6	44.1	84.2	58.5	97.2	62.7	28.2	27.7
EmaxR (2:1) No Skip	46.9	90.7	44.8	83.7	57.2	96.7	62.3	27.1	29.2
EmaxR (3:1)	45.1	94.4	44.3	82.1	57.1	96.8	64.3	29.5	21.6
EmaxR (3:1) No Skip	45.7	90.4	47.1	83.8	55.9	97.2	66.5	31.0	21.5
EmaxR (4:1)	43.3	94.6	45.1	86.8	57.1	98.1	66.5	29.1	19.1
EmaxR (4:1) No Skip	44.0	88.4	48.2	89.2	58.8	98.0	69.2	29.7	19.5
Proportion of toxicity responses
EmaxR (2:1)	17.5	9.9	19.7	22.9	18.2	8.1	25.2	22.5	25.8
EmaxR (2:1) No Skip	17.6	9.5	19.9	22.5	17.8	8.5	24.8	22.4	25.5
EmaxR (3:1)	19.1	10.2	21.2	24.4	18.8	8.3	25.8	23.6	25.5
EmaxR (3:1) No Skip	18.8	9.9	21.5	24.3	18.9	7.9	25.6	23.1	25.5
EmaxR (4:1)	20.4	10.7	23.2	26.8	20.9	8.6	27.6	25.1	26.9
EmaxR (4:1) No Skip	19.5	10.3	23.1	26.0	20.7	8.9	26.8	24.4	26.8
$\hat{\mathcal{P}}$
EmaxR (2:1)	73.5	11.4	63.1	71.9	68.3	19.9	7.8	17.6	4.8
EmaxR (2:1) No Skip	69.1	12.1	63.8	69.6	69.1	20.0	8.5	16.7	4.6
EmaxR (3:1)	74.7	16.3	66.7	71.5	68.0	24.2	12.2	29.3	6.6
EmaxR (3:1) No Skip	75.6	12.2	68.7	72.8	71.3	24.1	11.0	27.2	7.7
EmaxR (4:1)	76.3	11.6	64.5	72.1	72.3	20.5	11.9	24.4	9.9
EmaxR (4:1) No Skip	74.5	14.4	69.2	69.9	66.1	20.9	13.1	23.1	9.9
NMSE
EmaxR (2:1)	1.7	2.1	3.4	3.3	3.2	1.7	2.4	1.6	3.6
EmaxR (2:1) No Skip	1.7	2.1	3.3	3.5	3.5	1.7	2.5	1.6	3.5
EmaxR (3:1)	1.7	2.0	2.2	3.5	3.1	1.5	2.5	1.7	3.7
EmaxR (3:1) No Skip	1.6	2.0	3.4	3.3	3.2	1.7	2.4	1.6	3.7
EmaxR (4:1)	1.6	1.9	3.2	3.4	3.3	1.5	2.5	1.7	3.7
EmaxR (4:1) No Skip	1.7	2.0	3.1	3.6	3.4	1.6	2.6	1.6	3.6

Comparing models with the same randomization ratio, the no skipping constraint has a minor influence on all characteristics of EmaxR (2:1), EmaxR (3:1), and EmaxR (4:1) in all scenarios. In the majority of cases, the difference does not exceed 4%. The only exception in terms of the PCS is EmaxR (4:1) in scenario 2 with the MTC being the last combination. Given a large cohort size ($c=5$), the total sample size $N=48$ is not always enough to test all combination and reach the MTC, which results in less accurate MTC selections by 6%. For EmaxR (4:1) no skipping constraint also increases $\hat{\mathcal{P}}$ by 5% in scenario 3 with the MTC being at the beginning of the curve and decreases it by 6% in scenario 5 with the MTC being the medium combination.

Considering different randomization ratios with combinations skipping, all approaches result in comparable PCS in all scenarios with EmaxR(2:1) in scenario 9 being the exception. The increased number of patients on the control arm (16 patients) gives a higher PCS in this scenario when the MTC is the SoC itself. At the same time, a smaller number of patients on the investigational arm might result in a substantial loss in the ability to detect the ET, for example, in scenario 8, $\hat{\mathcal{P}}$ decreased by 10% compared to EmaxR (3:1). Note that EmaxR(4:1) also results in a less accurate ET detection in this case as there are fewer patients on the control. Moreover, ratio 4:1 results in more toxic responses as more patients are assigned to the control. There is no noticeable difference in the NMSE across difference ratios.

Overall, while the operating characteristics remain unchanged by the randomization ratio in the majority of scenarios, a ratio of 3:1 seems to be a reasonable trade-off between the MTC selection and a plateau detection.

6. Discussion

The Bayesian model-based dose-escalation methods incorporating randomization to the control arm into the Emax model were studied and compared with widely used model-based alternatives with and without the control group. It was found that the inclusion of the randomization to the control group and the Emax model allow for the identification of the ET between the toxicity risk at the optimal dose combination and the single agent more frequently than alternative approaches. The cost of randomization to the control arm is presented by a small reduction in the proportion of correct selections in scenarios with the MTC being in the middle of a combination-toxicity curve. The sensitivity analysis suggests that the design is robust to the prior for reasonable choices and to the allocation ratio between investigation and control arms.

It is important to mention that the question of using the control group in Phase I trials is widely debated in the literature (e.g. Saad et al., 2017). The main argument against is the monotonicity assumption that makes comparing toxicity levels unnecessary and results in considering the control group as a “waste of resources” (not contributing to the MTC selection). However, this assumption has been found to be inappropriate for many immunotherapies. Moreover, it is found that the randomization between MTC and the control allows establishing that an effect observed for a combination is associated with the interaction of treatments and not only due to one of the therapies in the combination. Another advantage of the inclusion of the control is that it allows assessment for preliminary signs of activity which will help to make decisions about taking a treatment forward to the next phase of the development. Therefore, we argue that the control arm should be included in drug-combination trials. Importantly, the control group (patients allocated to the control arm) receives the SoC which makes the randomization an ethically viable option.

The proposed design uses a fixed randomization ratio during the trial which is specified prior to the trial. However, one can also consider adaptive randomization schemes which would change the randomization ratio based on the values of $\mathbb{P}\left(p_{j^{\star }}-p_0\ge \tau |{\mathbb{Y}}_n\right)$ estimated while the trial progress. For example, if $\mathbb{P}\left(p_{j^{\star }}-p_0\ge \tau |{\mathbb{Y}}_n\right)$ is large, then one has a high confidence that the ET is present. In the setting when a clinician does not tolerate an increase in toxicity risk by more than $\tau$, this might be a sign to change the randomization ratio and make it more balanced randomization for overdose control. Similarly, if $\mathbb{P}\left(p_{j^{\star }}-p_0\ge \tau |{\mathbb{Y}}_n\right)$ is small, more patients can be randomized to the MTC.

Evidently, the proposed method has the same limitations as the classical Bayesian CRM design. Among others, a large variance of a number of patients on the MTC as indicated by Oron and Hoff (2013). Therefore, it is hard to guarantee the specific level of power in the frequentist hypothesis that could also be of interest in the considered framework. Consequently, the design that could maximize the power of such test might be also of a great interest. Moreover, further investigation is needed to explore the conditions when randomization is able to overcome the fitting problem of the four-parameter model with a small number of patients which might be a major issue of the proposed approach. Finally, while Phase I clinical trials are considered in this paper only, a proper comparison of the optimal combination and control effectivenesses in Phase II trials is also essential (Sharma and Allison, 2015). In this case, a proper way to conduct this comparison should be investigated.

Acknowledgments

The authors acknowledge the insightful and constructive comments made by two reviewers. This project has received funding from the European Union‘s Horizon 2020 research and innovation program under the Marie Sklodowska-Curie grant agreement No 633567. Xavier Paoletti is partially funded by the Institut National du Cancer (French NCI) grant SHS-2015 Optidose immuno project. This report is independent research arising in part from Prof Jaki’s Senior Research Fellowship (NIHR-SRF-2015-08-001) supported by the National Institute for Health Research. The views expressed in this publication are those of the authors and not necessarily those of the NHS, the National Institute for Health Research or the Department of Health.

Additional information

Funding

This work was supported by the H2020 Marie Skłodowska-Curie Actions [Marie Sklodowska-Curie grant agreement No 633567]; Institut National Du Cancer [Grant SHS-2015 Optidose immuno project]; National Institute for Health Research [Prof Jaki’s Senior Research Fellowship NIHR-SRF-2].