For PK modeling, we assume a one-compartmental model for simplicity, which treats the human body as a single, uniform compartment and administers drugs via intravenous (IV) injection. Furthermore, we assume a first-order process to describe the rate of drug absorption and elimination [17].

Let c i ( t ) denote the drug concentration of patient i ∈ { 1 , … , N } at time t after drug administration. The rate of a first-order process is d c i ( t ) / d t = − k i c i ( t ), with initial concentration c i ( 0 ) = d i / V i . This yields the drug plasma concentration for patient i at time t as (2.1) c i ( t | d i , V i , k i ) = d i V i e − k i t where d i is one of the D drug dosages administered to patient i, d i ∈ { 1 , … , D }, k i denotes the patient-specific elimination rate constant, and V i is the patient-specific volume of distribution. We also assume V i > 0 and k i > 0 follow some prior distributions (e.g., gamma or log-normal) with hyperparameters α and λ. In our earlier work PEDOOP [26], the priors for V i and k i can only take gamma distributions, a restriction we relax here.

Let X i j denote the observed drug concentration at time t j ∈ { 1 , … , T } for patient i. We assume X i j follows a log normal distribution: (2.2) log ( X i j ) = log ( c i ( t j ∣ d i , V i , k i ) ) + ϵ i j , and ϵ i j ∼ N ( 0 , σ 2 ) where ϵ i j represents random error.

As described in PEDOOP, we compute the area under the concentration-time curve (AUC) for patient i by integrating the drug concentration c i ( t | d i , V i , k i ) from 0 to ∞: (2.3) A U C i ( d i , V i , k i ) = ∫ 0 ∞ c i ( t | d i , V i , k i ) d t = d i V i k i .

Next, we incorporate the AUC of patient i as a covariate in modeling their toxicity probability. Specifically, let p i denote the toxicity probability for patient i and assume (2.4) p i ≜ p i ( V i , k i , β 0 , β 1 ) = l o g i t − 1 ( β 0 + β 1 × log A U C i ( d i , V i , k i ) ) = l o g i t − 1 ( β 0 + β 1 × log d i V i k i ) , where l o g i t − 1 ( x ) = 1 / ( 1 + e − x ). Let Y i denote the binary DLT outcome for patient i and Y i follows a Bernoulli distribution given by (2.5) Y i | V i , k i , β 0 , β 1 ∼ B e r n o u l l i ( p i ( V i , k i , β 0 , β 1 ) ) . As a result, the joint posterior distribution of all parameters θ = ( V i , k i , α V , λ V , α k , λ k , β 0 , β 1 , σ 2 ) given the observed data D I = { ( X i j , Y i , d i ) }, becomes π ( θ | D I ) ∝ ∏ i N ∏ j T 1 σ 2 π e { − ( log ( X i j ) − log ( c i ( t j ∣ d i , V i , k i ) ) ) 2 2 σ 2 } × ∏ i N p i Y i ( 1 − p i ) 1 − Y i × ∏ i N g ( V i | α V , λ V ) × g ( k i | α k , λ k ) × h ( β 0 ) × h ( β 1 ) × h ( α V ) × h ( λ V ) × h ( α k ) × h ( λ k ) × h ( σ 2 ) . Here, g ( · ) and h ( · ) denotes the corresponding prior distribution functions. We adopt the hyper-prior distributions from PEDOOP, given by (2.6) α V ∼ G a m m a ( 4 , 1 ) , λ V ∼ G a m m a ( 1 , 1 ) α k ∼ G a m m a ( 3 , 1 ) , λ k ∼ G a m m a ( 1 , 1 ) σ ∼ G a m m a ( 3 , 3 ) , β 0 ∼ N ( − 3 , 100 ) , β 1 ∼ log N ( − 1 , 2 ) . These priors have been shown to lead to weak information on the p i ’s [26]. Using the Markov chain Monte Carlo (MCMC) algorithm, we can generate posterior samples of θ = ( V i , k i , α V , λ V , α k , λ k , β 0 , β 1 , σ 2 ).

For the first c cohorts, each with a fixed number of patients, we follow the conventional dose allocation rules, that is, assigning the same dose level to patients within the same cohort. Starting from the lowest dose level, we use the dose-level toxicity probability p ˜ d (defined next) to determine the appropriate dose for the next cohort via posterior predictive inference.

Assume n patients have been enrolled in the trial so far. Let D I represent the cumulative data of stage I. Let Y n e w denote the DLT outcome of a new patient in the next cohort treated at the dose level Z n e w = d. The predictive toxicity probability p ˜ d is computed as (2.7) p ˜ d ≜ P r ( Y n e w = 1 | D I , Z n e w = d ) = ∫ P r ( Y n e w = 1 | V n e w , k n e w , β 0 , β 1 , Z n e w = d ) × π ( V n e w , k n e w , β 0 , β 1 | D I ) d V n e w d k n e w d β 0 d β 1 = ∫ l o g i t − 1 ( β 0 + β 1 × log d V n e w k n e w ) × π ( V n e w , k n e w , β 0 , β 1 | D I ) d V n e w d k n e w d β 0 d β 1 ≈ 1 M ∑ m = 1 M l o g i t − 1 ( β 0 ( m ) + β 1 ( m ) l o g d V n e w ( m ) k n e w ( m ) ) , where notation ( m ) denotes the mth sample from the MCMC simulation for the joint posterior π ( θ | D I ), and V n e w ( m ) = ∑ i = 1 n V i ( m ) n and k n e w ( m ) = ∑ i = 1 n k i ( m ) n for the new patients in the next cohort. When the trial just starts, D I is an empty set and π ( V n e w , k n e w , β 0 , β 1 | D I ) becomes the prior distribution of V n e w and k n e w .

Following the CRM [14], we choose the dose d ˜ for the next cohort where (2.8) d ˜ = arg min d | p ˜ d − p T | and p T is a prespecified target toxicity rate (e.g., 0.3).

In addition, we apply the following rules to enhance the efficiency and safety of the trial. •

Speed-up rule: If no DLT outcomes are observed in any tested dose levels, escalate to the next higher dose for the next cohort to accelerate exploration.

Once all patients in the first c cohorts have been treated, we can select the MTD for the training stage. To determine the MTD, we suggest comparing the posterior summary statistic of the dose-level toxicity probability p ˆ d with the target toxicity probability p T , and selecting the dose that minimizes the difference. As described in Safety rule 2, assuming a beta prior and a binomial distribution for DLT outcomes, the posterior mean of the toxicity probability at dose d is: p ˆ d = Y d + 0.05 n d + 0.1

Note that, to determine the MTD, p ˆ d is transformed via isotonic regression to p ˆ ∗ d using the pool adjacent violators algorithm [16], which guarantees the monotonicity of the dose-toxicity relationship.

Consequently, we select the dose as the MTD for stage I that satisfies (2.10) d M T D = arg min d | p ˆ ∗ d − p T | from among the tested doses that are deemed safe using criterion (2.9). Since p ˆ d ∗ can be tied, we choose the MTD by the following rule:

(1) If p ˆ d M T D > p T , the MTD is the lowest dose among tied doses.

(2) If p ˆ d M T D ≤ p T , the MTD is the highest dose among tied doses.

In precision dose-finding, each patient receives a predicted MTD using the posterior information learned from stage I. We recommend that 30% of the total number of patients be treated under the precision dose-finding scheme, or that the training sample size be set to at least 21, which is typical for phase I trials. Instead of following the dose escalation rules introduced in Section 2.2 to determine dosage, we administer the dose based on each patient’s individual toxicity profile.

Specifically, for each patient, we estimate the predicted toxicity probability of patient i if s/he is treated at dose d by (2.11) p ˆ i , d = l o g i t − 1 ( β ˆ 0 + β ˆ 1 × log d V ˆ i k ˆ i ) , where β ˆ 0 and β ˆ 1 are the posterior means of β 0 and β 1 from the first stage dose-finding, respectively. ( V ˆ i , k ˆ i ) represent pre-dose predicted PK parameters for patient i. In precision dose-finding, we assume that individual PK parameters can be either directly measured or reliably predicted. In practice, one can obtain ( V ˆ i , k ˆ i ) from a population PK (popPK) nonlinear mixed-effects model fit to Stage I PK data, combined with baseline covariates (e.g., body size and renal/hepatic function). Further details are provided in Section 4.

In addition, the model parameters will continue to be updated as new data becomes available. Then patient i is assigned to dose d i ∗ given by (2.12) d i ∗ = arg min d | p ˆ i , d − p T | . As a result, each new patient receives the dose that satisfies both (2.12) and (2.9). This continues until the maximum number of patients for the trial is reached. Note that, as this precision dose-finding approach is implemented, selecting a single MTD for the overall trial becomes optional, since each patient is associated with their own MTD d i ∗ .

We conducted simulation studies of the PDF design under different scenarios and compared the results with the CRM design. Following the guidance of Goodman et al. [3], we performed 1,000 trials across five distinct scenarios. Each scenario varied the true coefficient parameters β 0 and β 1 , under the assumption that one of the dose levels is the true MTD. The true coefficient pairs used in these scenarios were { ( − 3 , 1.5 ) , ( − 4 , 1.5 ) , ( − 4.5 , 1.5 ) , ( − 2.5 , 0.6 ) , ( − 1 , 1.2 ) }.

For θ = ( V i , k i , α V , λ V , α k , λ k , β 0 , β 1 , σ 2 ), we assume the prior distributions described in (2.6). In this simulation, V i and k i were assumed to follow gamma distributions, i.e., V i ∼ G a m m a ( 4 , 1 ) , k i ∼ G a m m a ( 3 , 1 ), which is the same as PEDOOP. However, later in Section 3.3.1 we also conducted a sensitivity analysis using different prior distributions.

The target toxicity rate was set to p T = 0.3. Simulated trials assumed five discrete dose levels d ∈ { 15 , 30 , 60 , 90 , 120 } and measured the drug concentration of each patient at time points t j ∈ { 1 , 3 , 5 , 7 , 12 , 24 } hours after injection. A total of 30 patients was enrolled for each trial, with the first 21 receiving doses using the traditional dose-finding algorithm for stage I, and the remaining 9 patients following the precision dose-finding approach for stage II. We set s ∗ = 0.95 in (2.9).

Note that due to the structure of our model, the true toxicity probability at each dose level did not exist. Instead, we generated 2,000 random pairs of V i and k i from their prior distributions and plugged them into (2.4), along with true β 0 and β 1 , to simulate the true average toxicity probability for each dose level, averaged across the 2,000 simulated patients. These simulated averages were reported as “True avg p d ” in the first line of each scenario in Tables 3.1 and 3.2.

The CRM design was based on the model p d = ϕ d exp ( − β ) , d = 1 , … , 5, where the skeleton ϕ = ( ϕ 1 , … , ϕ 5 ) was a set of prior guesses on toxicity probabilities at each dose level, and β followed a normal prior distribution with mean 0 and standard deviation 2. For fair comparison, the CRM design sample size was set to 21.

The simulation results of our approach are presented in Tables 3.1 and 3.2, reporting the results of the first 21 patients and the last 9 patients, respectively. “SC” denotes scenarios from 1 to 5; “True avg p d ” represents the true average toxicity rates at each dose level; “ Y d / n d ” is the empirical toxicity rate at dose d; “ n d ” is the average number of patients assigned to each dose level over 1 , 000 simulated trials; finally, “Sel%” is the proportion of simulated trials in which each dose is selected as the MTD. “No MTD” indicates early termination or failure to recommend a safe MTD.

In scenarios 1 and 2, the true MTDs are dose levels 2 and 3, respectively. The proposed PDF design demonstrates strong performance in both scenarios, achieving the highest Sel% at the correct doses. Moreover, the highest n d values also occur at the correct MTDs, showing efficient patient allocation. In addition, PDF allocates fewer patients and selects MTDs less frequently at overly toxic dose levels.

Both PDF and CRM designs assign the most patients to dose level 3 in scenario 3, where the true MTD is dose level 4. However, the PDF design assigns a slightly higher average number of patients ( n d ) to this dose level compared to the CRM. Additionally, both the PDF design and the CRM correctly identify the true MTD. Notably, our approach assigns fewer patients to, and selects less frequently, the overly toxic dose level 5 compared to CRM.

Scenario 4 presents a challenging case in which dose levels 3 − 5 are all close to the target p T . PDF and CRM perform comparably, allocating patients to mostly dose levels 3 and 4. With larger sample sizes, both designs may eventually recognize dose level 5 as the true MTD. In scenario 5, all doses are overly toxic. PDF assigns slightly fewer patients to toxic doses with much lower incorrect MTD selection.

Table 3.2 presents the simulation results for the last nine patients, who underwent precision dose-finding. Row 2 of each scenario presents Y d ∗ / n d ∗ , where n d ∗ and Y d ∗ are the numbers of patients among the last nine that are assigned to dose level d and that experience DLT at dose level d, respectively. We can see that most Y d ∗ / n d ∗ values are around or smaller than p T = 0.3, regardless of the true average p d values. This shows the benefit of precision dose-finding, since each patient is assigned based on his/her PK-implied optimal dose, regardless of what the expected toxicity rate p d is.

In the test stage, each patient’s dose is determined by his or her own PK profile. Because these PK values are predictions rather than directly observed quantities, they are subject to estimation error. To model this, we consider adding errors to V i and k i and assume V ˆ i = V i + e i , e i ∼ Truncated Normal ( 0 , σ V 2 9 , − V i , ∞ ) k ˆ i = k i + g i , g i ∼ Truncated Normal ( 0 , σ k 2 9 , − k i , ∞ ) where each error term has mean 0, variance equal to 1 / 9 of the true variance, and is truncated below at the negative of its true value to prevent negative PK parameters. When dosing, we compute p ˆ i , d by substituting V ˆ i and k ˆ i into the model. Table 3.3 summarizes the simulation results under this framework. Although toxicity rates increased slightly, the overall pattern remains unchanged, demonstrating that the precision dose-finding stage is robust to realistic PK-prediction errors.

Figure 3.1 displays the dose chosen for the last patient (# 30) in each scenario. In general, the largest bubble (count) increases in SC 1 to 4, whereas SC5 keeps the patient at the lowest dose level. The DLT rates in each scenario are all below p T = 0.3, implying that the precision-based dosing is safe.

We also present outcomes of six randomly selected trials from scenario 1 (Figure 3.2). It displays dose escalation patterns for the first 21 patients and varying dose assignments for the last 9 patients, who participated in the precision dose-finding. Notably, the last 9 patients tend to receive doses centered near the true MTD. The plots further reveal that some patients received the same dose but experienced different DLT outcomes. For instance, patient 22, 23, and 30 in “SC1 - Trial3” had PK profiles ( V i , k i ) of ( 6.53 , 3.18 ) , ( 1.80 , 1.65 ) , and ( 1.51 , 1.34 ). Other trials, not shown here, also suggest that patients with a higher V i / k i ratio may be more likely to experience DLT.

To further illustrate the precision dose-finding, in Figure 3.3 we report the predicted MTD for a patient with different values of V i and k i , given by d i = V i k i exp ( l o g i t ( p i ) − β ˆ 0 β ˆ 1 ), rounded up to the nearest tested dose, where p i = 0.3 and the coefficients are set to their posterior means. It can be seen that as k i increases, the predicted MTD also increases, implying a fast elimination rate helps increase tolerance to toxicity. Similarly, when V i increases, the predicted MTD also increases. For example, given k i = 4, the predicted MTD is dose level 3, 4, and 5 for V i value of 2.95, 4.01, and 6.5, respectively. This is also sensible if the toxicity is inversely related to V i as indicated in model (2.4).