Modeling synergism in early phase cancer trials with drug combination with continuous dose levels: is there an added value?

In parametric Bayesian designs of early phase cancer clinical trials with drug combinations exploring a discrete set of partially ordered doses, several authors claimed that there is no added value in including an interaction term to model synergism between the two drugs. In this paper, we investigate these claims in the setting of continuous dose levels of the two agents. Parametric models will be used to describe the relationship between the doses of the two agents and the probability of dose limiting toxicity and efficacy. Trial design proceeds by treating cohorts of two patients simultaneously receiving different dose combinations and response adaptive randomization. We compare trial safety and efficiency of the estimated maximum tolerated dose (MTD) curve between models that include an interaction term with models without the synergism parameter with extensive simulations. Under a selected class of dose-toxicity models and dose escalation algorithm, we found that not including an interaction term in the model can compromise the safety of the trial and reduce the pointwise reliability of the estimated MTD curve.


Introduction
The primary goal of cancer phase I trials is to study the safety of new cytotoxic, biologic, or immunotherapy drugs or combination of existing agents with possibly radiation therapy and to determine the maximum tolerated dose (MTD) for use in future efficacy studies. These trials enroll advanced stage cancer patients who have exhausted conventional treatment (Roberts et al., 2004). For ethical considerations and to minimize the number of subjects exposed to toxic doses, patients are accrued to these trials sequentially and dose allocation is adaptive and depends on the doses given to current and previously treated patients and their toxicity outcomes. For treatments where only one drug is allowed to vary during the trial, and assuming a non-decreasing relationship between the dose levels and the probability of dose limiting toxicity (DLT), the MTD is defined as the maximum dose γ that will produce DLT in a prespecified proportion θ of patients (Gatsonis & Greenhouse, 1992) Typically, DLT is a binary indicator of toxicity and is defined as a grade 3 or 4 non-hemathologic or grade 4 hematologic toxicity, although models with quasi-continuous (Chen et al., 2010(Chen et al., , 2012 and ordinal toxicity grades (Diniz et (Clertant & O'Quigley, 2017) and nonparametric dose-finding methods, e.g., the modified toxicity probability interval (Ji & Wang, 2013), the Bayesian optimal design (Yuan et al., 2016), and the nonparametric overdose control method (Lin & Yin, 2017a).
To reduce treatment resistance to therapy, a combination of cytotoxic, biologic, and immunotherapy drugs is often used in order to target different signaling pathways simultaneously and hence achieve higher response rates by using additive or synergistic agents. When all these drugs are allowed to vary during the trial, a challenging problem is to identify a subset of these combinations among a larger set of permissible doses that will produce the same DLT rate θ since the MTD is no longer unique. Model-based designs for drug combinations that recommend a single combination for use in future efficacy studies have been proposed by Riviere et al. (2014), Shi & Yin (2013), Wages et al. (2011), and Yin & Yuan (2009a,b) among others. However, since different MTDs may have different levels of efficacy, we argue for the use of drug combination designs that recommend more than one MTD and study their efficacy in future randomized or response adaptive designs. To this end, and for discrete dose levels, Thall et al., 2003 proposed a six-parameter dose-toxicity model and a two-stage algorithm that recommends three MTDs at the conclusion of the trial. Other parameteric and nonparametric designs that can recommend more than one MTD include the three-parameter regression model of Wang & Ivanova, 2005, a Bayesian hierarchical model by Braun & Wang, 2010, the product of independent beta probabilities (PIPE) (Mander & Sweeting, 2015), the two-dimensional Bayesian optimal interval design (2dBOIN) (Lin & Yin, 2017b), and the fully nonparametric approach of Razaee et al., 2022. In cancer treatment using novel agents, the use of continuous dose levels is very common, particularly when the drugs are delivered as infusions intravenously. Tighiouart et al., 2014 introduced a dose-finding method for continuous dose levels based on conditional EWOC and recommended an MTD curve at the conclusion of the trial. The method was extended to three drugs (Tighiouart et al., 2016) In the above model-based designs for drug combinations, several authors argued for the unnecessary use of an interaction term in the model due to the inherent small sample size of these trials and inability to estimate the synergistic parameter of interest with good precision (Lyu et al., 2019;Riviere et al., 2014;Wang & Ivanova, 2005). In fact, Mozgunov et al., 2021 carried out simulation studies to assess the effect of omitting the interaction term on the percent of dose recommendation in the logistic model of Riviere et al., 2014. In general, they found that only marginal benefits are seen when using an interaction term and in the absence of prior information regarding the extent of synergism between the two drugs, there is no added benefit in including the interaction coefficient. In this manuscript, we investigate these claims in the setting of drug combinations with continuous dose levels using the design of Tighiouart, 2022; Tighiouart et al., 2017. We show using simulation studies that when the two drugs are highly synergistic, safety of the trial is compromised when omitting the interaction term from the model and efficiency of the estimated MTD curve as assessed by the pointwise bias and percent selection is greatly reduced.
The article is organized as follows. In Section 2, we describe the dose-toxicity model for drug combinations and review the dose-escalation algorithm based on the conditional EWOC scheme. Simulation studies comparing the performance of models with and without an interaction term are found in Section 3 and a discussion of these findings is included in Section 4.

Phase I Design
We consider the setting of two drugs A and B with continuous dose levels x and y in the intervals [X min , X max ] and [Y min , Y max ], respectively. We further standardize the doses to be in the interval [0, 1] using the transformations h 1 (x) = (x-X min )/(X max -X min ) and h 2 (y) = (y-Y min )/(Y max -Y min ).

Dose-Toxicity Model
We assume a class of dose-toxicity models for dose allocation and MTD curve estimation of the form where T is the indicator of DLT, T = 1 if a patient given the dose combination (x, y) exhibits DLT, and T = 0 otherwise, the interaction parameter β 3 ≥ 0, and F is a known cumulative distribution function (cdf ). Finally, we assume that that the probability of DLT is non-decreasing with the dose of any one of the agents when the other one is held constant. A necessary and sufficient condition for this property to hold is to assume β 1 , β 2 > 0. This is a common assumption for cytotoxic, biologic, immunotherapy drugs or radiation treatment. By definition, the MTD is any dose combination (x * , y * ) such that P(T = 1|x * , y * ) = θ. The target probability of DLT θ is pre-specified by the clinician and depends on nature and severity of expected DLTs.
A convenient reparameterization of model (2) in terms of parameters clinicians can interpret and to facilitate prior distribution elicitation is to use ρ 10 , the probability of DLT when the levels of drugs A and B are 1 and 0, respectively, ρ 01 , the probability of DLT when the levels of drugs A and B are 0 and 1, respectively, and ρ 00 , the probability of DLT when the levels of drugs A and B are both 0. Under model assumption (2), the MTD is the set of dose combinations Let π(ρ 00 , ρ 01 , ρ 10 , η 3 ) be a prior distribution of the model parameters and . . , k} be the data after enrolling k patients to the trial. The posterior distribution is π(ρ 00 , ρ 01 , ρ 10 ,
Here, Π Γ A|B=y is the posterior cumulative distribution function of the MTD of drug A given that the dose of drug B is y. This posterior distribution is easily obtained using the MCMC samples of the model parameters since Γ A|B=y and Γ B|A=x are functions of ρ 00 , ρ 01 , ρ 10 , and η 3 .
In step 2 of the above algorithm, a dose for drug A is selected in such a way that the posterior probability of exceeding the MTD of drug A given the current dose level of drug B is bounded by a feasibility bound α. For example, after resolving the DLT status of the first two patients, the third patient receives dose combination (x 3 , y 3 ) where y 3 = y 1 and x 3 is the α-th percentile of Π Γ A|B=y 1 , i.e., the largest dose such that the posterior probability of exceeding the MTD of A given that B = y 1 is no more than α. This is the defining property of EWOC.
Stopping rule. For ethical considerations, after n 1 patients are evaluable for toxicity, enrollment to the trial is suspended if there is statistical evidence that the minimum dose combination is too toxic, i.e., P(P(T = 1|(x, y) = (0, 0)) > θ + ξ 1 |D n 1 ) > ξ 2 . Here, ξ 1 , ξ 2 are design parameters chosen to achieve good operating characteristics under various scenarios for the location of the true MTD curve.
At the conclusion of the trial, an estimate of the MTD curve in (3) is given bŷ whereρ 00 ,ρ 01 ,ρ 10 ,β 3 are the posterior medians given D n .

Simulation Studies
In this section, we compare the performance of this design in estimating the MTD curve between the working model F(β 0 +β 1 x+β 2 y+β 3 xy) and the working model that does not use an interaction term F(β 0 + β 1 x + β 2 y), with F(·) being the logistic function F(u) = (1 + e -u ) -1 .

Scenarios and prior distributions
We consider four scenarios for the location of the true MTD curve shown by a solid black curve in Figure 1. The target probability of DLT is θ = 0.33. Scenario 1 is a case where the MTD of drug A when dose level of drug B is 0, Γ A|B=0 , and the MTD of drug B when dose level of drug A is 0, Γ B|A=0 , are within the range of doses available in the trial and are around 0.8. Therefore, the true values of ρ 01 and ρ 10 must be greater than θ and are taken to be 0.6. Convexity of the true MTD curve is controlled by the interaction parameter β 3 , taken to be 40. A small value for ρ 00 = 0.01 is selected. Scenario 2 is a case where Γ A|B=0 is within the range of doses of drug A but Γ B|A=0 is outside the range of doses of drug B. Therefore, we took ρ 10 to be higher than θ and ρ 01 to be smaller than θ.
The interaction parameter was selected to be much higher than the value in scenario 1. Scenario 4 is a case of a model misspecification away from the family of link functions F(β 0 + β 1 x + β 2 y + β 3 xy).
The true model parameters have been selected in such a way that Γ A|B=0 ≈ Γ B|A=0 ≈ 1.0 and the MTD curve has an inflection point. The true parameter values (ρ 00 , ρ 00 , ρ 00 , β 3 ) of the true logistic and six-parameter models are shown in the first column of Table 1.
For scenarios (1-3), estimates of the MTD curves are obtained using true and misspecified models for the dose-toxicity relationship. These are the logistic, probit F(u) = Φ(u), and complementary log-log F(u) = 1 -e -e u link functions, where Φ(·) is the cdf of the standard normal distribution. The misspecified models were chosen so that they have the same MTD curve as the logistic model by shifting the intercept coefficient. For example, when F(·) is the logistic function, points (x ⋆ , y ⋆ ) on the MTD curve satisfy F -1 (θ) = β 0 + β 1 x ⋆ + β 2 y ⋆ + β 3 x ⋆ y ⋆ . The corresponding MTD curve for the probit model is any dose combination ( For the two models to have the same MTD curve, we set β i = β ′ i , i = 1, 2, 3 and β ′ 0 = β 0 +F -1 (θ)-Φ -1 (θ). Figure 2(a,b) show the extent of departure of the true logistic model from the models with probit and complementary log-log link functions under scenario 3. In scenario 4, a six-parameter model is used as the true model P(T = 0|x, y) = (1 + α 1 x β 1 + α 2 y β 2 + α 3 (x β 1 y β 2 ) β 3 ) -1 . The true MTD curve and contour plots under this model are shown in Figure 2 considerations, dose escalation of either drug during the trial is restricted to be no more than 20% of the dose range of that drug.

Operating characteristics
Summary statistics for trial safety are reported as the percent of DLTs across all patients and all simulated trials and the percent of trials with DLT rates exceeding θ + 0.05 and θ + 0.1. We also report an estimate of the MTD curvē C = (x * , y * ) : y * = (F -1 (θ) -F -1 (ρ 00 )) -(F -1 (ρ 10 ) -F -1 (ρ 00 ))x * (F -1 (ρ 01 ) -F -1 (ρ 00 )) +β 3 where F(·) is the logistic function andρ 00 ,ρ 01 ,ρ 10   When the working model is the logistic link function without an interaction term, the estimatē C in (6) is a decreasing line. Efficiency of the estimate of the MTD will be assessed by the pointwise average bias and percent selection. For j = 1, . . . , m, let C j be the estimated MTD curve (based on n patients) and C true be the true MTD curve. For every dose combination (x, y) ∈ C true , consider the relative minimum distance of (x, y) to C j where y ′ is such that (x, y ′ ) ∈ C j . The sign of d (j) (x,y) depends on whether (x, y) is above or below the estimated MTD curve C j . The pointwise average bias is Analogous to the percent selection for discrete dose combinations, we define the pointwise percent selection for continuous dose levels as follows. For each dose combination (x, y) ∈ C true , let ∆(x, y) be the Euclidean distance from the minimum dose combination to (x, y). The percent selection is the fraction of trials C j , j = 1, . . . , m falling inside the circle with center (x, y) and radius p∆(x, y), 0 < p < 1, where p is a tolerance probability. This is a measure of percent of trials with MTD recommendation within (100 × p)% of the true MTD. Figure 2(d) shows the variability of the radius of tolerance circles for various combinations on the MTD curve under two scenarios when the tolerance probability is p = 0.1. Using the distance formula in (7), the pointwise percent selection with tolerance probability p can be expressed as where 1(·) is the indicator function. Table 1 shows that under scenarios 1 and 2, the average DLT rate across all simulated trials is close to the target θ = 0.33 under the true and misspecified models when the working logistic model uses an interaction term. However, when omitting the interaction coefficient for the logistic working model under the true logistic model, this average DLT rate increases to 0.36 and 0.37 under these scenarios. Furthermore, the percent of trials with DLT rate exceeding θ + 0.5 can be at least 15% higher when using a working model without an interaction term under scenario 1 and up to 20% higher under scenario 2. A similar trend is observed under scenarios 3 and 4 but the trials are still safe with respect to the average DLT rate and percent of trials with a DLT rate exceeding θ + 0.5. In all cases, the percent of trials with an excessive rate of DLT (greater than θ + 0.1) is negligible. These results show that safety of the trial is compromised when the underlying synergism between the two drug is high and a model without an interaction term is used for dose finding.  Figure 1 shows that the estimated MTD curves using the logistic, probit, complementary loglog working models are much closer to the true MTD curve relative to a logistic model without an interaction term, for the majority of dose combinations along the true MTD curve. This is also reflected in the pointwise average bias displayed in Figure 3. Similar to the findings in (Diniz et al., 2017;Jiménez et al., 2019;Tighiouart et al., 2017Tighiouart et al., , 2014, the pointwise average bias is less than 5% of the dose range except perhaps near the edges of the true MTD curve. However, omitting an interaction term leads to a substantial increase in this average bias by as much as 20% of the dose range under scenario 1 for standardized doses of drug A around 0.1 or 0.7 and around 0.05 and 0.5 under scenario 2. Under scenario 3, when using the logistic model without the interaction coefficient, the absolute pointwise average bias varies between 10% and 15% of the dose range when the dose level of drug A varies between 0.8 and 1.0 compared to 2.5% and 5% of the dose range for models that include an interaction term. For scenario 4, the pointwise average bias for the logistic model with and without an interaction term are less than 7.5% of the dose range for almost all dose combinations.

Results
In Figure 4, the pointwise percent selection using a tolerance probability p = 0.1 is very poor when the two drugs are highly synergistic (β 3 = 40, 100) under scenarios 1 and 2 for the logistic model without an interaction coefficient. Under scenario 1, this percent selection varies between 80% and 100% when using all three models with interaction terms -logistic, probit, and comple-mentary log-log -compared to a 0% to 60% range when using a logistic model with β 3 = 0. A similar poor performance is observed under scenario 2 and to a lesser extent for scenario 3 where there is an overlap between the models when the dose level of drug A is less than 0.5. For dose combinations above 0.5, this pointwise percent selection drops considerably relative to models that include the interaction term. Under scenario 4 where the true dose-toxicity relationship is the six-parameter model, the pointwise percent selection using the logistic model with an interaction term is still higher relative to the model without an interaction parameter for most dose combinations and the extent of difference can be as high as 30% near dose level 0.15 of drug A. Therefore, we conclude that the efficiency of the estimate of the MTD curve as measured by the pointwise bias and percent selection is greatly reduced when omitting an interaction coefficient in the model.
For trial designs that recommend a single MTD, usually taken as the last dose combination given to the last patient in the trial or the dose that would have been given to the (n + 1)th hypothetical patient in the trial (Braun & Wang, 2010;Wages et al., 2011), Figure 5 shows the two-dimensional density plots of the last dose combination from all m = 2000 simulated trials under each scenario, using the logistic model with and without an interaction term. These plots show that using a model with an interaction term results in a much wider coverage area along the MTD curve relative to a model without an interaction coefficient. For example, under scenario 1, starting with density level 0.2, the coverage area for the logistic model with an interaction term is [0.1, 0.45] × [0.08, 0.32] compared to the dose combination range of [0.15, 0.32]×[0.14, 0.28] when omitting the interaction coefficient. The later model seems to recommend a dose that is concentrated near the middle of the curve with high probability. Therefore, if dose combinations away from the middle of the true MTD curve have high probability of efficacy relative to doses in the middle of the curve, then using models with no interaction term will more likely result in a failed phase II efficacy trial compared to models that include an interaction coefficient.

Conclusions
The use of drug combinations in cancer treatment allows targeting of different signaling pathways simultaneously, may reduce treatment resistance to therapy, clonal heterogeneity, and suppress cellular mechanisms associated with adaptive resistance. When several dose levels of these combinations are explored in early phase trials, more than one maximum tolerated doses can exist with possibly different levels of efficacy. This multiplicity problem is further amplified by the inherent small sample size of these trials, and this has led many authors to use parsimonious models for dose finding, especially in the absence of prior information about each drug when used individually. In fact, several authors recommended dropping the interaction term modeling synergism between the two drugs without affecting the performance of the design. They argued that precise estimation of the interaction coefficient is not possible due to the small sample size and omitting this term has negligible impact on the percent of dose recommendations at the end of the trial. These claims have been assessed for discrete dose combinations under one particular logistic model and a design that recommends one MTD by Mozgunov et al. (2021) and were found to be consistent with the previous recommendations.
We evaluated these assertions in the setting of drug combinations with continuous dose levels with a design that recommends an MTD curve at the end of the trial. Contrary to the findings of the previous authors, we found that when the two drugs are highly synergistic, not including an interaction term can result in a significantly higher percent of trials with a DLT rate exceeding θ + 0.05 and hence compromises the safety of the trial. In addition, efficiency of the estimated MTD curve as assessed by the pointwise average bias and percent selection is greatly reduced relative to a design that includes the interaction coefficient. We note that these findings were obtained when the true dose-toxicity model belongs to the family F(β 0 + β 1 x + β 2 y + β 3 xy) for different link functions and the six-parameter model (1 + α 1 x β 1 + α 2 y β 2 + α 3 (x β 1 y β 2 ) β 3 ) -1 , and the working model uses the logistic link function. In addition, the algorithm proposed by Tighiouart et al. (2017) was used for dose escalation and MTD curve estimation. Therefore, these conclusions may not apply to other forms of dose-toxicity models and dose allocation algorithms. We further showed that even for some designs that recommend a single MTD at the end of the trial, omitting an interaction term may result in a dose recommendation at the end of the trial with low probability of efficacy. In practice, for either discrete or continuous dose levels, we recommend that the statistician derives the operating characteristics in collaboration with the clinician using a model with and without an interaction term under a large class of scenarios that assume weak and very strong synergism between the two drugs.

Software
The repository https://github.com/jjimenezm1989/synergism-early-phase-continuous-dose-levels contains all the R codes used to produce the results presented in this paper.