 Methodology
 Open access
 Published:
A flexible doseresponse modeling framework based on continuous toxicity outcomes in phase I cancer clinical trials
Trials volume 24, Article number: 745 (2023)
Abstract
Background
The past few decades have seen remarkable developments in dosefinding designs for phase I cancer clinical trials. While many of these designs rely on a binary toxicity response, there is an increasing focus on leveraging continuous toxicity responses. A continuous toxicity response pertains to a quantitative measure represented by real numbers. A higher value corresponds not only to an elevated likelihood of side effects for patients but also to an increased probability of treatment efficacy. This relationship between toxicity and dose is often nonlinear, necessitating flexibility in the quest to find an optimal dose.
Methods
A flexible, fully Bayesian dosefinding design is proposed to capitalize on continuous toxicity information, operating under the assumption that the true shape of the dosetoxicity curve is nonlinear.
Results
We conduct simulations of clinical trials across varying scenarios of nonlinearity to evaluate the operational characteristics of the proposed design. Additionally, we apply the proposed design to a realworld problem to determine an optimal dose for a molecularly targeted agent.
Conclusions
Phase I cancer clinical trials, designed within a fully Bayesian framework with the utilization of continuous toxicity outcomes, offer an alternative approach to finding an optimal dose, providing unique benefits compared to trials designed based on binary toxicity outcomes.
Background
Phase I cancer clinical trials are a critical first step in the study of novel cancer therapeutic approaches. One of the primary goals of the phase I studies is to determine the dose of a new drug or therapeutic agent for use in subsequent phase II trials [1,2,3,4]. In these trials, one of the fundamental assumptions is that toxicity is a precondition for antitumor activity to eliminate fastgrowing cancer cells [5]. This means that patients must endure some degree of treatmentrelated toxicity to have a reasonable chance of a favorable response. More precisely, the purpose of the cancer phase I clinical trial is to estimate the maximum tolerated dose (MTD) of a new drug associated with an acceptable level of doselimiting toxicity (DLT) [6].
The use of modelbased adaptive clinical trial designs in phase I clinical trials has received much attention because it allows adaptations of trials and statistical designs of ongoing clinical trials [7,8,9]. Most designs developed in the literature consider binary toxicity outcomes, but recently, there has been increasing recognition of the need to identify the MTD by incorporating continuous toxicity information. One motivation for this is that some binary toxicity outcomes are obtained by dichotomizing continuous data, inevitably incurring losses of information and statistical power in clinical trial designs [10,11,12]. This convention imposes an additional burden of finding an optimal cutpoint for the continuous data and may require more subjects than if the endpoint variable were utilized in its original continuous form [13].
The other motivation is that various continuous measures of toxicity arise during phase I studies, which can be beneficial for examining doseresponse relationships more precisely. For example, if one aims to consider all grades of toxicities from multiple adverse events when allocating doses, it may be necessary to derive a continuous toxicity score, such as the normalized equivalent toxicity score [14, 15], or use a weighted average form as discussed by [16], and utilize such scores as the continuous toxicity response.
Additionally, in pursuit of efficiency and speed in drug development, trialists are increasingly relying on the use of realworld data to assess the potential risks and benefits of new drugs, and much of this data may be originally measured on a continuous scale [17]. One example is the drug concentration in plasma from patients [18, 19]. It is widely known that overly high anticancer drug concentrations may be a risk factor for many side effects, such as cytokine release syndrome [20,21,22], and therefore, the drug exposure may represent a toxic reaction, depending on what is known about the mechanism of action of the drug. Another example is a logarithmic transformation of the number of white blood cells, which serves as a continuous toxicity response [23]. Furthermore, any quantity measured in real numbers, as a measurable indicator of the severity of some toxicity status related to new drug exposure, may be used as a surrogate endpoint, provided that a higher value of the quantity leads to a higher probability of side effects for patients as a price for a higher chance of treatment effect [24, 25].
Over the past few decades, there has been a remarkable development in adaptive dosefinding designs for phase I studies using binary toxicity outcomes, such as the continual reassessment method [26], the escalation with overdose control (EWOC) [3], and the Bayesian logistic regression model [27], along with their several extensions [28,29,30,31,32]. Refer to [4] for a survey of these methods. A theoretical framework using the binary toxicity response can be found in [33].
However, the utilization of continuous toxicity outcomes for optimal dose finding has garnered relatively little attention compared to its counterpart based on binary toxicity outcomes. Since the earliest research work by Eichhorn in the 1970s [34], only a handful of research studies have been published [35], and there is no established theoretical framework available to develop a modelbased design in the literature. Recently, Chen et al. [15] proposed a variant version of EWOC, named EWOCNETS, based on a pseudoBernoulli likelihood [36], where the binary toxicity response is replaced with a continuous fractional response derived from a toxicity score system. One drawback of EWOCNETS is that the dosefinding does not follow the fully Bayesian paradigm, thus failing to describe the uncertainty of the MTD in a fully Bayesian manner. More recently, Lee et al. [16] introduced a fully Bayesian design based on constrained linear regression, called the twoparameter linear dosefinder. In this design, the authors aimed to leverage all grade information of toxicities according to the Common Toxicity Criteria for Adverse Events (CTCAE) [37] within a fully Bayesian framework. Although the design provides a fully Bayesian dosefinding algorithm, it falls short in describing the nonlinear shape of the dosetoxicity curve, and its application is confined to analyzing grade information based on CTCAE.
In this paper, our objective is to bridge the existing gap by introducing a novel, fully Bayesian dosefinding algorithm that provides flexibility in describing the dosetoxicity curve and wider applicability, based on continuous outcomes. The incorporation of a continuous response within our framework accommodates various scenarios, such as using a continuous toxicity score [16], measuring drug concentration in plasma from patients [18], a biomarker response from a molecularly targeted agent [38], and so on.
Methods
This section presents a general modeling framework and performance evaluation metrics for dosefinding designs that utilize continuous toxicity responses, which are widely applicable in the context of phase I cancer trials. Following this, we introduce a new fully Bayesian design aiming to estimate the MTD or an optimal dose under a specific setting where the doseresponse curve is nonlinear.
Continuous toxicity response
We start by characterizing some basic concepts that define a continuous toxicity response in the doseresponse modeling framework. To that end, we first need to “clinically” understand two fundamental concepts: MTD and DLT, which are central to dosefinding problems in phase I trials. The mathematical definitions of MTD and DLT will be discussed in the next subsection.
The aims of typical phase I oncology trials are to determine MTD, assess the safety and tolerability, and investigate the pharmacokinetics and pharmacodynamics of a new drug. The United States National Cancer Institute defines MTD as the highest dose of a drug that does not cause unacceptable side effects (Visit the website www.cancer.gov/ for the definition). The determination of the MTD is based on the occurrence of DLT. DLT refers to druginduced toxicity associated with side effects of a drug that are serious enough to prevent an increase in dose. A clinical trial protocol must specify the criteria for DLT, often defined as any severe or lifethreatening adverse event [4].
Throughout the paper, x represents a dose of a new anticancer drug. We denote Y(x) to represent the continuous toxicity response of a patient against dose x. Technically, the dose x is an input with a positive real number (that is, \(x\in (0,\infty )\)), and the response Y(x) is an output assuming a real number. Allowing the response Y(x) to attain a negative real number is important for the generality of the dosefinding problems because, in many pharmaceutical applications, Y(x) may represent a change from baseline or a logtransformation of some continuous measurement. Eventually, for each patient, the outcome of the trial is represented by the ordered pair \((x,Y(x)) \in (0,\infty ) \times \mathbb {R}\).
We denote an open interval \((x_{\text {min}}, x_{\text {max}}) \subset (0,\infty )\) to represent an admissible dose range that clinicians want to explore during the trial. The dose range is determined by clinicians and remains fixed during the trials. Although clinicians expect that MTD would belong to the interval \((x_{\text {min}}, x_{\text {max}})\) (i.e., \(\text {MTD} \in (x_{\text {min}}, x_{\text {max}})\)), MTD is an unknown quantity, and we do not know whether it will indeed fall within this interval (i.e., it could happen that \(\text {MTD} \notin (x_{\text {min}}, x_{\text {max}})\)). Statistically, MTD is the parameter of main interest in phase I studies. Given no toxicity information from any patient, all we can do is widen the range of the interval \((x_{\text {min}}, x_{\text {max}})\) to increase the probability of MTD falling within the dose range. This comes at the price of decreasing the power of a model, as commonly encountered in many statistical problems.
Phase I trial designs for cytotoxic agents are based on the assumptions that (a) the clinical benefit of the agent increases with increasing dose, (b) the toxicity of the agent increases with increasing dose, and (c) there is a dose with acceptable toxicity that offers clinical benefit [39]. With these assumptions in mind, we posit five assumptions that define the continuous toxicity response Y(x). Later on, we will see that these assumptions are essential in bringing the clinical concepts of MTD and DLT to a statistical modelbased design:

A1. A side effect caused by the toxicity response Y(x) is negligible for doses \(x\in (0,x_{\text {min}})\).

A2. A side effect caused by the toxicity response Y(x) is mild at the dose \(x_{\text {min}}\).

A3. A side effect caused by the toxicity response Y(x) becomes more and more serious as the dose increases over the interval \((x_{\text {min}}, x_{\text {max}})\).

A4. A side effect caused by the toxicity response Y(x) at the dose \(x_{\text {max}}\) is lifethreatening or close to death.

A5. A side effect caused by the toxicity response Y(x) is too fatal for doses \(x \in (x_{\text {max}},0)\).
The assumptions above are common features describing a toxicity response used in many phase I clinical trials [5, 27, 40]. Many binary toxicity outcomes may be obtained by dichotomizing a continuous toxicity outcome (\(\text {DLT if }\) \(Y(x) > \eta ;\) , \(\text {nonDLT otherwise}\), provided a constant \(\eta\)) [35]. Regarding the statements in the assumptions, the terms “mild,” “lifethreatening,” and “death,” indicating the severity of toxicity, can refer to the United States National Cancer Institute’s Common Toxicity Criteria (CTCAE) [37]. However, these terms can be generalized depending on the context of the therapeutic area or the safety guidelines from regulatory agencies for drug approval, as long as the monotonic relationship between the response Y(x) and dose x holds.
The following are some important implications based on the assumptions. The first assumption (A1) means that patients are expected to have no treatment effect when they are assigned a dose \(x\in (0,x_{\text {min}})\). Thus, doses \(x\in (0,x_{\text {min}})\) are not going to be explored in phase I clinical trials. The second assumption (A2) implies that the dose \(x_{\text {min}}\) is safe enough for many patients; hence, dose \(x_{1} = x_{\text {min}} + \varepsilon \in (x_{\text {min}},x_{\text {max}})\) with some small value \(\varepsilon >0\) can be used as an initial dose for the first patient. The value of \(\varepsilon\) may depend on the unit of the dose. The third assumption (A3), called a monotonic dosetoxicity assumption, is the basic principle underlying cytotoxic anticancer agents or a combination of a biologic with a cytotoxic drug being developed for cancer chemotherapy. The maximum dose \(x_{\text {max}}\) in the fourth assumption (A4) is the supreme dose of the range, and dose \(x_{\text {max}}  \varepsilon \in (x_{\text {min}},x_{\text {max}})\) with some small value \(\varepsilon >0\) can be explored as the highest dose in phase I studies. Similar to the first assumption, the fifth assumption (A5) means that the interval \((x_{\text {max}},\infty )\) is out of the range of the studies.
The maximum tolerated dose
In this subsection, we aim to mathematically define MTD and DLT based on a continuous toxicity response Y(x) assuming (A1)–(A5). This inevitably requires some probabilistic statements. The definition of MTD based on continuous toxicity response was first conceptualized by authors [34], and recently, Lee et al. [16] modernized the definition for a more practical phase I clinical trial, accommodating multiple adverse events in determining MTD.
For cytotoxic anticancer agents, one fundamental assumption is usually made: the probability of a toxicity response increases with dose. (This is associated with the third assumption (A3) of Y(x).) A mathematical description of this assumption is as follows. We assume that clinicians have some prior knowledge about a threshold value \(\eta > 0\), such that any value Y(x) greater than \(\eta\) may cause the occurrence of DLT in patients assigned with dose x. Intuitively, what clinicians want to control is the “probability” of the occurrence of DLT at dose x, denoted as \(\textbf{Pr}[Y(X) \ge \eta  X = x]\). Clinically, this probability represents the proportion of patients who experience DLT at dose x. With these formulations, the fundamental assumption implies that the probability \(\textbf{Pr}[Y(X) \ge \eta  X = x]\) monotonically increases with dose x. For the safety of patients, in most cases, this probability is controlled by upperbounding it, so that the number of patients with DLT can be probabilistically restricted. The following defines the MTD:
Definition 1
Given prespecified values \(\eta >0\), \(\gamma \in (0.5, 1)\), and an open interval \((x_{\text {min}}, x_{\text {max}})\) \(\subset (0,\infty )\), the MTD is defined to be the largest value of x which satisfies the following inequalities
Throughout the paper, we shall denote the MTD as \(\xi\). Then, the MTD \(\xi\) based on Definition 1 can be rewritten as follows:
Here, the variables \(\eta > 0\), \(\gamma \in (0.5, 1)\), and \((x_{\text {min}},x_{\text {max}})\) are prespecified by clinicians at the planning stage of the design.
The medical interpretations of the variables are as follows:

Maximum toxicity level \(\eta >0\): the toxicity level that defines the DLT (i.e., \(Y(X) \ge \eta\)) and nonDLT (i.e., \(Y(X) < \eta\)). The determination of the threshold value \(\eta\) depends on the specific applications.

Homogeneity constant \(\gamma \in (0.5,1)\): the degree of clinicians’ prior belief in the homogeneity of patients’ toxic reaction against doses to be assigned during phase I clinical trials. A higher value of \(\gamma\) implies a stronger homogeneity of the toxic reactions of patients.

Dose range \((x_{\text {min}},x_{\text {max}}) \subset (0,\infty )\): The minimum dose \(x_{\text {min}}\) and the maximum dose \(x_{\text {max}}\). These are typically inferred from preclinical studies.
We detail some modeling considerations concerning the variables. The clinical meaning of the maximum toxicity level \(\eta\) depends on the context of the applications using the continuous toxicity response Y(x). For example, if Y(x) is based on a continuous toxicity score as discussed by [16], then \(\eta\) represents the maximum toxicity score. If Y(x) represents the pharmacokinetic exposure of a new drug, then \(\eta\) is associated with the maximum tolerated concentration [18, 41].
On the other hand, the homogeneity constant \(\gamma\) may have a closer relationship with the characteristics of patients undergoing treatment with a specific agent. Specifically, it indicates whether these patients exhibit homogeneity or heterogeneity in their response to a dose x of the new agent [16]. Mathematically, as inferred from the inequality (2), the value \(\theta = 1\gamma\) serves as the least upper bound for the probability of DLT occurrences among patients across the dose range \((x_{\text {min}},x_{\text {max}})\). This constant, denoted as \(\theta = 1\gamma\) in [42], is known as the target toxicity level [43]. Clinically, target toxicity level \(\theta\) is the target probability of DLT at MTD, representing the acceptable likelihood of a patient experiencing a DLT at MTD [42]. Normally, \(\theta\) is set relatively high when the DLT is reversible or nonfatal condition, and low if it is lifethreatening [40]. As a default value, we recommend \(\theta = 0.01\) (or equivalently, \(\gamma = 0.99\)) when dealing with a completely novel agent [16, 34]. With this choice, the safety of patients is prioritized, allowing at most one out of a hundred patients to exhibit DLT. If the agent under consideration has been used before but with variations in schedule, route of administration, or concomitant drugs, the default value for \(\gamma\) is often not explicitly defined. Its determination becomes highly dependent on the specific therapeutic area. Some authors suggest that \(\gamma\) typically falls within the range of 0.6 to 0.9 (or equivalently, \(\theta\) ranging from 0.1 to 0.4). For more details, refer to the Chapter on Phase I Trials in [44], or page 40 in [45].
The dose range \((x_{\text {min}},x_{\text {max}})\) depends on the unit of the drug (e.g., gram, milligram, or microgram) and is mostly inferred from animal studies, metaanalysis, or previous clinical studies with similar drug molecules, etc.
Figure 1 provides a visual depiction of the MTD \(\xi\) (2) from two distinct perspectives of toxicity outcomes: continuous (left panel) and binary (right panel). In both panels, the xaxis represents the dose x, while the yaxis corresponds to the continuous and binary toxicity responses in the left and right panels, respectively. The binary responses are obtained by dichotomizing the continuous responses using the threshold value of the maximum toxicity level \(\eta\). On the left panel, the curve illustrates the dosetoxicity relationship, represented by a monotonic function. The right panel presents the probability of DLT at a given dose x, denoted as \(\textbf{Pr}[Y(X) \ge \eta  X = x]\) in Eq. (1). In reality, the precise shapes of these curves remain unknown. However, a crucial consideration in phase I cancer studies is that, due to the assumption of a monotonic relationship for anticancer drugs, the true curves are expected to monotonically increase over the dose range \((x_{\text {min}},x_{\text {max}})\). The greenshaded intervals represents safe doses that fulfill the condition on the righthand side of Eq. (2). In contrast, the red intervals consists of highly toxic doses that could potentially induce DLT in patients. Ultimately, the primary objective of phase I clinical trials is to estimate the MTD \(\xi\), which corresponds to the upper limit of the green interval.
A fully sequential adaptive design
The present subsection aims to provide an algorithmic description of a framework of adaptive dosefinding design for finding MTD \(\xi\) (2). We explain fully sequential adaptive design (FSAD), which is the default setting in this paper. FSAD is characterized as follows:

1. Patients are introduced to the trials individually and sequentially.

2. Each patient is assigned an optimal estimate of MTD \(\xi\) (2) based on the accumulated patients’ information at interim.
To operate FSAD, some essential ingredients are needed. Suppose that we have a total of N patients who can participate in a phase I clinical trial. Let \((x_i, y_i)\) denote the ordered pair (dose, continuous toxicity response against the dose) for the ith patient (\(i=1,\cdots ,N\)). Let \(\mathcal {F}_n = \{(x_i, y_i) \}_{i=1}^{n}\) represent the accrued information from n patients (\(n=1,\cdots ,N\)). Due to the accumulation of patients’ information during the trial, it holds that \(\mathcal {F}_1 \subset \mathcal {F}_2 \subset \cdots \subset \mathcal {F}_N\), where the notation \(\subset\) represents the subset relationship.
Finally, we need a dosefinding rule (also called a design adaptation rule [4]), which is defined as a mapping from the information space \(\varvec{\mathcal {F}}\) to the dose range \((x_{\text {min}}, x_{\text {max}})\):
Dosefinding rule \(\mathcal {D} (\cdot )\) (3) receives the cumulative information set \(\mathcal {F}_n\in \varvec{\mathcal {F}}\) from the first n patients as input, and prints out the dose \(x_{n+1} = \mathcal {D}(\mathcal {F}_n)\) for the \((n+1)\)th patient that belongs to the interval \((x_{\text {min}} , x_{\text {max}})\). The output \(x_{n+1}\) is an optimal dose for the \((n+1)\)th patient. Generally, good operating characteristics of adaptive designs is determined by dosefinding rule \(\mathcal {D} (\cdot )\) (3), which is the drive engine of adaptive clinical trial designs.
Algorithm 1 describes four steps to implement FSAD. A pictorial description is displayed in Fig. 2. Starting from the initial dose \(x_1 \in (x_{\text {min}}, x_{\text {max}})\), we observe the continuous toxicity response \(y_1 = Y(x_{1})\) for the first patient, which then forms the information set \(\mathcal {F}_1=\{(x_1,y_1) \}\). Based on \(\mathcal {F}_1\), we select an optimal dose for the second patient by \(x_2 = \mathcal {D}(\mathcal {F}_{1})\), completing the first cycle. The second cycle starts by introducing the dose \(x_2\) to the second patient, and we record his or her response \(y_2\), leading to the accrued information \(\mathcal {F}_2=\{(x_1,y_1),(x_2,y_2)\}\). We then find an optimal dose for the third patient, that is, \(x_3 = \mathcal {D}(\mathcal {F}_{2})\). This cycle continues until we reach the Nth patient, who will be assigned the dose \(x_N = \mathcal {D}_{\alpha }(\mathcal {F}_{N1})\) based on the accrued information \(\mathcal {F}_{N1} = \{(x_i, y_i) \}_{i=1}^{N1}\). The final estimate \(x_N\) is then regarded as an optimal estimate of the MTD \(\xi\) (2).
Normally, for safety, the initial dose \(x_{1}\) is set to a dose very close to the minimum dose \(x_{\text {min}}\) (i.e., \(x_{1} = x_{\text {min}} + \epsilon\) with sufficiently small \(\epsilon >0\)), and we expect that the sequence of doses \((x_{n})\) slowly converges to the targeted MTD \(\xi\) (2) as the trials progress. Here, the dose sequence \((x_{n})\) does not need to be monotonic (that is, the inequalities \(x_{1} \le x_{2} \le \cdots \le x_{N}\) are not required) because we may deescalate the dose if some previous doses are overly toxic.
In practice, patients are often treated in a cohort size of three. In this case, FSAD (Algorithm 1) can be can be easily generalized to cohort sequential adaptive design (CSAD). Suppose that there are in total N cohort groups and patients are enrolled in a cohort of size of C, where C is a fixed positive integer (if \(C = 1\), then CSAD reduces to FSAD; \(C = 3\) is often used for CSAD). Therefore, we have in total \(N\cdot C\) patients who can participate in a phase I clinical trial. Algorithm for CSAD can be obtained by (1) assigning dose \(x_{n}\) to C patients in the nth cohort; (2) recording Cdimensional vector of continuous toxicity responses \(y_{n,1:C}=(y_{n1},\cdots ,y_{nC})\) from C patients in the nth cohort; (3) formulating the information set \(\mathcal {F}_N = \{(x_i, y_{i,1:C}) \}_{i=1}^{N}\); and (4) finding an optimal dose for the \((n+1)\)th cohort based on the accumulated cohorts’ information at interim by \(x_{n+1} = \mathcal {D} (\mathcal {F}_n)\).
Diagnosis of a dosefinding rule
Now, the fundamental question is this: given a dosefinding rule \(\mathcal {D} (\cdot )\) (3) under FSAD (Algorithm 1), how can we evaluate the algorithm’s utility for actual phase I clinical trials? Due to the smallsample nature of the problem, perhaps the best approach is to examine this through conducting clinical trial simulations as follows. The basic idea here is to assess the “selfsequential learning ability” by replicating a large number of phase I clinical trials and observing their clinical operating characteristics. To achieve this, we begin by specifying three variables \((N, x_{1}, \xi _{0})\) as follows: setting the sample size N to be small enough (e.g., \(N=20\) or 25 patients); determining the initial dose \(x_{1}\) to be sufficiently close to the minimum dose \(x_{\text {min}}\); and fixing the MTD to a “true” value denoted as \(\xi _{0}\) (2). The subscript “0” is used to indicate the “truth.”
Next, we define a data generating distribution for simulating toxicity responses \(y_{n}\) given doses \(x_{n}\) (\(n=1,\cdots ,N\)). Conceptually, this distribution takes the form of a conditional distribution, denoted as \(P_{\xi }(yx)\), parameterized by the MTD \(\xi \in (x_{\text {min}} , x_{\text {max}})\). With this assumption, we can implement Step 2 in Algorithm 1 by recording the response \(y_{n} \sim P_{\xi _{0}}(yx_{n})\) (\(n=1,\cdots ,N\)). Here, we assume that the N responses (\(y_{1},\cdots ,y_{N}\)) are conditionally independent given the parameter \(\xi\) and covariates (\(x_{1},\cdots ,x_{N}\)). Steps 3 and 4 in Algorithm 1 remain unchanged. These modifications to Algorithm 1 operate on a single simulated phase I clinical trial, and the outcome will be the information set \(\mathcal {F}_N = \{(x_i, y_i) \}_{i=1}^{N}\) obtained from N patients. The dose assigned to the last patient, \(x_N\), is then considered the final estimate of the MTD.
The distribution \(P_{\xi }(yx)\) describes the relationship between toxicity and dose. In most modelbased designs, the underlying distribution \(P_{\xi }(yx)\) is assumed to be parsimonious; otherwise, the doses for the first few patients could be suboptimal due to the small sample size. This may violate individual ethics—doing what is best for current patients in the trial. For that reason, the Bernoulli distribution is conventionally used if y represents a binary toxicity response [3, 26, 27], and the Gaussian distribution is often used if y is a continuous toxicity response [16, 34], while allowing only few number of parameters to describe some dynamics (e.g., slope, intercept, curvature, interpatient variability, etc.) of the dosetoxicity relationship. Therefore, one may rewrite the distribution \(P_{\xi }(yx)\) more concretely as \(P_{(\xi ,\phi )}(yx)\), where \(\phi\) represents some statistical nuisance parameter(s) describing such dynamics. However, for readability, we will keep using the notation \(P_{\xi }(yx)\) until this subsection.
The eventual success of modelbased dosefinding design relies on a dosefinding rule \(\mathcal {D}(\cdot )\) (3) that enables the design to generate the sequence of doses \((x_{n})\) by cycling through the steps in Algorithm 1. We hope that the sequence \((x_{n})\) converges to the true MTD \(\xi _{0}\) as n grows. Inappropriate choices of the dosefinding rule will lead to very slow convergence of the sequence or, in the worst scenario, it will never converge to the truth.
With cnth patient, \(\mathcal {F}_n = \{(x_i, y_i) \}_{i=1}^{n}\), we can generate the following four sequences of performance metrics to assess the utility of a proposed dosefinding rule \(\mathcal {D} (\cdot ): \varvec{\mathcal {F}} \longrightarrow (x_{\text {min}} , x_{\text {max}})\):

Number of patients with DLT (NPD): \(\text {NPD}(n) = \sum _{i=1}^{n} {\textbf {1}}( y_{i} \ge \eta )\).

Number of patients overdosed (NPO): \(\text {NPO}(n) =\sum _{i=1}^{n} {\textbf {1}}(x_{i}>\xi _{0})\).

Bias to MTD (BTM): \(\text {BTM}(n) =x_{n}  \xi _{0}\).

Square root of the relative mean squared error (RMSE): \(\text {RMSE}(n) = x_{n}  \xi _{0}/\xi _{0}\).
The notation \({\textbf {1}}(\cdot )\) represents the indicator function. NPD and NPO are intended to measure the safety of a treatment as drug dose, whereas BTM and RMSE evaluate estimation accuracy of a dosefinding rule. Specifically, NPD, NPO, and BTM are metrics tailored to phase I studies, while RMSE is generally reported when the parameter of interest is positive [46]. Generally, smaller values of NPD and NPO indicate better clinical safety. As for BTM, a smaller negative value is desired because we normally expect \(x_{n}\) to be smaller than the true MTD \(\xi _{0}\) but close enough to \(\xi _{0}\). Lastly, the smaller the square root of the RMSE, the better the estimation accuracy. (In this paper, we will simply denote the metric as “RMSE” instead of “SRMSE” to avoid lengthy notation.)
By evaluating the four metrics  \(\text {NPD}(n)\), \(\text {NPO}(n)\), \(\text {BTM}(n)\), and \(\text {RMSE}(n)\)  at each value of \(n=1,2,3,\cdots ,\infty\), we can access asymptotic aspects of a dosefinding rule \(\mathcal {D} (\cdot )\). According to the law of large numbers, \(\text {NPD}(n)/n\) and \(\text {NPO}(n)/n\) converge almost surely to the probability of DLT (that is, \(\textbf{Pr}[Y > \eta ]\)) and the probability of overdosing event (that is, \(\textbf{Pr}[X > \xi _{0}]\)) as n goes to infinity, respectively. By using the law of total expectation, we have \(\textbf{Pr}[Y> \eta ] = \mathbb {E}[{\textbf {1}}(Y> \eta )]= \mathbb {E}[\mathbb {E}[{\textbf {1}}(Y(X)> \eta )X]] = \mathbb {E}[\textbf{Pr}[Y(X) > \eta ]X]] \le \mathbb {E}[1 \gamma X] =1 \gamma\). Therefore, it holds \(\text {lim}_{n \rightarrow \infty } \text {NPD}(n)/n \le 1 \gamma\) almost surely. This implies that the proportion of patients with DLT is asymptotically controlled by the homogeneity constant \(\gamma\), or equivalently, target toxicity level \(1\gamma\). For example, by specifying \(\gamma = 0.9\), maximally one patient out of ten patients may show DLT in the long run. Later, we will see that the probability of overdosing event, \(\textbf{Pr}[X > \xi ]\), can be also controlled by clinicians by using “feasibility bound,” but unlike the controlling mechanism of the probability of DLT, \(\textbf{Pr}[Y > \eta ]\), done by \(\gamma\), the probability of overdosing event can be controlled in a Bayesian way. As for the two accuracy metrics, the sequences \((\text {BTM}(n))\) and \((\text {RMSE}(n))\) should converge to 0 as n grows; a dosefinding rule with this property is said to have a consistency property [34].
Due to the small sample nature of phase I studies, a practically useful diagnosis of a dosefinding rule \(\mathcal {D} (\cdot )\) can be performed by crosssectional analyses of the four metrics evaluated at a certain number \(n=N\), where N is a small natural number, say \(1530\), and by checking the metrics via replicated numerical experiments under diverse scenarios by varying design parameters \(\gamma\), \(\eta\), and \((x_{\text {min}},x_{\text {max}})\) that were used in defining the MTD and DLT. This way, we can assess whether the rule \(\mathcal {D} (\cdot )\) is robust enough to handle variations in environmental variables and can be used in real phase I studies across various therapeutic areas. Particularly, a dosefinding rule \(\mathcal {D} (\cdot )\) with a positive BTM (that is, \(x_{N} > \xi _{0}\)) is inappropriate for actual phase I clinical trials because it is highly likely that overly toxic doses would be suggested by the underlying model.
Nonlinear dosetoxicity curve
Let Y(x) denote continuous toxicity response, assuming the validity of five assumptions (A1)–(A5). Additionally, we assume that the patients’ responses at the minimum dose \(x_{\text {min}}\) have been appropriately transformed and centered around zero, ensuring \(\mathbb {E}[Y(x)x = x_{\text {min}}]=0\). In this paper, we further assume a nonlinear correlation between toxicity and dose, expressed as follows:
where \(\beta\), \(\nu\), and \(\sigma\) are positive real numbers. Here, \(\mathcal {N}(\mu ,\sigma ^{2})\) denotes the Gaussian distribution with mean \(\mu\) and standard deviation \(\sigma\). The parameters in the nonlinear Eq. (4) have specific interpretations. Specifically, \(\beta\) represents the slope of the dosetoxicity curve, and when \(\nu =1\), it represents the rate of the increment of toxicity response per unit of dose. On the other hand, \(\nu\) represents the nonlinearity parameter, with a higher value indicating a greater curvature of the dosetoxicity curve. Finally, \(\sigma\) represents the standard deviation of the toxic reaction of patients at a given dose x, describing interpatient variability.
Note that the nonlinear regression (4) reduces to a simple linear regression when \(\nu = 1\). This simple linear regression was studied by [16, 34] for the application of phase I clinical trials. Specifically, Eichhorn and Zacks [34] discussed a Bayesian estimation of the slope parameter \(\beta\), while fixing the standard deviation \(\sigma\). On the other hand, Lee et al. [16] aimed to estimate both parameters in a fully Bayesian way. Estimating \(\sigma\) is crucial because phase I cancer trials might enroll terminal cancer patients with different types of malignant tumors at various disease stages; hence, the patient population is usually heterogeneous [47]. Both of the previous studies assumed that the toxicity response is linear in dose, lacking the flexibility to describe the toxicitydose curve when the true curve is nonlinear. In this paper, we extend the previous research by introducing the nonlinear parameter \(\nu\) and estimating it in a fully Bayesian way.
A closedform expression of MTD \(\xi\) (2) can be derived by providing some constraints on the parameters \((\beta ,\nu ,\sigma )\):
Theorem 1
Consider a nonlinear regression to describe the relationship between toxicity and dose: \(y \sim \mathcal {N}(\beta (x  x_{\text {min}})^{\nu },\sigma ^{2})\). Suppose that the slope and nonlinearity parameters and standard deviation are restricted by \(\beta > (\eta  \sigma \Phi ^{1}(\gamma ))/(x_{\text {max}}  x_{\text {min}})^{\nu }\), \(\nu > 0\), and \(0< \sigma < \eta /\Phi ^{1}(\gamma )\). Then, MTD \(\xi\) (2) is given by
Proof
See the Appendix.
Figure 3 is a pictorial description of the mean part of the dosetoxicity curve (4), that is, \(f(x) = \mathbb {E}[Y(X)X=x] = \beta (x  x_{\text {min}})^{\nu }\). The curve is convex if \(\nu > 1\) (as shown in Fig. 3), line if \(\nu = 1\), and concave if \(0<\nu <1\). The maximum toxicity level \(\eta\) divides the dosetoxicity plane into the DLT region \(\{(x,y)  x \in (x_{\text {min}},x_{\text {max}} ), y \ge \eta \}\) and the nonDLT region \(\{(x,y)  x \in (x_{\text {min}},x_{\text {max}} ), 0< y < \eta \}\). The two regions, \(\{(x,y)  x \in (0, x_{\text {min}}], y \le 0 \}\) and \(\{(x,y)  x \in [x_{\text {max}},\infty ), y \ge \beta (x_{\text {max}}  x_{\text {min}} )^{\nu }\}\), represent subtherapeutic and overly toxic areas, respectively, associated with assumptions (A1) and (A5). Normally, clinicians believe that MTD \(\xi\) (5) belongs to the dose range \((x_{\text {min}},x_{\text {max}})\). By using elementary calculus, the solution of the equation \(f(x) = \eta\) is \(x_{\eta } = x_{\text {min}} + (\eta /\beta )^{1/\nu }\), represented as the red vertical dashed line in the panel. Therefore, MTD \(\xi\) (5) is located on the left side of the point \(x_{\eta }\) since it always holds \(\sigma \Phi ^{1}(\gamma ) > 0\) and \(\gamma \in (0.5,1)\): see the violet vertical dashed line.
Threeparameter nonlinear dosefinder (3PND)
Due to the smallsample, sequential nature of phase I clinical trials, the Bayesian framework is preferred for accurately estimating the MTD \(\xi\) (5) than frequentist framework in the literature [4, 48, 49]. From a Bayesian viewpoint, the parameters \(\beta\), \(\nu\), and \(\sigma\) are considered random variables, making MTD \(\xi\) a random variable as well. This means that the uncertainty underlying the estimation of the MTD \(\xi\) can be probabilistically described once an appropriate prior \(\pi (\beta , \nu , \sigma )\) is specified. This advantage of Bayesian methodology in phase I clinical trials does not exist in traditional rulebased designs [50].
We propose a threeparameter nonlinear dosefinder (3PND) as a fully Bayesian modelbased design. Given the information set \(\mathcal {F}_n = \{(x_i, y_i) \}_{i=1}^{n}\), the hierarchy of 3PND is given as follows:
where \(l(\sigma ,\nu ) = \{\eta  \sigma \Phi ^{1}(\gamma )\}/(x_{\text {max}}  x_{\text {min}})^{\nu }\) and \(u(\sigma ,\nu ) =\{\eta /(x_{\text {max}}  x_{\text {min}})^{\nu }\} + \sigma \Phi ^{1}(\gamma )\). \(\Phi (z)\) is the cumulative distribution function of the standard normal distribution. Notation \(\mathcal {U}(l,u)\) represents the uniform distribution supported on open interval (l, u). \(\mathcal {C}^{+}(0,1) \mathcal {I}_{(0, d)} \propto \{1/(1+z^2) \}\cdot \mathcal {I}_{(0, d)}\) represents the unit scaled halfCauchy distribution truncated on the interval (0, d). \(log \mathcal {N}(0,\delta ^{2})\) represents the lognormal distribution with a unit median and a scale \(\delta >0\). See the Appendix for the detail of posterior computation.
The prior distribution \(\pi (\beta , \nu , \sigma )\) (7)–(8) is conceptually weakly informative, imposing minimal restrictions on the prior, while assuring that the MTD \(\xi\) lies in the dose range \((x_{\text {min}}, x_{\text {max}})\) with probability 1. Specifically, the prior is constructed as follows: first, we assume a flat prior for \(\beta\) given \(\nu\) and \(\sigma\) (i.e., \(\pi (\beta \nu , \sigma ) \propto 1\), supported on \((\infty , \infty )\)), a unitscaled halfCauchy prior for \(\sigma\) (i.e., \(\pi (\sigma ) \propto 1/(1+\sigma ^2)\), supported on \((0, \infty )\)), and a lognormal prior for \(\nu\) (i.e., \(\pi (\nu ) = \log \mathcal {N}(0,\delta ^{2})\), supported on \((0, \infty )\)), centered around one, with the scale hyperparameter \(\delta > 0\). After that, we restrict the support of the joint prior \(\pi (\beta , \nu , \sigma )\) to ensure that the MTD \(\xi\) (5) belongs to the dose range \((x_{\text {min}}, x_{\text {max}})\). Note that the uniform and halfCauchy distributions are considered noninformative and weakly informative priors, respectively. Lognormal distribution is a subexponential distribution, a class of heavytailed distributions studied by Lee [51] for a smallsample problem. Its tailheaviness provides the flexibility of nonlinearity of the dosetoxicity curve. Default values for \(\delta\) are \(\delta = 0.1\) or 0.5. In real applications, a plausible value for \(\delta\) can be also chosen via sensitivity analysis. The higher the value of \(\delta\), the stronger the prior guess on the nonlinearity of the dosetoxicity relationship. Under the prior formulation, the values of \(\delta\), \(\eta\), \(\gamma\), \(x_{\text {min}}\), and \(x_{\text {max}}\), used in defining the MTD (1), are introduced as the hyperparameters of the 3PND.
Figure 4 displays a directed asymmetric graphical (DAG) model representation of the 3PND. Following the grammar of the graphical model (Chapter 8 of [52]), the circled variables indicate stochastic variables, while the observed ones are additionally colored in gray. Nonstochastic quantities are uncircled. The arrows indicate the conditional dependency between the variables.
Some difficulty in prior elicitation is briefly discussed. As noted in the guidance for the use of Bayesian statistics provided by the Food and Drug Administration [53], special care is needed in incorporating appropriate prior information in Bayesian modeling. Specifically, the prior information should allow the Bayesian model to be flexible and efficient in identifying any pattern during trials [9]. This implies that we should prevent a suggested Bayesian model from being overly dominated by the prior information, and to that end, imposing weaker prior information may be more reasonable. However, given the FSAD setting, patients are introduced to the trials sequentially and individually. Hence, allowing too weak prior information may lead to unstable parameter estimations, and it is highly likely that the first few patients are suboptimally dosed. The CSAD setting also suffers from this issue. This may contradict the individual ethics mentioned earlier. Therefore, a good model should retain a reasonable balance between these two competing requirements.
Theoretical guarantee of MTD estimation using 3PND
The fundamental principle of a dosefinding study design is to allocate each included subject to the current best estimate of MTD. The central question we may ask here is, “Will the proposed dosefinding algorithm discover the true MTD given a large number of patients, such as 100 patients or even more, say, 1,000 patients?” Although enrolling such a significant number of patients is not practical, one might seriously doubt the proposed design if it fails to detect the MTD even with such an extensive sample size.
In the present subsection, we provide a theoretical demonstration that the updated knowledge regarding the MTD, denoted as \(\xi\), progressively becomes more accurate and precise as the number of patients increases. This idealistic phenomenon is established through the concept of posterior consistency [54]. To elaborate, let \(\xi _0\) represent the true value of MTD \(\xi\). Our objective is to establish that, as long as the true value \(\xi _0\) resides within the dose range \((x_{\text {min}}, x_{\text {max}})\), the posterior distribution of MTD \(\xi\) becomes increasingly concentrated around \(\xi _0\) with an expanding sample size. In this case, we say that, “the posterior distribution of MTD \(\xi\) is consistent at \(\xi _0\)” [55]. To establish posterior consistency, we leverage Doob’s theorem [56] to provide a sufficient condition for the posterior consistency of the 3PND. This theorem asserts that “for any prior \(\pi\), the posterior is consistent at every value in the parameter space except, possibly, on a set of \(\pi\)measure zero.”
Theorem 2
Within the hierarchy of the 3PND (6) – (9), the posterior distribution of the MTD \(\xi\) is consistent at any value of \(\xi _0\) within the dose range \((x_{\text {min}},x_{\text {max}})\).
Proof
See the Appendix.
Theorem 2 guarantees that as long as the dose range \((x_{\text {min}}, x_{\text {max}})\) encompasses the targeted MTD \(\xi _0\), the posterior distribution of MTD \(\xi\) remains consistent at \(\xi _0\). In reality, since the true MTD \(\xi _0\) is unknown, it is advisable to select a dose range that is sufficiently broad to encompass the true value \(\xi _0\). However, the interval should not be overly wide as it may require a larger sample size than otherwise.
A dosefinding rule using 3PND
A dosefinding rule \(\mathcal {D}(\cdot )\) (3) is the driving engine to operate an adaptive dosefinding design (Algorithm 1). In the present subsection, we derive a dosefinding rule based on the 3PND (6) – (9) to satisfy the desiderata.
Suppose that we have accrued information from n patients \(\mathcal {F}_n = \{(x_i, y_i)\}_{i=1}^{n}\), which is the input of the mapping \(\mathcal {D}(\cdot )\) (3). The goal is now to select a dose for the \((n+1)\)th patient (denoted as \(x_{n+1}\)), having observed \(\mathcal {F}_n\). In the dose assigning procedure, our aim is to select an optimal dose \(x_{n+1}\) while controlling the following two quantities:

(i) The probability of the occurrence of DLT from patients

(ii) The posterior probability of the event of overdosing to patients
The controlling mechanism of the probability of DLT has already been taken into consideration in the formula of MTD \(\xi\) (5), which can be controlled by changing the homogeneity constant \(\gamma \in (0.5,1)\) (see Fig. 1). Noting from Eq. (2), \(1  \gamma\) represents the maximum proportion of patients with DLT at MTD. A reasonable value for \(\gamma\) may depend on the specific therapeutic area. If the anticancer agent is entirely novel, a recommended value for \(\gamma\) is 0.99 so that probabilistically, at most one patient out of a hundred patients shows the DLT status. Otherwise, the default value for \(\gamma\) is often not explicitly defined, as it highly depends on a certain therapeutic area. For example, values may range from 0.6 to 0.9.
On the other hand, to control the posterior probability of the overdosing event, we introduce a new variable, called (Bayesian) feasibility bound denoted as \(\alpha\) in the literature [16, 34]. To describe the concept, let us assume that \(\Pi _n(x) = \textbf{Pr}[\xi \le x  \mathcal {F}_n]\) denotes the posterior cumulative distribution function of MTD \(\xi\) given \(\mathcal {F}_n\). With a value \(\alpha \in (0,1)\) set by clinicians, the selected dose \(x_{n+1}\) is the posterior \(\alpha\)quartile for the MTD \(\xi\):
Inequality (10) implies that the posterior probability of the event of overdosing to the \((n+1)\)th patient is bounded above by the constant \(\alpha\). It is important to note that, unlike the homogeneity constant \(\gamma\), the feasibility bound \(\alpha\) is neither a hyperparameter nor a variable introduced to define MTD and DLT. Essentially, \(\alpha\) is a pure Bayesian apparatus to control the convergence rate of the dose sequence \((x_{n})\) toward the MTD. This quantilebased dose selection scheme has also been adopted in EWOC [40] and its extension [15].
In practice, the proper value of \(\alpha\) can be chosen through extensive simulations, depending on factors such as the disease, patients’ characteristics, and the number of patients to be enrolled. A lower value for \(\alpha\) leads to a more conservative dose escalation, allowing for a smaller jump size from the current dose \(x_n\) to the next dose \(x_{n+1}\) during the trials. On the other hand, a higher value for \(\alpha\) results in a more aggressive jump size from \(x_n\) to \(x_{n+1}\). We recommend using \(\alpha = 0.001\), 0.01, 0.05, or 0.1 as default values of the feasibility bound. Sometimes, we can use a variable feasibility bound \(\alpha\), starting with some small value of \(\alpha\), and as the trial progresses, \(\alpha\) increases in small increments [42]. In practice, an optimal value for the \(\alpha\) can be determined by extensive simulation experiments by trying different value of \(\alpha\), given values of \(\gamma\) and sample size N.
Now, we describe how to obtain \(x_{n+1}\) which satisfies the inequality (10). Because \(\xi\) is a continuous random variable, we have the theoretical expression \(x_{n+1} = \Pi _{n}^{1}(\alpha )\), where \(\Pi _{n}^{1}(\cdot )\) is the inverse function of \(\Pi _{n}(\cdot )\). However, the function \(\Pi _n(\cdot )\) is not known in a closedform distribution; hence, it is difficult to obtain \(x_{n+1}\) in a closedform solution. Instead, we resort to a Markov chain Monte Carlo (MCMC) algorithm [57] to approximate \(x_{n+1}\). Algorithm 2 describes the steps to produce the next dose \(x_{n+1}\) based on accrued patients’ information \(\mathcal {F}_{n}\) and feasibility bound \(\alpha\). More technically, the next dose \(x_{n+1}\) is an output of the function \(\mathcal {D}_{\alpha }(\cdot )\) evaluated at the input \(\mathcal {F}_{n}\): that is, \(x_{n+1} = \mathcal {D}_{\alpha }(\mathcal {F}_{n}) \in (x_{\text {min}}, x_{\text {max}})\) \((n=1,2, \cdots , N1)\). In the notation \(\mathcal {D}_{\alpha } (\cdot )\), the Greek letter \(\alpha\) is subscripted to emphasize that \(\alpha\) is fixed during trials.
Algorithm 2 represents the vanilla version, where the dosesearching process relies solely on the quantile Eq. (10). The development of an efficient sampling algorithm in Step 1 of Algorithm 2 is pivotal to the success of the dosefinding algorithm. We have devised a sampling algorithm that combines the Gibbs sampler [58], the slice sampler [59], and the elliptical slice sampler [60]. Detailed information can be found in the Appendix.
Results
To assess the operating characteristics and utility of 3PND in phase I cancer clinical trials, in this section, we conduct extensive simulation experiments and apply the proposed design to an optimal dosefinding problem using a doseresponse dataset resampled from Friedman et al. [38].
Simulation experiments
Outline of simulation experiments
To assess the method’s average behavior, we perform a simulation study. The general setup adopted here is similarly designed to the simulation experiments from [16, 34]. As we are mainly interested in evaluating the operating characteristics of the design, experiments are conducted based on the vanilla dosefinding algorithm described in Algorithm 2.
We compare four dosefinding algorithms that mainly differ in the enrollment schedule of patients with the same number of total patients. Additionally, we want to explore the “ExplorationExploitation Dilemma” [61]. As normally encountered in many sequential learning problems, there is a tradeoff between the exploration of new knowledge about MTD and the exploitation of old knowledge assuring patients’ safety. Optimal performance usually requires some balance between exploratory and exploitative behaviors. One of the benefits of the doseresponse modeling framework suggested in this paper is that this balance can be controlled by clinicians. We demonstrate this via simulation experiments.
The four algorithms are denoted as (1) FSAD (acc); (2) FSAD (std); (3) CSAD (acc); and (4) CSAD (std). The first two algorithms are fully sequential adaptive designs with accelerated and standard dosing strategies, respectively, whereas the third and fourth ones are cohort sequential adaptive designs with a cohort size of three \((C=3)\), with accelerated and standard dosing strategies, respectively.
Simulation setup
To simulate a phase I clinical trial, experimenter needs to specify the following five categories of variables: (1) the total number of patients N; (2) the maximum toxicity level \(\eta\), the homogeneity constant \(\gamma\), and the dose range \((x_{\text {min}},x_{\text {max}})\); (3) the feasibility bound \(\alpha\); (4) the initial dose \(x_1\) for the first patient; and (5) true values of MTD \(\xi _{0}\), standard deviation \(\sigma _{0}\), and nonlinearity parameter \(\nu _{0}\). Note that once the values of \(\eta\), \(\gamma\), \(x_{\text {min}}\), \(\xi _{0}\), \(\sigma _{0}\), and \(\nu _{0}\) are determined, then the true value of the slope parameter \(\beta _{0}\) is automatically derived through the formula (5), that is, \(\beta _{0} = (\eta  \sigma _{0} \Phi ^{1}(\gamma ))/(\xi _{0}  x_{\text {min}})^{\nu _{0}}\).
After specifying the aforementioned variables, we generate a continuous toxicity response \(y_{i} = Y(x_{i})\) according to the dosetoxicity curve (4), that is, \(y_{i} \sim P_{(\xi _{0},\nu _{0},\sigma _{0})}(yx_{i}) = \mathcal {N}(y\beta _{0} (x_{i}  x_{\text {min}})^{\nu _{0}}, \sigma _{0}^{2})\), and follow the procedure described the previous section to examine the operating characteristics of dosefinding algorithms. Provided a sample size N, the eventual outcome based on FSAD (and similarly for CSAD) is the sequence of doses \((x_{1},x_2, \cdots ,\) \(x_n, \cdots , x_N)\), and we evaluate the four metrics, \(\text {NPD}(N)\), \(\text {NPO}(N)\), \(\text {BTM}(N)\), and \(\text {RMSE}(N)\).
We assume that the scale of the continuous toxicity response \(y =Y(x)\) is aligned with the toxicity grade information in CTCAE with adjustment as similarly done by [16]. More specifically, on average, the values \(y = 0, 1,2,3\), and 4 indicate mild, moderate, severe, lifethreatening, and deathrelated toxicity of an adverse event against a dose x. These values are corresponding to CTCAE Grade 1, 2, 3, 4, and 5, respectively; see [37] for more detail about the generic symptom of adverse events. To check smallsample performance, we set the number of total patients to be \(N=18\), 30, and 45, and replicate 1000 times. Eventually, we report median values of the four metrics obtained from replications.
We experiment with two dosing strategies, accelerated and standard dosing strategies. To that end, we set the values of the homogeneity constant \(\gamma\) and the feasibility bound \(\alpha\) accordingly while fixing other variables. Accelerated design is more inclined to produce a fast and accurate estimation of MTD with a higher risk of DLT and overdosing, and standard design indicates the the acceleration is not used.
To summarize, the following are the variables we set for the simulation experiment:

Number of patients: \(N\in \{18, 30, 45\}\).

Variables specific to accelerated designs: \(\gamma =0.9\) and \(\alpha = 0.01\).

Variables specific to standard designs: \(\gamma =0.99\) and \(\alpha = 0.001\).

Variables describing the degree of nonlinearity of dosetoxicity curve: \(\nu _{0}\in \{0.6, 1, 1.4\}\).

Variables shared in accelerated and standard designs: \(\delta = 0.1\), \(\eta =3\), \((x_{\text {min}},x_{\text {max}})=(5,80)\), \(x_{1}=6\), \(\xi _{0}=50\), and \(\sigma _{0}=0.1\).
Under the above specification, a patient shows DLT if the patient’s toxicity response y is greater than \(\eta =3\) (that corresponds to CTCAE Grade 4), and a patient is overdosed if assigned dose x is greater than the true MTD \(\xi _{0}\).
Simulation results
Table 1 presents the results of simulation experiments based on the four designs. The second column in the table indicates the curvature of the true dosetoxicity curve: linear (\(\nu _{0}=1\)), strictly convex (\(\nu _{0}=1.4\)), and strictly concave (\(\nu _{0}=0.6\)). The bestperforming algorithm for each metric is highlighted in bold within each row. Overall, all designs exhibit excellent consistency in MTD estimation across various sample sizes and curvatures, as indicated by the presence of small negative values for BTMs alongside small positive values for RMSEs. It is noteworthy that accuracy improves with increasing sample size, leading to smaller RMSEs as N increases. This consistency is theoretically guaranteed, as demonstrated in Theorem 2.
Regarding safety, the standard designs provide a safer dosing strategy in terms of NPD and NPO when compared to the accelerated designs. Notably, all standard designs yield a median value of zero for both NPD and NPO. In terms of the accuracy of MTD estimation, the accelerated designs outperform the standard designs, as evidenced by their smaller RMSEs. The results underscore that 3PND possesses a desirable property for balancing safety (exploitation) and accuracy (exploration) through the homogeneity constants \(\gamma\) or feasibility bound \(\alpha\).
Application to find an optimal \(\text {O}^{6}\)BG dose
Outline of application
We have illustrated the use of continuous outcomes to indicate the degree of severity of toxicity in the context of oncology trials for finding the MTD for chemotoxic agents. One typical example is a continuous toxicity score measured based on grade information from multiple adverse events [14,15,16]. However, our modeling framework can be generalized to molecularly targeted agents that have little or no toxicity in the therapeutic dose range. Typical examples include biomarker responses of molecularly targeted agents (see [62] for more details). For example, such a biomarker response might be based on the level of a molecular target, or the change in the level of a target that suggests clinical promise. Basically, this generalization is possible because our modeling framework is based on the principle of the monotonic relationship between a continuous response and a dose.
In the following, we develop a redesign of a dose escalation trial based on a molecularly targeted endpoint, utilizing the doseresponse dataset resampled from Friedman et al. [38]. In their study, the authors conducted a phase I trial involving a molecularly targeted agent, \(\text {O}^{6}\)benzylguanine (\(\text {O}^{6}\)BG). The escalation strategy in this trial was grounded on the reduction of the target enzyme \(\text {O}^{6}\)alkylguanineDNA alkyltransferase (AGT) activity. The patient responses were initially recorded as continuous measurements—specifically, tumor AGT activity measured in fmol/mg—corresponding to discrete dose levels of 40, 60, 80, 100, and 120 \(\text {mg}/\text {m}^{2}\) (refer to Table 1 in [38]). The plan involved treating up to 13 patients at each dose level. Despite the continuous nature of the patients’ outcomes, the authors dichotomized these outcomes to facilitate the dose escalation procedure in the trial. Authors employed a rulebased design based on the depletion of the target AGT activity and concluded that a 100 \(\text {mg}/\text {m}^{2}\) dose of \(\text {O}^{6}\)BG is an optimal dose that will be used in another phase I trial.
Figure 5a displays the values of tumor AGT activity from 24 patients: 3, 3, 9, and 9 patients at doses 40, 60, 80, and 100 \(\text {mg}/\text {m}^{2}\), respectively. No patients in the trial were assigned to the dose level of 120 \(\text {mg}/\text {m}^{2}\). These values are as follows:

Tumor AGT activity (fmol/mg) of patients assigned to dose 40 \(\text {mg}/\text {m}^{2}\): 26.35, 42.00, 15.00.

Tumor AGT activity (fmol/mg) of patients assigned to dose 60 \(\text {mg}/\text {m}^{2}\): 23.00, 13.50, 11.00.

Tumor AGT activity (fmol/mg) of patients assigned to dose 80 \(\text {mg}/\text {m}^{2}\): 31.67, 8.00, 9.00, 14.50, 11.50, 7.00, 11.70, 9.03, 8.00.

Tumor AGT activity (fmol/mg) of patients assigned to dose 100 \(\text {mg}/\text {m}^{2}\): 4.07, 5.00, 8.70, 2.50, 4.07, 6.13, 3.60, 5.00, 5.00.
Among the 24 patients, there are 4 patients whose responses are smaller than 5 fmol/mg, resulting in undetectable tumor AGT activities. These values will later be used in the redesign based on 3PND and EWOC to construct a hypothetical dataset that resembles the actual trial.
In the following, we employ 3PND and EWOC to determine the optimal dose of \(\text {O}^{6}\)BG, resulting in undetectable tumor AGT levels (< 5 fmol/mg), as similarly researched by [35].
Finding an optimal \(\text {O}^{6}\)BG dose using 3PND
In order to apply 3PND for determining an optimal dose of \(\text {O}^{6}\)BG, we perform a transformation on the actual values of tumor AGT activity, yielding “60  Tumor AGT activity,” which will be used as the continuous response in 3PND. This transformation is crucial to establish a monotonic relationship between the continuous response y and the dose x to implement a modelbased design.
In Fig. 5b, we present the transformed dataset alongside the fitted curve (represented by a red dashed line) derived using the leastsquares method. Our assumption is that the true doseresponse function adheres to the form \(y = \beta \cdot (x  20)^{\nu }\). The estimated parameter values are \(\hat{\beta } = 10.57\) and \(\hat{\nu } = 0.37\). The nonlinear nature of the curve underscores the necessity of incorporating flexibility into the doseresponse relationship within our model design.
We apply a doseescalation rule based on 3PND (Algorithm 2), incorporating two options: “Discrete Dose Selection” and “Monotonicincreasing Dose Sequence.” (Details on the options are described in Discussion.) To simulate a hypothetical dataset for the doseescalation procedure and leverage the reported data from [38], we introduce perturbation errors sampled from a uniform distribution with specified lower and upper bounds: ±8, ±6, ±4, and ±2 (fmol/mg) corresponding to the dose level of 40, 60, 80, and 100 \(\text {mg}/\text {m}^{2}\). (Observing the actual dataset from [38], it is evident that the standard deviation of responses decreases as the dose increases.) These perturbations are correspondingly added to the resampled tumor AGT activity responses. This process constructs a dataset used to assess and refine our dosefinding strategy. We construct CSAD with cohort size of three with accelerated and standard doseescalation. Variables are specified as follows:

Number of patients and cohort size: \(N = 30\) and \(C = 3\).

Variables specific to accelerated designs: \(\alpha = 0.1\).

Variables specific to standard designs: \(\alpha = 0.05\).

Variables shared in accelerated and standard designs: \(\delta = 0.5\), \(\gamma = 0.6\), \(\eta =55\), \((x_{\text {min}},x_{\text {max}})=(20,140)\), and \(x_{1}=40\).
Figure 6 displays the results of the doseescalation based on 3PND designs. Panels a (i)–(ii) and panels b (i)–(ii) correspond to accelerated and standard designs, respectively. Blue curves in panels represent the posterior means of the dosetoxicity curve based on 3PND. Continuous responses and doses are reported in Table 2. The accelerated design ends up finding an optimal dose of 100 \(\text {mg}/\text {m}^{2}\), consistent with the findings from [38] and [35]. On the other hand, the standard design selects 80 \(\text {mg}/\text {m}^{2}\), which is more conservative than the accelerated design. No patients were overdosed in either design. In the accelerated design, 3 out of 30 patients reported tumor AGT activity smaller than 5 fmol/mg. In contrast, none of the patients in the standard design reported tumor AGT activity smaller than 5 fmol/mg.
The result of a redesign proposed by Wang and Ivanova [35] led to three patients being assigned a dose level of 120 \(\text {mg}/\text {m}^{2}\), resulting in overdosing. Since no patients were overdosed using the 3PND, our designs are safer than the similar research conducted previously. Overall, we observe that the 3PND with accelerated doseescalation shows promise in identifying an optimal \(\text {O}^{6}\)BG dose for the current scenario.
Comparison with finding an optimal \(\text {O}^{6}\)BG dose using EWOC
We proceed to apply EWOC [40] in order to utilize the binary responses. For the application of EWOC, we dichotomize the continuous responses (i.e., 60  tumor AGT activity) using a threshold value of 55 (fmol/mg). As a result of this dichotomization, the binary responses yield \(y=1\) when tumor AGT activity is less than 5 fmol/mg, and 0 otherwise. The result of dichotomization process on the continuous response is visually represented in Fig. 7.
To implement EWOC, we use a builtin R function ewoc_d1classical(.) within an R package library(ewoc), with the following setup: target rate of patients with undetectable tumor AGT activity to be theta = 0.4, the minimum dose to be min_dose = 20, the maximum dose to be max_dose = 140 , the initial dose to be first_dose = 40, and feasibility bound for accelerated design to be alpha = 0.1 and feasibility bound for standard design to be alpha = 0.05, with the specification of uniform prior distributions for the parameters involved in EWOC (See help(ewoc_d1classical) for details about the variable setup). To simulate a hypothetical dataset with number of patients 30 with cohort size of 3 (\(N=30\) and \(C=3\)), we use the same method adopted in the 3PND redesign, followed by dichotomization with threshold values of 55, as illustrated above. Eventually, this will lead to two EWOC designs with accelerated and standard dosefinding rules, such that the general setup will be similar to the two 3PND designs with accelerated and standard dosefinding rules, respectively, with the crucial difference being the data form (continuous versus binary responses).
Figure 8 displays the results of dosefinding based on EWOC designs. Panels a (i)–(ii) and panels b (i)–(ii) correspond to accelerated and standard designs, respectively. The blue curve in the panels represents the posterior mean of the dosetoxicity curve based on EWOC design. Binary responses and doses are reported in Table 3. Both the accelerated and standard designs end up identifying optimal doses of 100 \(\text {mg/m}^{2}\) and 80 \(\text {mg/m}^{2}\), respectively, consistent with the outcomes using 3PND. In the accelerated design, 4 out of 30 patients reported tumor AGT activity smaller than 5 fmol/mg. In the standard design, 3 out of 30 patients exhibited tumor AGT activity smaller than 5 fmol/mg. These patients are indicated by “1” in the binary response column of Table 3.
Overall, we observe that EWOC produces an optimal dose matching that of 3PND. However, the results show that EWOC required a larger number of patients with undetectable tumor AGT activity compared to 3PND. Our results imply that EWOC using binary responses might necessitate a larger patient sample than 3PND using continuous response to effectively learn the doseresponse curve.
Discussion
Modern dosefinding studies and designs are highly specific to individual clinical settings. For example, clinicians may wish to estimate the MTD more precisely by accommodating a particular shape of the dosetoxicity curve or by incorporating information from multiple groups in adaptive designs, etc. Additionally, implementing a dosefinding design for actual phase I cancer clinical trials involves several practical requirements. For instance, it is conventional to have discrete dose levels while allowing only one dose level increment for each patient or cohort. The current section presents several extensions of the basic modeling framework to accommodate such complexities and requirements, thus making the dose selection schemes more realistic.
Different nonlinear dosetoxicity curve
The likelihood part (6) of our design assumes that the relationship between dose and toxicity response can be captured by a power function. The mean function part can be replaced and generalized to a more complex growth curve shape, such as the Richards growth curve [63] or the Gompertz growth curve [64], each having a greater number of parameters than a power function used in 3PND. Such curves can describe a plateau on the dosetoxicity curve so that higher doses may not improve clinical benefit and toxicity does not necessarily increase with increased doses. While these curves would more dynamically describe the dosetoxicity curve than using the power curve, specifying prior distributions requires extensive research to ensure plausibility under small sample problems. Related research based on binary response can be found in [43].
Accommodation of multiple groups
The model discussed so far assumes that all patients are grouped together, implying that the MTD is expected to be the same for all of them. However, there are circumstances where patients need to be treated in different groups, leading to different MTDs across these groups. Simultaneously, given the limited sample size inherent to dosefinding problems in phase I studies, it is desirable to estimate these MTDs using a single model that allows information to be shared across different groups. The objective within this formulation is to maximize the therapeutic effect of a treatment for individual patients/groups, which is referred to as personalized maximum tolerated doses (MTDs) [65].
Our framework can be extended to estimate MTDs for multiple groups. To illustrate, we consider two groups, labeled as A and B, with the goal of estimating MTDs denoted as \(\text {MTD}_{\text {A}}\) and \(\text {MTD}_{\text {B}}\). For simplicity, we assume that the curvature parameter \(\nu\) is the same and known between the two groups. One approach is to allow the slope parameter to differ while maintaining a shared standard deviation between the groups. Consequently, the likelihood of the extended model will be \(y = \beta _{\text {A}}(x  x_{\text {min}})^{\nu } + \sigma \epsilon\) for group A and \(y = \beta _{\text {B}}(x  x_{\text {min}})^{\nu } + \sigma \epsilon\) for group B. Conditional priors for \(\beta _{\text {A}}\) and \(\beta _{\text {B}}\) given \(\sigma\), and the prior for \(\sigma\), remain the same as Eqs. (7) to (8), respectively. This implies \(\beta _{\text {A}}, \beta _{\text {B}}\sigma \sim \mathcal {U}(l(\sigma ,\nu ),u(\sigma ,\nu ))\) and \(\sigma \sim \pi (\sigma ) = \mathcal {C}^{+}(0,1) \mathcal {I}_{(0,\eta /\Phi ^{1}(\gamma ))}\). Sharing the same standard deviation \(\sigma\) between the two groups is crucial; otherwise, the estimation of MTDs becomes parallel, and information borrowing does not occur. With this suggested extension, explicit formulae for the MTDs are \(\text {MTD}_{\text {A}} = \xi (\beta _{\text {A}}, \sigma ) = x_{\text {min}} + [(\eta  \sigma \Phi ^{1}(\gamma ))/\beta _{\text {A}}]^{1/\nu }\) and \(\text {MTD}_{\text {B}} = \xi (\beta _{\text {B}}, \sigma ) = x_{\text {min}} + [(\eta  \sigma \Phi ^{1}(\gamma ))/\beta _{\text {B}}]^{1/\nu }\).
Discrete dose selection
Let \(\Omega = \{ d_{k} \in (x_{\text {min}},x_{\text {max}})  d_1< \cdots < d_K, \, k=1,\cdots ,K \}\) be the set of ordered dose levels selected for a trial. To adapt the vanilla algorithm (Algorithm 2) to incorporate discrete dose selection, we modify the algorithm’s third step by setting \(x_{n+1} = \mathcal {D}_{\alpha } (\mathcal {F}_{n}) = \text {argmin}_{k =1,\cdots ,K }  d_{k}  \Pi _{n}^{1}(\alpha ) \in \Omega\).
Stopping rule
A stopping rule can be integrated into Algorithm 1 by adding a break statement. This break statement will terminate the forloop within the algorithm. The specific condition for the break depends on the context of the toxicity response Y(x). For instance, if Y(x) represents a continuous toxicity score, clinicians can implement a break condition whenever the number of patients experiencing a CTCAE Grade \(\ge 4\) exceeds a prespecified threshold. This condition ensures that the trial is halted due to safety concerns.
Monotonicincreasing dose sequence
This option aims to ensure that the dose sequence \((x_{n})\) is monotonically increasing: \(x_{1} \le x_{2} \le \cdots \le x_{N}\). To achieve this, the third step in Algorithm 2 is replaced with \(x_{n+1} = \mathcal {D}_{\alpha } (\mathcal {F}_{n}) =\text {max} \{x_n, \Pi _{n}^{1}(\alpha )\}\). It is important to note that this option disables the deescalation of the dose. Therefore, the choice of the feasibility bound value \(\alpha\) becomes crucial.
Upper bounding dose increment
To safeguard patients from potential overdose, it is possible to impose an upper bound on the dose increment. For this purpose, in the third step of Algorithm 2, we employ \(x_{n+1} = \mathcal {D}_{\alpha } (\mathcal {F}_{n}) = \text {min}\{x_n + M, \Pi _{n}^{1}(\alpha )\}\), where \(M > 0\) is a constant. This design ensures that the dose increment remains constrained within an upper limit of M.
Repeated measurement of toxicity response
In cases where each patient is allowed to be assigned multiple doses, our modeling framework can be easily extended to a nonlinear mixedeffect modeling framework, or more generally, Bayesian hierarhical modeling framework [66]. Using 3PND as an example, the likelihood will now take the form of \(y_{ij}\sim \mathcal {N}(\beta _{i}(x_{ij}  x_{\text {min}})^{\nu _{i}},\sigma ^{2}),\) \((i=1,\cdots ,N;j = 1,\cdots ,M_{i})\), where \(y_{ij}\) represents the continuous toxicity response of the ith patient to the jth dose \(x_{ij}\). The parameters \(\beta _{i}\) and \(\nu _{i}\) denote the slope and nonlinearity parameter specific to the ith patient, while the standard deviation \(\sigma\) is shared across all N patients. A joint prior distribution for these parameters can be similarly given as the prior of 3PND (7) – (9) to ensure that the individualized MTD belong to the dose range, while individual model parameters \((\beta _{i},\nu _{i}), (i=1,\cdots ,N)\) follow a populationlevel model. Under this formulation, the objective is not only to estimate the individualized MTDs for each patient but also to estimate the population MTD, representing the MTD for all patients.
Conclusions
In this article, we reexamined the dosefinding problem, delving into its foundational aspects and utilizing continuous outcomes. We presented a comprehensive dosefinding analysis using our modeling framework and introduced the novel 3PND dosefinding algorithm. This algorithm, which estimates the nonlinearity parameter through a fully Bayesian approach, benefits from a modern sampling technique that enhances the efficiency of posterior computation (refer to the Appendix “Posterior computation” for details). Simulation experiments yielded results suggesting that our modeling framework demonstrates favorable trial operating characteristics and is wellsuited for actual phase I studies, ensuring patient safety and efficient convergence of dose sequences to identify the MTD even with reasonably small sample sizes.
Our modeling framework is not only applicable to determining the MTD for cytotoxic agents in firstinhuman studies but also proficient in identifying optimal doses for noncytotoxic agents and animal studies. This adaptability stems from its foundation in the regression framework.
One specific context where our modeling framework finds relevance is in the domain of molecularly targeted agents. Advances in these agents, characterized by minimal or no toxicity within the therapeutic dose range, have spurred investigations into optimal phase I trial designs [67]. In such scenarios, incorporating pharmacodynamic biomarkers becomes imperative to monitor the drug’s biological effects, particularly in trials involving molecularly targeted agents. Here, the determination of the administered dose is guided by target inhibition considerations rather than toxicity concerns, as elucidated in [62]. As demonstrated in the previous section, the 3PNDbased doseescalation rule can identify optimal doses by accounting for the nonlinearity of the doseresponse curve.
Another pertinent application lies in preclinical studies, particularly those involving repeated doses in individual patients [68]. In such instances, our modeling framework seamlessly extends to a nonlinear mixedeffect modeling approach [66]. This extension accommodates multiple doses for each patient, empowering the determination of personalized MTDs for every individual, as outlined in Discussion.
Availability of data and materials
The research used simulated data.
Code availability
Details of posterior computation are described in the Appendix.
Abbreviations
 MTD:

Maximum tolerated dose
 DLT:

Doselimiting toxicity
 CTCAE:

Common toxicity criteria for adverse events
 NPD:

Number of patients with doselimiting toxicity
 NPO:

Number of patients overdosed
 BTM:

Bias to MTD
 RMSE:

Square root of the relative mean squared error
 3PND:

Threeparameter nonlinear dosefinder
 EWOC:

Escalation with overdose control
References
Potter DM. Phase I studies of chemotherapeutic agents in cancer patients: a review of the designs. J Biopharm Stat. 2006;16(5):579–604.
Legedza AT, Ibrahim JG. Heterogeneity in phase I clinical trials: prior elicitation and computation using the continual reassessment method. Stat Med. 2001;20(6):867–82.
Tighiouart M, Rogatko A, Babb JS. Flexible Bayesian methods for cancer phase I clinical trials. Dose escalation with overdose control. Stat Med. 2005;24(14):2183–96.
Sverdlov O, Wong WK, Ryeznik Y, et al. Adaptive clinical trial designs for phase I cancer studies. Stat Surv. 2014;8:2–44.
O’Quigley J, Pepe M, Fisher L. Continual reassessment method: a practical design for phase 1 clinical trials in cancer. Biom. 1990;46(1):33–48.
Ratain MJ, Mick R, Schilsky RL, Siegler M. Statistical and ethical issues in the design and conduct of phase I and II clinical trials of new anticancer agents. JNCI: J Natl Cancer Inst. 1993;85(20):1637–43.
Berry DA. Adaptive clinical trials in oncology. Nat Rev Clin Oncol. 2012;9(4):199–207.
Chow SC, Chang M. Adaptive design methods in clinical trialsa review. Orphanet J Rare Dis. 2008;3(1):1–13.
Chow SC, Chang M. Adaptive design methods in clinical trials. Boca Raton: Chapman and Hall/CRC; 2011.
Schmitz S, Adams R, Walsh C. The use of continuous data versus binary data in MTC models: a case study in rheumatoid arthritis. BMC Med Res Methodol. 2012;12(1):1–17.
Royston P, Altman DG, Sauerbrei W. Dichotomizing continuous predictors in multiple regression: a bad idea. Stat Med. 2006;25(1):127–41.
Wason JM, Mander AP, Eisen TG. Reducing sample sizes in twostage phase II cancer trials by using continuous tumour shrinkage endpoints. Eur J Cancer. 2011;47(7):983–9.
Faraggi D, Simon R. A simulation study of crossvalidation for selecting an optimal cutpoint in univariate survival analysis. Stat Med. 1996;15(20):2203–13.
Yuan Z, Chappell R, Bailey H. The continual reassessment method for multiple toxicity grades: a Bayesian quasilikelihood approach. Biometrics. 2007;63(1):173–9.
Chen Z, Tighiouart M, Kowalski J. Dose escalation with overdose control using a quasicontinuous toxicity score in cancer Phase I clinical trials. Contemp Clin Trials. 2012;33(5):949–58.
Lee SY, Munafo A, Girard P, Goteti K. Optimization of dose selection using multiple surrogates of toxicity as a continuous variable in phase I cancer trial. Contemp Clin Trials. 2022;113:106657.
CorriganCuray J, Sacks L, Woodcock J. Realworld evidence and realworld data for evaluating drug safety and effectiveness. JAMA. 2018;320(9):867–8.
Hutchinson TH, Bögi C, Winter MJ, Owens JW. Benefits of the maximum tolerated dose (MTD) and maximum tolerated concentration (MTC) concept in aquatic toxicology. Aquat Toxicol. 2009;91(3):197–202.
Ursino M, Zohar S, Lentz F, Alberti C, Friede T, Stallard N, et al. Dosefinding methods for phase I clinical trials using pharmacokinetics in small populations. Biom J. 2017;59(4):804–25.
ShimabukuroVornhagen A, Gödel P, Subklewe M, Stemmler HJ, Schlößer HA, Schlaak M, et al. Cytokine release syndrome. J Immunother Cancer. 2018;6(1):1–14.
Tvedt THA, Vo AK, Bruserud Ø, Reikvam H. Cytokine release syndrome in the immunotherapy of hematological malignancies: the biology behind and possible clinical consequences. J Clin Med. 2021;10(21):5190.
Lin SK, Su SF, Pan CH. Higher plasma drug concentration in clozapinetreated schizophrenic patients with side effects of obsessive/compulsive symptoms. Ther Drug Monit. 2006;28(3):303–7.
Mick R, Ratain MJ. Modelguided determination of maximum tolerated dose in phase I clinical trials: evidence for increased precision. J Natl Cancer Inst. 1993;85(3):217–23.
Renwick AG, Walton K. The use of surrogate endpoints to assess potential toxicity in humans. Toxicol Lett. 2001;120(1–3):97–110.
Kemp R, Prasad V. Surrogate endpoints in oncology: when are they acceptable for regulatory and clinical decisions, and are they currently overused? BMC Med. 2017;15(1):1–7.
O’Quigley J, Shen LZ. Continual reassessment method: a likelihood approach. Biometrics. 1996;52(2):673–84.
Neuenschwander B, Branson M, Gsponer T. Critical aspects of the Bayesian approach to phase I cancer trials. Stat Med. 2008;27(13):2420–39.
Babb JS, Rogatko A. Patient specific dosing in a cancer phase I clinical trial. Stat Med. 2001;20(14):2079–90.
Rogatko A, Ghosh P, Vidakovic B, Tighiouart M. Patientspecific dose adjustment in the cancer clinical trial setting. Pharm Med. 2008;22(6):345–50.
Tighiouart M, Rogatko A. Number of patients per cohort and sample size considerations using dose escalation with overdose control. J Probab Stat. 2012;2012. https://www.hindawi.com/journals/jps/2012/692725/.
Tighiouart M, CookWiens G, Rogatko A. Incorporating a patient dichotomous characteristic in cancer phase I clinical trials using escalation with overdose control. J Probab Stat. 2012;2012. https://www.hindawi.com/journals/jps/2012/567819/.
Mauguen A, Le Deley M, Zohar S. Dosefinding approach for dose escalation with overdose control considering incomplete observations. Stat Med. 2011;30(13):1584–94.
O’Quigley J. Theoretical study of the continual reassessment method. J Stat Plan Infer. 2006;136(6):1765–80.
Eichhorn BH, Zacks S. Sequential search of an optimal dosage. I J Am Stat Assoc. 1973;68(343):594–8.
Wang Y, Ivanova A. Dose finding with continuous outcome in phase I oncology trials. Pharm Stat. 2015;14(2):102–7.
Gong G, Samaniego FJ. Pseudo maximum likelihood estimation: theory and applications. Ann Stat. 1981;9(4):861–9.
NIH. Common Toxicity Criteria for Adverse Events v5.0. 2017. https://ctep.cancer.gov/protocoldevelopment/electronic_applications/ctc.htm#ctc_50. Accessed 20 Nov 2023.
Friedman HS, Kokkinakis DM, Pluda J, Friedman AH, Cokgor I, Haglund MM, et al. Phase I trial of O6benzylguanine for patients undergoing surgery for malignant glioma. J Clin Oncol. 1998;16(11):3570–5.
Korn EL. Nontoxicity endpoints in phase I trial designs for targeted, noncytotoxic agents. Oxford University Press; 2004.
Babb J, Rogatko A, Zacks S. Cancer phase I clinical trials: efficient dose escalation with overdose control. Stat Med. 1998;17(10):1103–20.
Siriwatwechakul W. TemperatureSensitive Poly (Acrylamide) Hydrogels for Drug Delivery Applications. Sci Technol Asia. 2010;15(5):94–101.
Tighiouart M, Rogatko A. Dose finding with escalation with overdose control (EWOC) in cancer clinical trials. Stat Sci. 2010;25(2):217–26. https://projecteuclid.org/journals/statisticalscience/volume25/issue2/DoseFindingwithEscalationwithOverdoseControlEWOCinCancer/10.1214/10STS333.full.
Wheeler GM, Mander AP, Bedding A, Brock K, Cornelius V, Grieve AP, et al. How to design a dosefinding study using the continual reassessment method. BMC Med Res Methodol. 2019;19(1):1–15.
Armitage P, Colton T, et al. Encyclopedia of biostatistics. New York: J. Wiley; 1998.
Zhou T, Ji Y. Emerging Methods for Oncology Clinical Trials. CHANCE. 2020;33(3):39–48.
Fonseca TC, Ferreira MA, Migon HS. Objective Bayesian analysis for the Studentt regression model. Biometrika. 2008;95(2):325–33.
Chow SC, Liu Jp. Design and analysis of clinical trials: concepts and methodologies. vol. 507. Hoboken: John Wiley & Sons; 2008.
Pallmann P, Bedding AW, ChoodariOskooei B, Dimairo M, Flight L, Hampson LV, et al. Adaptive designs in clinical trials: why use them, and how to run and report them. BMC Med. 2018;16:1–15.
Rogatko A, Babb JS, Tighiouart M, Khuri FR, Hudes G. New paradigm in dosefinding trials: patientspecific dosing and beyond phase I. Clin Cancer Res. 2005;11(15):5342–6.
Thall PF, Simon R, Ellenberg SS. Twostage selection and testing designs for comparative clinical trials. Biometrika. 1988;75(2):303–10.
Lee SY. The use of a lognormal prior for the student tdistribution. Axioms. 2022;11(9):462.
Bishop CM, Nasrabadi NM. Pattern recognition and machine learning, vol. 4. New York: Springer; 2006.
Fda U. Guidance for the use of Bayesian statistics in medical device clinical trials. Guidance for industry and FDA staff US FDA Docket. 2010;2006D–0191:50.
Barron A, Schervish MJ, Wasserman L. The consistency of posterior distributions in nonparametric problems. Ann Stat. 1999;27(2):536–61.
Ghosal S. A review of consistency and convergence of posterior distribution. In: Varanashi Symposium in Bayesian Inference, Banaras Hindu University; 1997.
Doob JL. Application of the theory of martingales. Paris: Coll Int du CNRS; 1948. p. 22–8.
Robert C, Casella G. Monte Carlo statistical methods. New York: Springer Science & Business Media; 2013.
Lee SY. Gibbs sampler and coordinate ascent variational inference: a settheoretical review. Commun StatTheory Methods. 2021;51(6):1–21.
Neal RM, et al. Slice sampling. Ann Stat. 2003;31(3):705–67.
Murray I, Adams R, MacKay D. Elliptical slice sampling. In: Proceedings of the thirteenth international conference on artificial intelligence and statistics. JMLR Workshop and Conference Proceedings; 2010. p. 541–548.
BergerTal O, Nathan J, Meron E, Saltz D. The explorationexploitation dilemma: a multidisciplinary framework. PLoS ONE. 2014;9(4):95693.
Curtin NJ. DNA repair dysregulation from cancer driver to therapeutic target. Nat Rev Cancer. 2012;12(12):801–17.
Richards FJ. A flexible growth function for empirical use. J Exp Bot. 1959;10(2):290–301.
Gompertz B. On the nature of the function expressive of the law of human mortality, and on a new mode of determining the value of life contingencies. In a letter to Francis Baily, Esq. FRS &c. By Benjamin Gompertz, Esq. FR S. In: Abstracts of the Papers Printed in the Philosophical Transactions of the Royal Society of London. 2. The Royal Society London; 1833. p. 252–253.
Chen Z, Li Z, Zhuang R, Yuan Y, Kutner M, Owonikoko T, et al. Adaptive estimation of personalized maximum tolerated dose in cancer phase I clinical trials based on all toxicities and individual genomic profile. PLoS ONE. 2017;12(1):0170187.
Lee SY. Bayesian nonlinear models for repeated measurement data: an overview, implementation, and applications. Mathematics. 2022;10(6). https://www.mdpi.com/22277390/10/6/898.
Hunsberger S, Rubinstein LV, Dancey J, Korn EL. Dose escalation trial designs based on a molecularly targeted endpoint. Stat Med. 2005;24(14):2171–81.
Aston WJ, Hope DE, Nowak AK, Robinson BW, Lake RA, Lesterhuis WJ. A systematic investigation of the maximum tolerated dose of cytotoxic chemotherapy with and without supportive care in mice. BMC Cancer. 2017;17(1):1–10.
Casella G, George EI. Explaining the gibbs sampler. Am Stat. 1992;46(3):167–74.
Damlen P, Wakefield J, Walker S. Gibbs sampling for Bayesian nonconjugate and hierarchical models by using auxiliary variables. J R Stat Soc Ser B Stat Methodol. 1999;61(2):331–44.
Acknowledgements
We would like to thank two reviewers for their constructive comments.
Funding
The research received no funding.
Author information
Authors and Affiliations
Corresponding author
Ethics declarations
Ethics approval and consent to participate
Not applicable.
Consent for publication
Not applicable.
Competing interests
The authors declare that they have no competing interests.
Additional information
Publisher’s Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Appendix
Appendix
Proof of Theorem 1
We can simplify the inequality (1) as follows
where \(\Phi\) denotes the cumulative distribution function of standard normal distribution. By using elementary calculus on the last inequality (11), we have
The inequality (12) implies that the value \(\xi = x_{\text {min}} + [\{\eta  \sigma \Phi ^{1}(\gamma )\}/\beta ]^{1/\nu }\) is the largest value that satisfies the inequality (11). On the other hand, the following inequality holds due to the parameter constraints on the \((\beta ,\nu ,\sigma )\) (that is, \(\beta > (\eta  \sigma \Phi ^{1}(\gamma )/(x_{\text {max}}  x_{\text {min}})^{\nu }\), \(\nu > 0\), and \(0< \sigma < \eta /\Phi ^{1}(\gamma )\)), we have \(x_{\text {min}}< \xi < x_{\text {max}}\). The two derived inequalities implies that \(\xi = x_{\text {min}} + [\{\eta  \sigma \Phi ^{1}(\gamma )\}/\beta ]^{1/\nu }\) is the largest value that satisfies the inequality (2).
Proof of Theorem 2
To prove the theorem, we first show that the marginal prior density of MTD \(\xi\), denoted as \(\pi (\xi )\), is supported within the dose range \((x_{\text {min}}, x_{\text {max}})\), which encompasses the complete parameter space of MTD \(\xi\). To this end, we begin by defining \(Z = \xi \sigma ,\nu\) as the continuous random variable representing MTD \(\xi\), conditioned on the standard deviation \(\sigma\) and the nonlinearity parameter \(\nu\). We can derive the closedform expression of the cumulative distribution function of Z as follows:
It is straightforward to show that \(G_{Z}(x)\) has the following properties: (i) \(G_{Z}(x)\) is strictly increasing on the open interval \((x_{\text {min}}^{*}(\sigma ,\nu ) ,x_{\text {max}})\), (ii) \(G_{Z}(x_{\text {min}}^{*}(\sigma ,\nu )) = 0\), and (iii) \(G_{Z}(x_{\text {max}}) = 1\), where the left end point \(x_{\text {min}}^{*}(\sigma ,\nu )\) of the interval is given by
Therefore, the random variable \(Z = \xi \sigma ,\nu\) is supported on the interval \((x_{\text {min}}^{*}(\sigma ,\nu ) ,x_{\text {max}})\).
Note that the left end point \(x_{\text {min}}^{*}(\sigma ,\nu )\) is a function of \(\sigma\) and \(\nu\), whose support is given as \((0, \eta /\Phi ^{1}(\gamma )) \times (0,\infty )\). On the other hand, by using elementary calculus, we can show that the function \(x_{\text {min}}^{*}(\sigma ,\nu )\) (13) is monotonically deceasing and continuous on the interval \((0, \eta /\Phi ^{1}(\gamma ))\), and the following onesided limits hold: \(\text {lim}_{\sigma \rightarrow (\eta /\Phi ^{1}(\gamma ))^{}}x_{\text {min}}^{*}(\sigma ,\nu ) = x_{\text {min}} \text { and } \text {lim}_{\sigma \rightarrow 0^{+}}x_{\text {min}}^{*}(\sigma ,\nu ) = x_{\text {max}}\), for any fixed value of \(\nu >0\). Summing the results indicates that the marginal random variable \(\xi\) is supported on the open interval \((x_{\text {min}},x_{\text {max}})\).
Because the marginal prior of \(\xi\) is compactly supported on the open interval \((x_{\text {min}},x_{\text {max}})\), the measure zero set over the dose range \((x_{\text {min}},x_{\text {max}})\) with respect to the prior measure \(\pi (\xi )\) is the empty set. Thus, by the Doob’s theorem [56], the posterior distribution of \(\xi\) is consistent at any value \(\xi _0 \in (x_{\text {min}},x_{\text {max}})\).
Posterior computation
Gibbs sampler
Let \(\mathcal {F}_n = \{(x_i, y_i)\}_{i=1}^{n}\) denote the set of the n observations from the first n subjects, where \(y_i = Y(x_i)\) represents the continuous toxicity response of the ith patient against the assigned dose \(x_i\). We consider a vectorform of the 3PND (6) – (9): \({\textbf {y}}_n = {\textbf {z}}_n ^{\nu } \beta + \sigma \varvec{\epsilon }_{n}\), \(\varvec{\epsilon }_{n}\sim \mathcal {N}_{n}({\textbf {0}} ,{\textbf {I}}_{n})\) (6), where \({\textbf {y}}_{n} = (y_1, \cdots , y_n)^{\top }\) and \({\textbf {z}}_{n}^{\nu } = ((x_1  x_{\text {min}})^{\nu }, \cdots , (x_n  x_{\text {min}})^{\nu })^{\top }\) are ndimensional response and covariate vectors, respectively. Here, the \(\nu\) in the vector notation \({\textbf {z}}_{n}^{\nu }\) implies the componentwise exponent.
Now, our eventual goal is to simulate a sample from the joint posterior distribution
where \(l(\sigma ,\nu ) = \{\eta  \sigma \Phi ^{1}(\gamma )\}/(x_{\text {max}}  x_{\text {min}})^{\nu }\) and \(u(\sigma ,\nu ) =\{\eta /(x_{\text {max}}  x_{\text {min}})^{\nu }\} + \sigma \Phi ^{1}(\gamma )\).
The nominator of the joint density \(\pi (\beta , \nu , \sigma  \mathcal {F}_n )\) (14) is detailed as
where \(c(\beta ,\nu ) = ({\textbf {y}}_{n}  {\textbf {z}}_{n}^{\nu } \beta )^{\top } ({\textbf {y}}_{n}  {\textbf {z}}_{n}^{\nu } \beta )/2\). Because the joint density \(\pi (\beta ,\nu ,\sigma \mathcal {F}_n)\) cannot be obtained in a closedform distribution, we develop a MCMC algorithm to sample from the joint density. Algorithm 3 details a Gibbs sampler [58, 69] to sample from the density \(\pi (\beta ,\nu , \sigma \mathcal {F}_n)\) (15).
Note that the density \(\pi (\beta \sigma ,\nu , \mathcal {F}_{n})\) (16) in Step 1 is a truncated normal distribution supported on the interval \((l(\sigma ,\nu ),u(\sigma ,\nu ))\) with mean \(\mu _{n}(\nu )\) and variance \(\sigma _{n}^{2}(\nu )\). Using a statistical software R, one may use R package truncnorm to sample from it. To implement Step 2 and Step 3, one consideration is to use MetropolisHasting algorithm (MH) [57] due to the nonclosed form the full conditional densities, which may have an issue of sampling efficiency. Instead, in the next subsections we introduce modern MCMC sampling techniques to sample from the densities.
Slice sampler implementation in Step 2
We use the slice sampler [59] to implement Step 2. Slice sampler is popular for its ability to adapt the change of the step size based on the local property of the target density. A central idea of the slice sampler is to find a valid parameter expansion from the target density, and use the Gibbs sampler to the expanded target distribution [70]. See [59] for more detail.
To apply the slice sampler in the Step 2 (17), first, use the change of variable, \(\tau = \sigma ^{2}\), to get
where \(g(\tau ) = 1/(1 + \tau )\) and \(A(\beta ,\nu ) = ( \{h(\beta ,\nu )\}^2, \{\eta / \Phi ^{1}(\gamma )\}^2 )\). The notation \(\mathcal{I}\mathcal{G}(a,b)\) denotes the inverse gamma distribution with shape parameter a and scale parameter b. Note that the function \(u = g(\tau )\) is decreasing on the positive real line, and its inverse function is given by \(\tau = g^{1}(u) = (1u)/u\). Consider a joint density, \(\pi (\tau , u \beta ,\nu , \mathcal {F}_n ) \propto \mathcal{I}\mathcal{G} (\tau  n/2, c(\beta ,\nu ) ) \cdot \mathcal {I}_{ (0, g(\tau )) }(u) \cdot \mathcal {I}_{A(\beta ,\nu )}(\tau ).\). We can show that \(\int \pi (\tau , u \beta ,\nu , \mathcal {F}_n ) du = \pi (\tau  \beta ,\nu , \mathcal {F}_n)\), which implies that \(\pi (\tau , u \beta ,\nu , \mathcal {F}_n )\) is a valid parameter expansion of the density \(\pi (\tau  \beta ,\nu , \mathcal {F}_n )\) (19). Finally, execute the Gibbs sampler to the joint density \(\pi (\tau , u \beta ,\nu , \mathcal {F}_n )\), and disregard u, and transform \(\tau\) back to \(\sigma\) to get a realization from \(\pi (\sigma  \beta ,\nu , \mathcal {F}_n )\).
Algorithm 4 details the steps to implement the slice sampler:
Elliptical slice sampler implementation in Step 3
The current subsection illustrates how to sample from the conditional posterior density \(\pi ( \nu \beta ,\sigma , \mathcal {F}_{n})\) (18) by using the elliptical slice sampler (ESS) [60]. Conceptually, ESS and MH algorithms are similar in that both comprises two steps: proposal step and criterion step. A difference between the two algorithms arises in the criterion step. If a new candidate does not pass the criterion, then MH takes the current state as the next state: whereas ESS reproposes a new candidate until rejection does not take place, rendering the algorithm rejectionfree. To apply the ESS in the Step 3 [18], we start with the change of variable, \(\phi = \log \nu\), to get
where \(\mathcal {L}(\phi ) = \exp \ \{c(\beta ,e^{\phi })/\sigma ^2 \}\). Algorithm 5 details the steps to implement ESS:
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.
About this article
Cite this article
Lee, S.Y. A flexible doseresponse modeling framework based on continuous toxicity outcomes in phase I cancer clinical trials. Trials 24, 745 (2023). https://doi.org/10.1186/s13063023077930
Received:
Accepted:
Published:
DOI: https://doi.org/10.1186/s13063023077930