Evaluation of Phase I Clinical Trial Designs for Combinational Agents along with Guidance based on Simulation Studies

Here we define some common notations used in these methods. As most of the designs we selected apply to dual-agent dose finding only (Copula has been extended to handle more than two agents), we assume two agents $A$ and $B$, with $J$ and $K$ doses respectively. Define $\pi_{jk}$ to be the true toxicity probability of the dose combination $(j,k)$, $j=1,2,\dots,J$, $k=1,2,\dots,K$; define $p_{j}$ to be the true toxicity probability of agent $A$ when used as a monotherapy, $j=1,2,\dots,J$, and $q_{k}$ to be the true toxicity probability of agent $B$ when used as a monotherapy, $k=1,2,\dots,K$. Define $\phi$ to be the pre-specified target toxicity probability. Define $N$ to be maximum sample size in the trial. Define $n_{jk}$ to be number of subjects that received dose $(j,k)$ and $y_{jk}$ to be number of dose-limiting toxicities (DLTs) observed among those $n_{jk}$ patients.

This method is a Bayesian design that extends CRM to accommodate dose-finding in dual agents.

Dose-toxicity model:

$$\pi_{jk}(\mathbf{\theta})=1-(a_{j})^{\alpha}(1-b_{k})^{\beta+\gamma\log(1-a_{j})},$$

(1)

where $\mathbf{\theta}=(\alpha,\beta,\gamma)$ is a vector of unknown design parameters and restricts $\alpha>0,\beta>0,\gamma<0$ to satisfy the assumption of toxicity monotonicity; $0\leq a_{1}<\dots<a_{J}$ and $0\leq b_{1}<\dots<b_{K}$ are constants instead of actual doses of agents. If no interaction between two agents exists in Equation 1, the model becomes

$$\pi_{jk}(\mathbf{\theta})=1-(a_{j})^{\alpha}(1-b_{k})^{\beta}.$$

(2)

Start-up phase: the trial is initiated with a dose combination $(1,1)$. Next, escalate agent $A$ while agent $B$ is maintained at the lowest dose if no DLT is observed. If still no DLT is observed when agent $A$ reaches its maximum dose, agent $B$ is escalated to its second lowest dose combining with agent $A$’s $(J-2)^{th}$ dose. Then if no DLT is observed, we continue to a combination where agent $B$ is at its third lowest dose and agent $A$’s $(J-4)^{th}$ dose, namely the one used in previous combination minus 2. When agent $B$ reaches its maximum dose, if agent $A$ is at its $m^{th}$ dose, we evaluate all combinations from $(m,K),(m+1,K),\dots,(J,K)$. The start-up phase ends if at least one DLT is observed at any time.

Post start-up escalation / de-escalation dose set: if the current combination is $(j,k)$, I2D only considers doses $(j-1,k),(j+1,k),(j,k),(j,k+1),(j,k-1),(j+1,k-1),(j-1,k+1)$ to the next subject prohibiting diagonal moves.

Post start-up trial conduct: after the start-up phase ends, the working model Equation 1 or Equation 2 will be used to obtain toxicity estimates of all dose combinations. Due to safety concerns, the working model starts at the combination dose where agent $B$ is at its lowest dose and agent $A$ is at the dose that makes the combination’s estimated toxicity probability closest to $\phi$.

MTD determination: the dose combination whose posterior probability of toxicity is closest to $\phi$ will be selected as the MTD.

This method utilizes copula to model the dose-toxicity relationship because copula allows to link the joint distribution and marginal distributions via a dependence parameter.

Dose-toxicity model:

$$\pi_{jk}=1-\{(1-p_{j}^{\alpha})^{-\gamma}+(1-q_{k}^{\beta})^{-\gamma}-1\}^{-1/\gamma},$$

(3)

where $\alpha$ and $\beta$ are power parameters as in CRM to accommodate the uncertainty, and $\gamma>0$ represents the interaction between two agents. Intermediate informative prior distributions with prior mean 1 and a relatively small variance will be assigned to $\alpha$ and $\beta$ (e.g. Gamma (2,2)). A Gamma distribution with a large variance is usually chosen as the non-informative prior for $\gamma$. If only one agent is involved, this approach reduces to regular CRM.

Start-up phase: the start-up phase begins with the lowest dose combination $(1,1)$. It proceeds vertically $(1,2),(1,3),\dots,(1,K)$ until the first toxicity is observed, then it proceeds horizontally $(2,1),(3,1),\dots,(J,1)$ until the first toxicity is observed. Once one toxicity is observed in both directions, the formal design starts.

Post start-up escalation / de-escalation dose set: if the current combination is $(j,k)$, the dose escalation set is defined as ${\cal A}_{E}=\{(j+1,k),(j,k+1),(j+1,k-1),(j-1,k+1)\}$, dose de-escalation set is define to be ${\cal A}_{D}=\{(j-1,k),(j,k-1),(j+1,k-1),(j-1,k+1)\}$. As we only know the partial order of dose toxicity in combinational agents, we do not know whether dose combinations $(j+1,k-1)$ and $(j-1,k+1)$ are more/less toxic than dose combination $(j,k)$. Therefore, the authors included $(j+1,k-1)$ and $(j-1,k+1)$ in both escalation and de-escalation sets.

Post start-up trial conduct: this design involves two parameters $c_{e}$ and $c_{d}$ that represent the fixed probability cut-offs for dose escalation and de-escalation, respectively, and $c_{e}+c_{d}>1$. Detailed algorithm is laid out as below.

1. a

   if at current dose combination $(j,k)$, $P(\pi_{jk}<\phi)>c_{e}$, we will escalate to the dose that belongs to ${\cal A}_{E}$ and with the toxicity probability closest to $\phi$ and higher than that of $(j,k)$. If current dose is $(J,K)$, then stay at the same combination.

2. b

   If at current dose combination $(j,k)$, $P(\pi_{jk}>\phi)>c_{d}$, we will de-escalate to the dose that belongs to ${\cal A}_{D}$ and with the toxicity probability closest to $\phi$ and lower than that of $(j,k)$. If current dose is $(1,1)$, the trial is terminated.

3. c

   Otherwise, stays at the same dose combination.

MTD determination: after $N$ subjects are exhausted, the MTD is determined as the dose combination with the estimated probability of toxicity closest to $\phi$.

To solve the problem of only knowing partial order of dose toxicity, POCRM proposes to pre-specify a subset of possible orderings, then utilize the CRM on each of them. This way, two-dimensional dose-finding is reduced to a one-dimensional problem.

Dose-toxicity model: define $T$ as the total number of dose combinations, $T=J\times K$; $d_{n}$ as the dose assigned to subject $n$; $R(d_{n})$ as the true toxicity probability of $d_{n}$; $y_{n}$ as a binary indicator of whether subject $n$ has toxicity or not, $n=1,2,\dots,N$; and $\Omega_{n}$ as data collected after having $n$ subjects where $\Omega_{n}=\{d_{1},y_{1},\dots,d_{n},y_{n}\}$. Assume we have $M$ possible partial ordering in total. For a specific ordering $m$, $m=1,2,\dots,M$, $R(d_{n})$ is modeled as below, similar to the CRM

$$R(d_{n})=E(Y_{n}|d_{n})\doteq\psi_{m}(d_{n},a)$$

(4)

where $\psi_{m}$ is some working dose-toxicity model, $a$ is the model parameter. After having $n$ patients, the likelihood under partial order $m$ is

$$L_{m}(a|\Omega_{n})=\prod_{i=1}^{n}\psi_{m}^{y_{i}}(d_{i},a)\{1-\psi_{m}(d_{i},a)\}^{(1-y_{i})}.$$

(5)

Then, the estimate of parameter $a$ under ordering $m$, $\hat{a_{m}}$, could be obtained through maximizing Equation 5. Then we could obtain the posterior probability of partial order $m$:

$$p(m|\Omega_{n})=\frac{p(m)L_{m}(\hat{a_{m}}|\Omega_{n})}{\sum_{l=1}^{M}p(l)L_{l}(\hat{a_{l}}|\Omega_{n})}.$$

(6)

Start-up phase: when POCRM was first proposed, it was a single stage method (Wages et al., 2011a, ). Later it was extended to include a start-up phase (Wages et al., 2011b, ). The start-up phase partitions the dose combination matrix to different “zones” and starts the first cohort with zone 1, which is the lowest dose combination. If no DLT is observed, then it assigns next cohort to doses in zone 2. If there are multiple combinations in zone 2, it randomly selects one of them. If no DLT is observed, it continues to assign the next cohort to other combinations in the same zone. Moving to the next zone is only allowed when all the dose combinations have been explored in lower zones. The start-up phase ends when one DLT is observed. POCRM also allows users to specify their own scheme in the start-up phase as they see appropriate.

Post start-up escalation / de-escalation dose set: since POCRM pre-specifies a subset of possible orderings, the escalation and de-escalation dose are known within each ordering.

Post start-up trial conduct: after the start-up phase ends, the authors used weighted randomization to select the partial ordering with $p(m|\Omega_{n})$ from Equation 6 being the weight. After selecting the working partial ordering $m$, $\pi_{t}$ is estimated for all $t\in\{1,2,\dots,T\}$ through Equation 4 and assign the dose combination that minimizes $|\hat{\pi_{t}}-\phi|$ to the next subject. But for final MTD determination, the ordering with the maximum posterior probability will be chosen among all candidate orderings.

MTD determination: after $N$ subjects are exhausted, the MTD is determined as the dose combination with the estimated probability of toxicity closest to $\phi$, given the ordering with maximum posterior probability.

Dose-toxicity model: This method employs a hierarchical model:

$$\pi_{jk}\sim Beta(\alpha_{jk},\beta_{jk}),$$

(7)

$$log(\alpha_{jk})=\theta_{0}+\theta_{1}a_{j}+\theta_{2}b_{k},$$

(8)

$$log(\beta_{jk})=\phi_{0}+\phi_{1}a_{j}+\phi_{2}b_{k},$$

(9)

where $\mathbf{\theta}=\{\theta_{0},\theta_{1},\theta_{2}\}$ follows a multivariate normal distribution with mean $\mathbf{\mu}=\{\mu_{0},\mu_{1},\mu_{2}\}$ and variance covariance matrix $\sigma^{2}\cal{I}$ where $\cal{I}$ is a $3\times 3$ identity matrix; $\mathbf{\phi}=\{\phi_{0},\phi_{1},\phi_{2}\}$ follows a multivariate normal distribution with mean $\mathbf{\omega}=\{\omega_{0},\omega_{1},\omega_{2}\}$ and the same variance covariance matrix; $a_{j}$ and $b_{k}$ are “effective doses” instead of actual clinical values. This method omits the interaction effects between two agents.

The authors provided recommendations about selecting priors and methods to calculate “effective doses”. They used the fact that $\frac{K\tilde{\pi}_{11}}{K(1-\tilde{\pi}_{11})}=\frac{exp\{\mu_{0}\}}{exp\{\omega_{0}\}}$ to obtain the solutions for $\mu_{0}$ and $\omega_{0}$:

$$\mu_{0}=log(K\tilde{\pi}_{11}),\omega_{0}=log(K(1-\tilde{\pi}_{11})).$$

They suggested setting $\mu_{1}=\mu_{2}=\omega_{1}=\omega_{2}=2\sqrt{\sigma^{2}}$ and selecting $\sigma^{2}\in[5,10]$. They set $a_{1}=b_{1}=0$, then

$$a_{j}=(\mu_{1}+\omega_{1})^{-1}\log(\tilde{OR}_{j.}),$$

$$b_{k}=(\mu_{2}+\omega_{2})^{-1}\log(\tilde{OR}_{.k})$$

where

$$\tilde{OR}_{j.}=\exp\{\frac{\tilde{\pi}_{j1}/(1-\tilde{\pi}_{j1})}{\tilde{\pi}_{11}/(1-\tilde{\pi}_{11})}\},$$

$$\tilde{OR}_{.k}=\exp\{\frac{\tilde{\pi}_{1k}/(1-\tilde{\pi}_{1k})}{\tilde{\pi}_{11}/(1-\tilde{\pi}_{11})}\}.$$

Therefore, this design needs inputs of $\pi_{j1}$ and $\pi_{1k}$, $j=1,2,\dots,J,k=1,2,\dots,K$ from the clinicians.

Start-up phase: this method is a single stage design without a start-up phase.

Escalation / de-escalation dose set: if the current combination is $(j,k)$, acceptable dose combination set to the next subject $S$ is defined as $(j-1,k)$, $(j+1,k)$, $(j,k)$, $(j,k+1)$, $(j,k-1)$, $(j+1,k-1)$, $(j-1,k+1)$, $(j+1,k+1)$, $(j-1,k-1)$. As dose combinations $(j+1,k+1)$ and $(j-1,k-1)$ are included, this design allows simultaneous dose escalation or de-escalation of both agents.

Trial conduct:

1. a

   Compute a $95\%$ CI for overall toxicity rate from cumulative data of currently recruited subjects.

2. b

   If the lower bound of this CI is greater than $\phi$, terminate the trial.

3. c

   If the lower bound of this CI is less than or equal to $\phi$, use all previous information to obtain posterior mean of $\pi_{jk}$, $j=1,2,\dots,J,k=1,2,\dots,K$.

4. d

   Select a dose that belongs to set $S$ and is closest to $\phi$. Assign this dose combination to the next patient.

5. e

   Continue until all $N$ subjects are exhausted.

MTD determination: Use the outcomes and assignments of all $N$ subjects to derive posterior mean of $\pi_{jk}$, $j=1,2,\dots,J,k=1,2,\dots,K$. If the last subject received dose $(j^{\prime},k^{\prime})$, then the dose combination that is among set $S$ of $(j^{\prime},k^{\prime})$ and with estimated toxicity closest to $\phi$ will be selected as the MTD.

Dose-toxicity model: this method uses logistic regression to link doses and toxicities of the two agents:

$$\text{logit}(\pi_{jk})=\beta_{0}+\beta_{1}a_{j}+\beta_{2}b_{k}+\beta_{3}a_{j}b_{k},$$

(10)

where $a_{j}$ and $b_{k}$ are “effective doses” instead of actual clinical values, $\beta_{1}>0$, $\beta_{2}>0$, $\beta_{1}+\beta_{3}b_{k}>0$, $\beta_{2}+\beta_{3}a_{j}>0$ to ensure monotonicity. The authors define “effective doses” as $a_{j}=\log(\frac{p_{j}}{1-p_{j}})$, $b_{k}=\log(\frac{q_{k}}{1-q_{k}})$ and recommend a vague normal prior $N(0,1)$ for $\beta_{0}$ and $\beta_{3}$, an informative prior exp$\{1\}$ for $\beta_{1}$ and $\beta_{2}$.

Start-up phase: the start-up phase starts from dose $(1,1)$. If no toxicity is observed, escalate the dose along the diagonal until at least one agent reaches maximum dose. If still no toxicity is observed when one agent reaches its maximum dose, we increase the dose of the other agent until both agents reach maximum doses. The start-up phase ends once the first toxicity is observed and the model-based design starts.

Post start-up escalation / de-escalation dose set: if the current combination is $(j,k)$, the dose escalation set is defined as ${\cal A}_{E}=\{(j+1,k),(j,k+1),(j+1,k-1),(j-1,k+1)\}$, dose de-escalation set is define as ${\cal A}_{D}=\{(j-1,k),(j,k-1),(j+1,k-1),(j-1,k+1)\}$. As we only know the partial order of dose toxicity in combinational agents, we do not know whether dose combinations $(j+1,k-1)$ and $(j-1,k+1)$ are more/less toxic than dose combination $(j,k)$. Therefore, the authors included $(j+1,k-1)$ and $(j-1,k+1)$ in both escalation and de-escalation sets.

Post start-up trial conduct: in the model-based design part, the escalation and de-escalation rule is the same as in Copula design (Yin and Yuan, 2009a, ).

MTD determination: DFCOMB utilizes a different method to identify MTD after the trial is completed. The dose combination that has the largest posterior probability $P(\pi_{jk}\in[\phi-\delta,\phi+\delta])$ and is used to treat at least one cohort will be selected as the MTD. Parameter $\delta$ is the length around the target toxicity probability.

This method is another generalization of the CRM.

Dose-toxicity model: it uses proportional odds logistic regression to model the dose-toxicity relationship:

$$\text{logit}(\pi_{jk})=\alpha_{k}+\beta a_{j},$$

(11)

where $\alpha_{k}$ is agent $B$ specific intercept, $k=1,2,\dots,K$; $\beta$ is a common coefficient across models; $a_{j}$ is the “effective dose” of agent $A$, $j=1,2,\dots,J$. For example, if agent $B$ has three dose levels, then gCRM will need three models: $\text{logit}(\pi_{j1})=\alpha_{1}+\beta a_{j}$, $\text{logit}(\pi_{j2})=\alpha_{2}+\beta a_{j}$, and $\text{logit}(\pi_{j3})=\alpha_{3}+\beta a_{j}$. Later these “sub” models will be aggregated together through a joint prior distribution that forces correlation among $(\alpha_{1},\alpha_{2},\dots,\alpha_{k})$. As observed, gCRM assumes no interaction between two agents as well.

In terms of parameters $\alpha_{k}$ and $\beta$, their paper assumes that $\beta$ follows a Gamma distribution with mean $\mu_{\beta}$ and variance $\sigma^{2}_{\beta}$, $\alpha_{1}\sim N(\mu_{\alpha},\sigma^{2}_{\alpha})$, and defines $\Delta_{k}=\alpha_{k}-\alpha_{k-1}\sim N(\delta_{k},2\sigma^{2}_{\alpha})$ for $k=2,3,\dots,K$ so the joint distribution of $\mathbf{\alpha}=(\alpha_{1},\dots,\alpha_{K})^{T}$ is multivariate normal. If one assumes that $\text{logit}(\pi_{j1})=E(\alpha_{1})+E(\beta)a_{j}$, $a_{j}\approx[\text{logit}(\pi_{j1})-\mu_{\alpha}]/\mu_{\beta}$ can be obtained. One can approximately obtain $\delta_{k}=\text{logit}(\pi_{1k})-\text{logit}(\pi_{1,k-1})$. The authors recommend setting $\mu_{\alpha}=-8$, $\mu_{\beta}=1$, $\sigma^{2}_{\alpha}=\sigma^{2}_{\beta}=1$. Therefore, with the inputs of $\pi_{j1}$ and $\pi_{1k}$ for $j=1,2,\dots,J,k=1,2,\dots,K$ from clinicians, all parameters can be calculated.

Start-up phase: this method is a single stage design without a start-up phase.

Escalation / de-escalation dose set: if the current combination is $(j,k)$, define acceptable dose combination set to the next subject $S$ to be $(j-1,k)$, $(j+1,k)$, $(j,k)$, $(j,k+1)$, $(j,k-1)$, $(j+1,k-1)$, $(j-1,k+1)$, $(j+1,k+1)$, $(j-1,k-1)$. As dose combination $(j+1,k+1)$ and $(j-1,k-1)$ are included, this design allows simultaneous dose escalation or de-escalation of both agents.

Trial conduct:

1. a

   Treat the first patient at dose $(1,1)$.

2. b

   For patient $2,3,\dots,N$, compute $\hat{\pi}_{jk}$ from logit$(\hat{\pi}_{jk})=\hat{\alpha}_{k}+\hat{\beta}a_{j}$ where $\hat{\alpha}_{k}$ and $\hat{\beta}$ are posterior means.

3. c

   As the posterior distribution of $\pi_{11}$ will be updated constantly, check that whether the stopping rule of $P(\pi_{11}>\phi)>0.95$ has been reached. If yes, then terminate the trial; otherwise assign next subject to the dose in set $S$ where $\hat{\pi}_{jk}$ is closest to $\phi$.

4. d

   Continue until all $N$ subjects are exhausted.

MTD determination: if the last subject received dose $(j^{\prime},k^{\prime})$, then the dose combination that is among set $S$ of $(j^{\prime},k^{\prime})$ and with estimated toxicity closest to $\phi$ will be selected as the MTD.

Bootstrap aggregating CRM is similar to POCRM as they both use one-dimensional CRM to identify the MTD. However, bCRM keeps updating the toxicity ordering of dose combinations rather than pre-specifying them.

Dose-toxicity model: in bCRM, it assigns a beta prior to $\pi_{jk}$ and obtain its posterior mean $\bar{\pi}_{jk}$. Then bCRM applies two-dimensional pool-adjacent-violators algorithm (PAVA) (Bril et al., , 1984) on $\bar{\pi}_{jk}$ to obtain $\tilde{\pi}_{jk}$ to ensure that these estimates satisfy partial ordering. To avoid ties among $\tilde{\pi}_{jk}$, a term $r_{jk}\epsilon$ is added, where $r_{jk}$ is the rank of dose $(j,k)$ and $\epsilon$ is a small positive number. The resulted estimates are denoted as $\tilde{\pi}^{\dagger}_{jk}$ and, as a result, one can obtain a new ordering $\cal{O}$. As noted by the authors, such orderings could vary dramatically due to data sparsity. Therefore, they bootstrapped $B$ samples of data to obtain corresponding orderings ${\cal O}_{b}$, and toxicity probability estimates $\hat{\pi}^{\text{b}}_{jk}$, $b\in 1,2,\dots,B$. The final estimate of $\pi_{jk}$ is

$$\hat{\pi}^{\text{Bagging}}_{jk}=\sum^{B}_{b=1}P({\cal O}_{b}|D)\hat{\pi}^{\text{b}}_{jk},$$

(12)

where D = $\begin{bmatrix}t_{1}&\dots&t_{n}\\ d_{1}&\dots&d_{n}\end{bmatrix}$ represents cumulative data up to $n^{th}$ subject, $t_{i}$ indicates whether subject $i$ experienced DLT or not; $d_{i}$ indicates dose combination subject $i$ received, $i=1,2,\dots,n$.

Start-up phase: the start-up phase is similar to DFCOMB.

Post start-up escalation / de-escalation dose set: similar to POCRM, bCRM uses one-dimensional CRM in the dose-finding process, therefore, the escalation and de-escalation dose is certain within each ordering.

Post start-up trial conduct: the trial conduct procedures are similar to DFCOMB.

MTD determination: After the trial is completed, one could select the combination that has been administered to patients and has the largest posterior probability of falling into the $\varepsilon$-neighbourhood of $\phi$, where $\varepsilon$ is a small positive number.

Combinational BOIN is a model-assisted design that is generalized from the single agent BOIN design (Liu and Yuan, , 2015; Yuan et al., , 2016).

Dose escalation and de-escalation rule: BOIN mainly involves two important parameters $\Delta_{L}$ and $\Delta_{U}$ which are lower and upper cut-offs. At current dose $j$, the escalation and de-escalation rules are below:

• •

  if $\hat{p_{j}}\in(\phi-\Delta_{L},\phi+\Delta_{U})$, then the next cohort stays at current dose;

• •

  if $\hat{p_{j}}\leq\phi-\Delta_{L}$, then the next cohort escalates to dose $j+1$;

• •

  if $\hat{p_{j}}\geq\phi+\Delta_{U}$, then the next cohort de-escalates to dose $j-1$;

where $\hat{p_{j}}$ is the estimated toxicity probability of dose $j$ in single agent dose-finding and it is simply proportion of patients experiencing toxicities among those who receive dose $j$.

In two-dimensional dose-finding, $\hat{p}_{jk}$ is calculated the same way: $\hat{p}_{jk}=y_{jk}/n_{jk}$.

An important task of cBOIN is to determine $\Delta_{L}$ and $\Delta_{U}$. Through minimizing the probability of incorrect movement given data at current dose,

$$\Delta_{L}=\phi-\frac{\log\{\frac{1-\phi_{1}}{1-\phi}\}}{\log\{\frac{\phi(1-\phi_{1})}{\phi_{1}(1-\phi)}\}},\Delta_{U}=\frac{\log\{\frac{1-\phi}{1-\phi_{2}}\}}{\log\{\frac{\phi_{2}(1-\phi)}{\phi(1-\phi_{2})}\}}-\phi.$$

The authors suggested using $\phi_{1}=0.6\phi$ and $\phi_{2}=1.4\phi$ through their simulation calibration.

Start-up phase: this method does not have a start-up phase.

Escalation / de-escalation dose set: admissible dose escalation set is defined as ${\cal A}_{E}=\{(j+1,k),(j,k+1)\}$, admissible de-escalation set is define as ${\cal A}_{D}=\{(j-1,k),(j,k-1)\}$.

Trial conduct:

1. a

   Treat the first cohort at dose $(1,1)$.

2. b

   Suppose that current dose is combination $(j,k)$. If $\hat{p}_{jk}\leq\phi-\Delta_{L}$, escalate to the dose combination that belongs to ${\cal A}_{E}$ and has the largest $P[p_{j\prime k\prime}\in(\phi-\Delta_{L},\phi+\Delta_{U})|y_{j\prime k\prime}]$.

3. c

   If $\hat{p}_{jk}\geq\phi+\Delta_{U}$, de-escalate to the dose combination that belongs to ${\cal A}_{D}$ and has the largest $P[p_{j\prime k\prime}\in(\phi-\Delta_{L},\phi+\Delta_{U})|y_{j\prime k\prime}]$.

4. d

   Otherwise if $\phi-\Delta_{L}<\hat{p}_{jk}<\phi+\Delta_{U}$, stay at current dose.

5. e

   Dose combinations with $P(p_{jk}>\phi|y_{jk})\geq\lambda$ will be permanently excluded, where $\lambda$ is pre-specified threshold probability. If dose combination $(1,1)$ satisfies this stopping rule, the trial will be terminated early.

6. f

   Continue until all $N$ subjects are exhausted.

MTD determination: after the trial is completed, isotonic regression will be used on $\hat{p}_{jk}$ to obtain estimator $\tilde{p}_{jk}$ so that they satisfy monotonic dose-toxicity within one agent when fixing the other agent’s dose. The MTD is the dose combination with $\tilde{p}_{jk}$ closest to $\phi$.

Similar to combinational BOIN, combinational Keyboard is a model-assisted design as well. Combinational Keyboard design starts with specifying a target toxicity interval ${\cal J}_{target}=(\phi-\varepsilon_{1},\phi+\varepsilon_{2})$, where $\varepsilon_{1}$ and $\varepsilon_{2}$ are tolerable deviations from $\phi$. This interval ${\cal J}_{target}$ is called target key. Then a series of equal-width keys are identified along both sides of the target key.

Dose escalation and de-escalation rule: in the setting of the single agent design, the escalation and de-escalation rules are straightforward. Define ${\cal J}_{\text{max}}$ to be the strongest key based on the posterior distribution of current dose $j$,

• •

  if ${\cal J}_{\text{max}}\prec{\cal J}_{target}$, then next cohort escalates to dose $j+1$;

• •

  if ${\cal J}_{\text{max}}\equiv{\cal J}_{target}$, then next cohort stays at current dose;

• •

  if ${\cal J}_{\text{max}}\succ{\cal J}_{target}$, then next cohort de-escalates to dose $j-1$;

To address two-dimensional dose-finding, the authors define five strategies of admissible escalate and de-escalate sets. After simulations, the strategy whose admissible escalation and de-escalation sets are the same with combinational BOIN design is recommended.

Start-up phase: this method does not have a start-up phase.

Escalation / de-escalation dose set: several dose assignment algorithms have been proposed for Keyboard design and the authors recommend to define admissible dose escalation set to be ${\cal A}_{E}=\{(j+1,k),(j,k+1)\}$, admissible de-escalation set to be ${\cal A}_{D}=\{(j-1,k),(j,k-1)\}$.

Trial conduct:

1. a

   Treat the first cohort at dose $(1,1)$.

2. b

   Suppose that current dose is combination $(j,k)$. If ${\cal J}_{\text{max}}\prec{\cal J}_{target}$, escalate to dose combination that belongs to ${\cal A}_{E}$ and has the largest $P[p_{j\prime k\prime}\in{\cal J}_{target}|(n_{jk},y_{jk)}]$.

3. c

   If ${\cal J}_{\text{max}}\succ{\cal J}_{target}$, de-escalate to dose combination that belongs to ${\cal A}_{D}$ and has the largest $P[p_{j\prime k\prime}\in{\cal J}_{target}|(n_{jk},y_{jk})]$.

4. d

   Otherwise if ${\cal J}_{\text{max}}\equiv{\cal J}_{target}$, stay at current dose.

5. e

   Dose combinations with $P(p_{jk}>\phi|y_{jk})\geq\lambda$ will be permanently excluded, where $\lambda$ is pre-specified threshold probability. If dose combination $(1,1)$ satisfies this stopping rule, the trial will be terminated early.

6. f

   Continue until all $N$ subjects are exhausted.

MTD determination: after the trial is completed, isotonic regression will be used to identify the MTD.

A total of 15 scenarios are displayed in Table 2. In the first 10 scenarios, agent $A$ has 5 dose levels and agent $B$ has 3. In scenarios 11 to 15, both agents have 4 dose levels. Target toxicity rate 0.3 is bolded. Among the first ten $5\times 3$ matrices: scenario 1 contains multiple MTD locations that are in the middle of matrix and diagonally connected; scenarios 2 and 4 represent over-toxic situations while scenario 4 is more extreme; scenarios 3 and 5 represent over-conservative situations while scenario 5 is more extreme; scenarios 6 and 7 contain multiple MTD locations but those locations are more scattered; scenarios 8, 9, and 10 contain single MTD at different locations. Among the last five $4\times 4$ square matrices: scenario 11 contains multiple MTD locations that are in the middle of matrix and diagonally connected; scenarios 12 and 13 contain multiple but more scattered MTD locations; scenario 14 contains two MTD locations that are at the bottom left and top right; scenario 15 has single MTD location.

Four evaluation metrics are used: (1) correct MTD selection $S_{C}$, defined as proportion of simulation runs that correctly identified the MTD among all 2000 simulations; (2) over-toxic MTD selection $S_{OT}$, defined as the proportion of simulation runs that identified over-toxic doses as MTD among all 2000 simulations; (3) correct patient assignment $A_{C}$, defined as the average proportion of patients that were assigned to the MTD during the trial across all 2000 simulations; (4) over-toxic patient assignment $A_{OT}$, defined as the average proportion of patients that were assigned to over-toxic doses during the trial across all 2000 simulations. Metrics $S_{C}$ and $S_{OT}$ will be used to evaluate the performance of designs in terms of MTD selection. The larger the $S_{C}$ is, the more accurate the design is in selecting the correct MTD. The larger the $S_{OT}$ is, the more aggressive the design is in selecting MTD. Metrics $A_{C}$ and $A_{OT}$ will be used to evaluate the characteristics of designs during trial conduct. The larger the $A_{C}$ is, the more accurate patient assignment is during the trial. The larger the $A_{OT}$ is, the more aggressive the design is in dose escalation during the trial. Ideally, a design should show relatively large $S_{C}$ and $A_{C}$ but small $S_{OT}$ and $A_{OT}$.

For I2D, we implemented published R codes (Ezzalfani, , 2019). We set cohort size of start-up phase to be 1 based on suggestions from simulation studies when target toxicity rate is 0.3 (Ivanova et al., , 2003; Wang and Ivanova, , 2005), and interaction to be 0 so that it is consistent with the paper’s focus. The prior of parameters $(\alpha,\beta)$ is the product of two independent exponential distributions with mean 1, which is the same as the one used in the I2D study (Wang and Ivanova, , 2005).

For Copula, the website (http://www.blackwellpublishing.com/rss) where simulation programs were originally published is not accessible now, so we used the executable file on the website https://odin.mdacc.tmc.edu/~yyuan/index_code.html. The executable file used escalation and de-escalation probability boundaries fixed at 0.8 and 0.45, respectively.

Hierarchy was implemented via R codes from website http://www-personal.umich.edu/~tombraun/software.html. We set $\sigma^{2}$ to be 10 based on the suggestion in the paper (Braun and Wang, , 2010). Together with other recommendations from the authors, we could obtain priors for all involved parameters. Details are described in Section  2.5.

POCRM was implemented via R package pocrm. We utilized six possible partial ordering as suggested (Wages and Varhegyi, , 2013; Wages and Conaway, , 2013): across rows, across columns, up diagonals, down diagonals, up-down diagonals, and down-up diagonals. As we do not have information about which partial ordering is more likely, the prior probabilities of all 6 partial ordering were set to be equal. The skeleton required by the program was obtained using getprior function in package dfcrm from algorithm of Lee and Cheung (Lee and Cheung, , 2009) as suggested (Wages and Varhegyi, , 2013). For the start-up phase, we used the “zoning” method as suggested (Wages et al., 2011a, ; Wages et al., 2011b, ).

Design cBOIN was implemented via R package BOIN. For cBOIN, the interval boundaries were set to be 0.18 and 0.42 as suggested (Lin and Yin, , 2017).

Design cKeyboard was implemented via R package Keyboard.

DFCOMB was implemented via R package dfcomb. As recommended by the authors, a vague normal prior $N(0,1)$ for $\beta_{0}$ and $\beta_{3}$, an informative prior exp$\{1\}$ for $\beta_{1}$ and $\beta_{2}$ (Riviere et al., , 2014). For DFCOMB, we set the target toxicity boundaries as 0.18 and 0.42 to be consistent with cBOIN. As one of our reviewers suggested, we also tried setting these toxicity boundaries as 0.25 and 0.35.

Design gCRM was implemented via R codes from website http://www-personal.umich.edu/~tombraun/software.html. As the authors suggested, we used a Gamma prior with mean 1 and variance 1 for $\beta$, normal prior with mean -8 and variance 1 for $\alpha_{1}$, normal prior with mean $\text{logit}(\pi_{1k})-\text{logit}(\pi_{1,k-1})$ and variance 2 for $\delta_{k}$, where $k=2,3,\dots,K$ (Braun and Jia, , 2013).

Design bCRM was implemented via R codes from the authors. Its skeleton setting is the same as POCRM.

For most model-based designs, there are several design parameters involved. Some of these parameters have recommended specifications provided by the authors, for example, interval boundaries in design cBOIN are recommended to set as $0.6\phi$ and $1.4\phi$ where $\phi$ is the target toxicity probability (Lin and Yin, , 2017). However, some design parameters lack authors’ suggested specifications and their influences to design performances are not clear. Therefore, we list such design parameters and corresponding designs in Table 1 where column “Main setting” and column “Alternative setting” contain parameter specifications used in our simulation studies.

For all designs and scenarios, the maximum sample size was set to be 60. This number is widely used in other studies (Yin and Yuan, 2009a, ; Braun and Jia, , 2013; Lin et al., , 2016; Lin and Yin, , 2017; Pan et al., , 2020; Riviere et al., 2015a, ). To verify that our conclusions are still valid with a different sample size, we repeated all simulations with a maximum sample size of 30.

Table 3, Table 4, Table 5, and Table 6 display design performances of $S_{C}$, $S_{OT}$, $A_{C}$, and $A_{OT}$, respectively. In these tables, I2D with parameter $p_{j}$ and $q_{k}$ specification in column “Main setting” of Table 1 is denoted as “I2D”; I2D with parameter $p_{j}$ and $q_{k}$ specification in column “Alternative setting” is denoted as “I2D.pq”. Copula with parameter $p_{j}$ and $q_{k}$ specification in column “Main setting” is denoted as “Copula”; Copula with parameter $p_{j}$ and $q_{k}$ specification in column “Alternative setting” is denoted as “Copula.pq”. Hierarchy with $\pi_{1k}$ and $\pi_{j1}$ specification in column “Main setting” is denoted as “Hierarchy”; Hierarchy with $\pi_{1k}$ and $\pi_{j1}$ specification in column “Alternative setting” is denoted as “Hierarchy.pi”. POCRM with skeleton specification in column “Main setting” is denoted as “POCRM”; POCRM with skeleton specification in column “Alternative setting” is denoted as “POCRM.skeleton”. DFCOMB with $p_{j}$ and $q_{k}$, and escalation/de-escalation probability cutoff specifications in column “Main setting” is denoted as “DFCOMB”; DFCOMB with parameter $p_{j}$ and $q_{k}$ specification in column “Alternative setting” is denoted as “DFCOMB.pq”; DFCOMB with escalation/de-escalation probability cutoff specification in column “Alternative setting” is denoted as “DFCOMB.cut”; DFCOMB with $p_{j}$ and $q_{k}$, and escalation/de-escalation probability cutoff specifications in column “Main setting”, but with toxicity boundaries being 0.25 and 0.35 as suggested by one of our reviewers, is denoted as “DFCOMB.sensitive”. Design gCRM with $\pi_{1k}$ and $\pi_{j1}$ specification in column “Main setting” is denoted as “gCRM”; gCRM with $\pi_{1k}$ and $\pi_{j1}$ specification in column “Alternative setting” is denoted as “gCRM.pi”. Design bCRM with skeleton and escalation/de-escalation probability cutoff specifications in column “Main setting” is denoted as “bCRM”; bCRM with skeleton specification in column “Alternative setting” is denoted as “bCRM.skeleton”; bCRM with escalation/de-escalation probability cutoff specification in column “Alternative setting” is denoted as “bCRM.cut”. We marked designs metrics that are “outstandingly” poor as red, and those that are poor, but not “outstanding” from the others as magenta.

I2D shows unstable performances in MTD identification in most simulation scenarios. While under extreme conditions like scenario 4 and 5, its $S_{C}$ is among the best, under scenarios like scenarios 1, 3, 8, 9, and 15, its $S_{C}$ is among the worst. With alternative toxicity profile of individual agents ($p_{j}$ and $q_{k}$), some scenarios had improved performances while some had worse, but its overall characteristics remain the same. Part of the reason of the unstable performance is that, I2D always starts its model-based part with agent $B$’s lowest dose no matter what happens in the start-up phase. This way, the starting dose combination of model-based part could be very far from the location of true MTDs, which makes I2D harder to identify them. Overall I2D is not an aggressive design as its $S_{OT}$ and $A_{OT}$ are relatively small under most scenarios.

Copula performed poorly in MTD identification as its $S_{C}$ and $S_{OT}$ are among the worst in several scenarios. This indicates that Copula is quite aggressive as it is more likely to select higher dose combinations as the MTD. With alternative toxicity profile of individual agents ($p_{j}$ and $q_{k}$), Copula’s performances fluctuated in several scenarios, but overall poor performances and aggressiveness were still observed. One potential reason of the unsatisfactory performances is the limited flexibility on changing parameters like escalation and de-escalation probability cutoffs in the executable file. Therefore, it is reasonable to argue that the default values (0.8 and 0.45) are not optimal for some scenarios. On the other hand, as we do not know what the true dose toxicity matrix looks like in real life, obtaining a uniform parameter set to achieve best performances under all scenarios through simulation calibration is not feasible.

Hierarchy is quite aggressive overall in trial conduct as it has the worst $A_{OT}$ under several scenarios. We can observe that incorrect $\pi_{j1}$ and $\pi_{1k}$ performed much worse than using correct ones in most scenarios. A possible reason for the aggressiveness is that, unlike most other designs, Hierarchy allows simultaneous dose escalation of both agents during trial conduct. Another aspect we should emphasize is that despite a high proportion of patients assigned to over-toxic doses, Hierarchy did not outperform other designs in terms of $S_{C}$. Some features of Hierarchy like omitting the interaction effect between agents and no start-up phase may contribute to this poor performance as well.

POCRM shows satisfactory characteristics across all scenarios except low $S_{C}$ in scenario 5. Then we found that the alternative skeleton setting in Table 1 improved $S_{C}$ from 0.54 to 0.73 in scenario 5. However, the alternative skeleton setting led to much worse performance metrics in scenario 15. Another finding is that POCRM performed well under scenarios (e.g., scenario 9) when the underlying true toxicity orderings of dose combinations are not among any of the six orderings we used. Such results “validate” the idea of POCRM: in practice we do not need to specify the correct toxicity ordering in POCRM, providing orderings close to the correct one is sufficient.

DFCOMB performed poorly in terms of $S_{C}$ under several scenarios. But overall DFCOMB is not an aggressive design as its $S_{OT}$ and $A_{OT}$ are relatively small under most scenarios. With alternative escalation and de-escalation probability cutoffs, some scenarios had improved $S_{C}$ (e.g. scenarios 9 and 10) while some had worse $S_{C}$ (e.g. scenarios 6 and 8). In addition, we observed worse $A_{OT}$ and better $S_{OT}$ with alternative specification (0.6 and 0.6 as escalation and de-escalation probability cutoffs). Worse $A_{OT}$ is expected as the alternative cutoff pairs make DFCOMB easier for dose escalation but more difficult for de-escalation. With more patients assigned to over-toxic doses, the estimation of dose toxicity probabilities were more accurate, which leads to better $S_{OT}$. With alternative toxicity profile of individual agents ($p_{j}$ and $q_{k}$), a few scenarios showed obvious impact: $S_{C}$ in scenario 1 was improved from 0.54 to 0.69, $S_{C}$ in scenario 11 was improved from 0.44 to 0.57, $S_{C}$ in scenario 14 was improved from 0.18 to 0.35, but $S_{C}$ in scenario 10 and 15 were dramatically reduced from 0.47 to 0.14, and from 0.67 to 0.45, respectively. With the “sensitivity run” of target toxicity boundaries suggested by one of our reviewers, we observed slightly worse $S_{OT}$ in some scenarios. Such results imply that the optimal design parameters are scenario-dependent. Therefore, it is not feasible for us to calibrate the parameters through simulations in real life clinical trials.

Design gCRM performed well in most scenarios. Comparing results using correct $\pi_{j1}$ and $\pi_{1k}$ with incorrect ones, we observed that using incorrect inputs leads to worse operating characteristics under some scenarios, and similar performances under the other ones. Interestingly, although gCRM allows simultaneous dose escalation of both agents, it did not show much aggressiveness as its $S_{OT}$ and $A_{OT}$ are not among the largest under most scenarios. One possible reason could be that gCRM uses single patient as unit, instead of cohorts during trial conduct. Therefore, every time when an over-toxic dose combination is assigned, only one patient instead of a cohort of several patients will receive it. From this perspective, gCRM could be viewed as more “flexible” in the dose-finding process and such flexibility may dilute the aggressiveness.

Design bCRM is another one whose performances were unstable across different scenarios. Its $S_{C}$ in scenario 5 and $S_{OT}$ in scenario 4 are among the worst. Metrics in other scenarios are acceptable. Similar to DFCOMB, we observed worse $A_{OT}$ with alternative escalation and de-escalation probability cutoffs. However, $S_{OT}$ was not improved. Alternative skeleton setting has influences to bCRM as well as it improved the $S_{C}$ in scenario 5 and $S_{OT}$ in scenario 4, but worsened $S_{C}$ in scenario 10 and $S_{OT}$ in scenario 9 and 10. Therefore, unstable performances remains an issue even using alternative skeleton setting or escalation and de-escalation probability cutoffs.

Both designs cBOIN and cKeyboard performed well across all scenarios. They may not always be the top performers, but their operating characteristics are never among the worst. This is especially important in real life clinical trials, as we do not know which scenario could be the truth. Therefore, cBOIN and cKeyboard are able to guarantee satisfactory performances in practice.

Results using 30 as the maximum sample size are shown in Supplementary Materials Table 5 to Table 12. All the findings discussed above are observed when maximum sample size is reduced from 60 to 30. Additionally, we found that in Supplementary Materials Table 5, POCRM with alternative skeleton had unsatisfactory $S_{C}$ in scenario 15 given maximum sample size of 30. This indicates that parameter calibration is not feasible for POCRM as well.