Construction of Windows for Pharmacokinetic Sampling

This paper describes a method for the construction of pharmacokinetic sampling windows so that they are around the D-optimum time points. Here we consider the situation where a pharmacokinetic (PK) study is accompanied by a dose-finding study in phase I clinical trial. The D-optimal criterion is often used to determine the optimal time for collecting blood samples so that they provide maximum information regarding the population PK parameters. However, collecting blood samples at the D-optimal time points is often difficult. Instead, the sampling time point chosen from a suitable time interval or window can ease the process. The proposed method is conceptually simple and considers the average value and standard deviation of D-optimal time points up to create sampling windows. Random time points can be chosen from these windows then to collect blood samples from the next cohort. The nonlinear random-effects model has been used to model the PK data. Also, we employ the continual reassessment method for dose allocation to the patients. Comparisons of the accuracy and precision for the PK parameter estimates obtained at the D-optimal and random time points are also presented. The results are convincing enough to suggest the proposed method as a useful tool for blood sample collection.


Introduction
Clinical trials are prospective studies to evaluate the effect of interventions in humans under pre-specified conditions. Commonly classified into four phases, clinical trials have become an integral part of drug development. Following pre-clinical and animal studies, phase I trial is the first-in-human experiment, and it focuses on the safety, tolerability, and PK of a drug. Pharmacokinetics is generally defined as what the body does to the drug. It involves studying how drugs enter the body, distribute throughout the body, and leave the body. A drug's pharmacokinetics are determined by the processes of absorption, distribution, metabolism, and excretion. Also, concentrations at the site of action are determined by these processes. After administration of a dose, the PK mechanism transforms it into plasma concentration, and through the systematic circulation of blood, it reaches the site of action and produces an effect.
The theory of optimal designs can guide in PK sampling. More specifically, D-optimum design is often used to determine the optimal time points for blood sample collection. D-optimality

Methods
This section is covering the statistical methodology and models that are used for window construction. Since the plan is to collect blood samples with a dose-finding study's progress, we require a dose-finding design to be used. We consider the continual reassessment method, as described below, to determine the maximum tolerated dose (MTD). We use a nonlinear random-effects model for modeling the PK data, as shown in Section 2.2. The procedure of window construction is discussed afterward.

Continual reassessment method
The continual reassessment method (CRM) is a model-based approach for dose finding in phase I clinical trials that was proposed by O'Quigley, Pepe, and Fisher (1990). The de-sign aims to reduce the number of patients at sub-therapeutic doses and obtain a more accurate estimate of the MTD. For an experimental drug, assume that d ordered doses X = {x (1) , x (2) , . . . , x (d) } are available from pre-clinical studies, and we want to determine the MTD to be used in the next phase. Unlike the original article where the authors utilised a one-parameter logistic model, we use a two-parameter model as suggested by Whitehead and Brunier (1995) and Whitehead and Williamson (1998). That is, the dose-response model to be utilised is ψ(x, θ) = exp(θ 1 + θ 2 x) 1 + exp(θ 1 + θ 2 x) , where θ = (θ 1 , θ 2 ) is the vector of dose-response parameters and x is the dose given to a cohort of size c. If we are at the kth stage of a trial, then we have k cohorts treated so far with different doses from X . Let x be a k × 1 dose vector with components x l and r be a k × 1 outcome vector with r l as the l th row (l = 1, . . . , k) representing the toxic outcomes obtained from a cohort. Then the likelihood function at the kth stage is Since the maximum likelihood estimation is not possible until sufficient information is available, the Bayesian approach is employed to estimate the parameters. The posterior means of θ = (θ 1 , θ 2 ) at the kth stage are obtained aŝ where Θ is the parameter space and g(θ) is the prior distribution of the parameters. The choice of u 1 < θ 1 < u 2 and u 3 < θ 2 < u 4 gives a restricted parameter space asΘ = {θ : u 1 < θ 1 < u 2 , u 3 < θ 2 < u 4 } so that The probability of toxicity at the end of stage k is then updated at each dose aŝ That dose is chosen for the next cohort for which the absolute difference between the estimated probability of toxicity and the target toxicity rate γ is minimum. That is, A trial is continued until the fixed number of cohorts m is achieved, and the MTD is the dose that would be allocated to the cohort m + 1 if that were in the trial. Note that the choice of a value for γ in a trial solely depends on the likely nature of the investigational drug.

PK model
The management of a drug in the body is very complex, as processes like absorption, distribution, metabolism and elimination work to alter drug concentration in tissues and fluids.
To simplify this complexities, mathematical principles are helpful to the processes. The most basic mathematical description of drug distribution and elimination is the one-compartment model. It utilises a single central compartment and assumes that equilibrium is instantaneously achieved with tissues. The one-compartment PK model with first-order absorption of a drug is defined as . . , c and l = 1, . . . , n, where y il is the concentration of a drug in the blood for the ith individual observed at time t il , x is the dose received by the individual and θ i = β + b i is the vector of parameters with β = (V, k a , k e ) T and b That is, an additive model for random effects is considered. The il denotes the random error and n is the number of measurements taken on an individual. We assume that i ∼ N n (0, σ 2 I n ) and b i ∼ N 3 (0, Ω), where Ω = diag(σ 2 1 , σ 2 2 , σ 2 3 ). The vector of all the population parameters to be estimated is represented by Ψ = (β T , λ T ) T , where λ = (σ 2 1 , σ 2 2 , σ 2 3 , σ 2 ) T is the vector of variances. The sampling windows that we construct require D-optimal time points. Since obtaining Doptimal time points requires the Fisher information matrix (FIM), we need the above model's likelihood function. As the model is nonlinear in parameters and parameters are random, there is no closed-form for the likelihood function. Therefore, we linearise the model using a first-order Taylor series expansion to obtain the FIM. Details of the derivation of FIM are available in Appendix A. The response for ith patient is measured at the times, represented by ξ i = (t i1 , . . . , t in ). Also, we describe the population design by Ξ = {ξ 1 , . . . , ξ c }. The population FIM is defined as the sum of c elementary Fisher information matrices, that is, The derivation of M i (Ψ, ξ i ) is shown in Appendix A. For a single group of c individuals with identical designs, the population FIM is Unlike the linear models, the FIM here depends on the model parameters. The uncertainty in a set of parameter estimators can be expressed in terms of the volume of a confidence ellipsoid. The precision of the estimators increases as the volume decreases. The D-optimality criterion minimises the volume of the confidence ellipsoid, which is a function of the determinant of the covariance matrix (Atkinson et al. 2007). More specifically, a design ξ * D is called D-optimal if it minimises Φ D {M (Ψ, Ξ)} = log |M −1 (Ψ, Ξ)|. That is,

Sampling windows
A window is an interval for blood sample collection. The main intention is to take blood samples at some convenient time points from intervals instead of the D-optimal ones. Let t a1 , . . . t an be the D-optimal time points obtained at stage a of a trial (a = 1, . . . , m). At the kth stage, lett .1 , . . . ,t .n be the means obtained so thatt .b = k a=1 t ab /k, for b = 1, . . . , n. Also, let SD(t .1 ), . . . , SD(t .n ) be the standard deviations such that At the kth stage of a trial, a sampling window is constructed ast .b ±δ SD(t .b ), for b = 1, . . . , n.
Here δ can take positive values such as 1, 2, 3, etc. Since standard deviation cannot be computed unless two observations are available, we start constructing sampling windows from the third stage. That is, k = 3, . . . , m, when we construct a window. Following the construction of sampling windows, a uniform distribution is used to select the random time points from those windows. The uniform distribution is used to ensure that any time point from a window can get selected to collect the blood samples. That is, at the kth stage, for the bth time point, we assume that t kb ∼ U (c 1 , c 2 ), where c 1 and c 2 are the lower and upper limits of the constructed sampling window for bth time point. After the random time points are obtained from the windows at each stage of a trial, we collect blood samples from a cohort to measure the concentrations.

Simulation study
To investigate the performance of proposed window construction, we conduct a simulation study. The details of that study and numerical findings are described below.

Simulation setup
We use theophylline data (Beal and Sheiner 1992) for the simulation study. Following the oral doses of theophylline given to twelve subjects, serum concentrations were measured at 11-time points over the next 25 hours. Note that these data are available in R as Theoph. The parameter estimates of the model in (7), obtained for these data, are used for data generation in the simulation study. The dose range found in the theophylline data has been followed in our dose-response scenarios. Since theophylline data do not contain any information on the dose-response relationship, we have assumed four plausible relationships, as shown in Figure  1. In each scenario, we have the discrete doses as X = {3, 3.5, 4, 4.5, · · · , 6}. The scenarios differ only in terms of the steepness in the toxicity curve. The steepness is very high in Scenario 1, while very low in Scenario 4. Here we assume the target toxicity rate γ to be 0.33. The true MTDs in the successive scenarios are 3.5, 4, 5, and 6. Note that these MTDs are the doses at which the probability of toxicity is closest to the target toxicity rate 0.33.  The mean PK parameter estimates obtained from the model fitting are V = 0.020 L, K e = 0.080 and K a = 1.444. The variances associated with these successive parameters are σ 2 1 = 1.28E − 05, σ 2 2 = 5.46E − 05, σ 2 3 = 2.216E − 02, and the error variance is σ 2 = 1.43. This single PK profile is used against each of the four dose-response scenarios in the simulation study. Blood samples are collected at three-time points per patient; that is, the number of patient observations is n=3. It has been found that a 4-point design is more efficient than a 3-point design. Compared to a 4-point design, the gain in a 5-point design is not that substantial. Also, if a 5-point design is considered, some of the sampling time points are very similar. We choose a 3-point design as some of the optimal time points are close to those for a 4-point design. Each trial starts with the lowest dose assigned to a cohort of c=3 patients. A total of 15 cohorts in each trial is assumed, that is, m = 15. However, to assess the impact of sample size on parameter estimates, we also investigate m = 25.
Following the lowest dose given to a cohort, for each individual in the cohort, we generate b i from N 3 (0, Ω) to obtain the individual parameters as θ i = β + b i . Then we find the concentrations at the D-optimal time points. Since the PK model is nonlinear, the FIM depends on both design variable t and the model parameters. Therefore, to find the Doptimal time points, we need to plug in values for the parameters. Starting with some initial values for the parameters, we replace them by the updated parameter estimates at every stage. Note that the D-optimal time points are obtained by using PFIM3.2 (Bazzoli, Retout, and Mentré 2010). The random errors are then generated from N 4 (0, σ 2 I 4 ) to be added to the previously generated concentrations to obtain simulated PK responses for an individual. The same mechanism is used to generate responses for other individuals in the cohort. Since we want to explore the effect of time points from sampling windows on parameter estimation, we need the generation of PK responses for the same individual at those time points. Therefore, we use the same θ i and i to obtain the PK responses on different occasions for an individual. The dose-response outcomes are generated using the binomial distribution.
Below are the steps that we follow at each stage of a trial. As mentioned earlier, k represents the stage of a trial.
Step 1: Administer the best dose x k to the kth cohort. Initially, k = 1 and the best dose is the lowest dose.
Step 2: Determine the D-optimal time points for dose x k . Measure the drug concentration for the cohort at these time points if 1 ≤ k ≤ 2.
Step 3: Construct the sampling windows as explained in Section 2.3, and measure the concentration for the cohort at random time points taken from the windows if k ≥ 3.
Step 4: Observe whether the outcomes are toxic or not for the cohort. Obtain the dose-response parameter estimatesθ k following (3). Also, obtain the PK model parameter estimatesΨ.
Step 5: Determine the dose for the next cohort following (6).
Step 6: Set k = k + 1 and repeat the steps 1-6 until the trial reaches a total of m cohorts.
Step 7: If k = m, determine the MTD for further investigation in the next phase.
At each stage of a trial, we estimate both the PK and dose-response parameters. The PK parameters are estimated using the R package nlme. The dose-response parameters are estimated using the Bayesian approach. We consider a bivariate uniform prior distribution for the dose-response parameters θ. A single parameter space is chosen asΘ = {θ : −9 < θ 1 < −5, 0 < θ 2 < 3}. The dose-optimisation criterion in (6) is used to determine the most appropriate dose for the next cohort. Now for the new cohort, we generate PK responses and dose-response outcomes. Based on the first two cohorts' response, we update both the PK and dose-response parameter estimates and select the dose for the third cohort. Since sampling windows are utilised from the third stage onwards, the generated data are the same for an individual in the first two cohorts in all the data sets. We continue the process of dose selection for every next cohort until the trial reaches m. We have considered two, three, four, and six standard deviations for each scenario to construct the sampling windows, that is, δ = 2, 3, 4, 6. Consequently, we have five PK data sets for model fitting. One of these is for D-optimal time points, and the other four are for time points from the windows. Also, we collect samples at three random time points taken from the sampling region (0, 25) hours. We do not use δ = 1, as we have found it to provide too narrow windows for the time points. A large δ will provide a wider sampling window and consequently will provide flexibility in choosing a time point from it. Although we have found it not to be common, the lower limit of a sampling window may become negative, and in such a case, we assume it as 0. For a large δ, the sampling windows may overlap. As a consequence, the random time point from the first window may be greater than that of the second window. Similarly, the second time point may be greater than that from the third point. Since our target is to collect blood samples at distinct time points, it does not matter from which window a time point comes. Getting different time points is okay with the procedure. All the computations in this paper are conducted using a self-written code in R (R Core Team 2020).

Simulation findings
Now we present the numerical results obtained from simulation. We have conducted 1000 simulations for each of the scenarios to investigate windows' performance on parameter estimation. Table 1 gives the distribution of MTD selection and dose allocation for the assumed dose-response scenarios. Note that bold values in the table correspond to the true MTDs for the respective scenario. When we consider δ = 2 to get the sampling windows and m = 15, dose 3.5 is selected in 83% of the trials as MTD in Scenario 1. This dose is allocated to 55.44% of the cohorts during the trials. Dose 4 is the true MTD in Scenario 2. This dose is selected in 77.20% trials as the MTD and allocated to 54.62% of the cohorts. In Scenario 3, dose 5 is selected in 64.30% trials as the MTD and allocated to 41.18% of the cohorts during the trials. The correct MTD in Scenario 4 is 6, which is recommended in 79.10% trials. As m increases to 25, the accuracy of MTD identification increases for the scenarios. Also, more cohorts receive the actual MTD during a trial. Since we do not incorporate PK information in the dose-finding study, the MTD selection and dose allocation for the scenarios are found not to be affected by δ. Therefore, Table 1 presents results only for δ = 2. Along with the D-optimal time points, Table 6 provides the sampling windows and random time points taken from those windows at successive stages of a trial for Scenario 1. There are no sampling windows for the first two cohorts, as we start constructing windows from the third cohort. Note that blood samples are collected at the D-optimal time points for the first two cohorts when δ = 2, 3, 4, or 6. It is seen that larger values of δ give more expansive windows at the beginning. However, differences in the widths of respective intervals across δ decrease towards the ending cohorts.
The relative bias and coefficient of variation (CV) of PK model parameter estimates, obtained for different choices of δ, are available in Tables 2-5 along with those obtained for D-optimal and random time points. It would be useful to mention that the absolute bias is expressed as a percentage of a parameter's true value to obtain the relative bias. Similarly, the CV for a parameter expresses the standard deviation of estimates as a percentage of the mean estimates. Table 2 indicates that the estimates are close to the true values, as reflected through the relative biases for the mean PK parameters if m = 15. The relative bias for variance components is relatively higher than those for the mean parameters in Scenario 1. In most of the occasions, a high relative bias is found for the error variance. Generally speaking, the CV for a parameter increases as δ increases. The CV is high for σ 2 3 when samples are collected at random time points. There is not much difference in the estimates if we utilise time points randomly chosen from the windows, instead of the D-optimal time points. If we move to m = 25, then the CV for parameters decrease.     Tables 3-5 present the relative bias and CV of parameter estimates for Scenarios 2-4, respectively. The results in these cases are similar to those obtained earlier for Scenario 1. As δ increases, we do observe an increase in relative bias for some of the parameters. The CV increases if we increase δ. That is to say, accuracy and precision in the parameter estimates are affected to some extent if we consider proposed windows to collect blood samples instead of D-optimal time points. However, the loss in accuracy and precision is very negligible compared to the flexibility that windows open.

Conclusion
This paper presents a simple way of constructing windows for PK sampling in a phase I trial. Blood samples can be collected at time points chosen from windows when patients or volunteers have complexity and hardship for giving samples to medical personnel at the exact D-optimal time points. We compared the PK parameter estimates at D-optimum and randomly chosen time points from windows. It has been found that the estimates differ over the two approaches very little. Also, we have investigated three random time points taken from the design region. The estimates of variance component σ 2 3 are found to be affected often in such a case. Since the design region in our example is 0 to 25 hours, the random time points from sampling windows and the random time points from the design region are very similar. Therefore, the random time points case produces very competitive estimates. However, in case of a much broader design region than this example, any time point taken from the region is likely to produce less accurate and efficient estimates. Four dose-response scenarios have been studied for various values of δ and m. A large value of δ provides more flexibility in sampling but provides less accurate estimates of parameters. On the other hand, with the increase in m, the accuracy of both MTD recommendations and PK parameter estimates increases. Although the CRM has been used for dose-escalation, one can also use other dose-finding designs. Other PK models can be used to construct windows as well. Note that we do not use PK data to determine the MTD. The PK parameter estimates obtained at each stage of a trial are used to find the locally D-optimum time points. To conclude, we have shown how to construct windows to choose time points around the D-optimum time points. The comparison of estimates for each of the scenarios ensures the accuracy and precision of the proposed method.
Appendix A: Fisher information matrix for the PK model As in Alam (2015), we linearise the model in (7) using a first-order Taylor series expansion about φ i = (V, k e , k a , b V i , b ke i , b ka i ) T at φ 0 = (V 0 , k 0 e , k 0 a , 0, 0, 0) T , where V 0 , k 0 e and k 0 a are some prior values of the population mean parameters. Then we obtain an approximated linear mixed effects model of the form where µ i is a n × 1 vector of constants. The matrix H i is defined as The elements of H i are given as Since θ i = β + b i , we also have H i = L i , where It can be shown that the approximate FIM for individual i is where (W i ) lk for l, k = 1, 2, 3 are elements of W i = L T i V −1 i L i and (P i ) ll are diagonal elements of P i = L T i V −2 i L i . Note that V i = Var(y i ) ∼ = L i ΩL T i + σ 2 I n .