The role of pharmacokinetic (PK) and pharmacodynamic (PD) modeling and simulation is expanding in almost all drug development processes, from non-clinical to different phases of clinical trials.[1] The number of regulatory decisions, including new drug approval and labeling, that were effected by pharmacometric analysis increased from 45 submissions between 2000 and 2004 to 87 submissions between 2007 and 2008.[23] The population approach or pharmacometric analysis based on computational methods is being applied to bridging studies, proof of concept studies for go/no-go decision, simulation for dose selection or study design, and extension to other indications (drug repositioning), among many others.[1] PK/PD modeling and simulation is also used for individualized pharmacotherapy based on relevant demographic factors including race, age, sex, weight, height and genotype.[45]

Pharmacometric analysis usually relies on non-linear mixed effect models to explain and quantify time-varying PK/PD parameters, disease progression, and their relationships based on population data. Pharmacokinetic models describe the change of drug concentration over time and pharmacodynamic models quantify the relationship between concentration and effect. PK/ PD models link PK models to PD models using direct effect models, effect compartment models, or indirect response models.[67] This mathematical method is also used for constructing disease progression models comprising the natural history of a disease and drug effect.[8] E_max models play important roles in describing excitatory or inhibitory exposure-response relationship. Sigmoid E_max models, also known as Hill equations, have successfully explained pharmacodynamics of many drugs with Hill coefficients determining sigmoidicity or sensitivity between exposure and response.[9] However, drug exposure often does not cover the higher range of values required to estimate maximum effect (E_max) since higher concentrations can either result in toxicity or the testing doses in dose escalation study have little information on dose-response surface. This leads to potentially biased or imprecise PD parameter estimates or overly-simplified models such as linear or log-linear models that are not useful for extrapolation.

The purpose of this simulation study is to explore the limitation of the population PK/PD analysis using data from a clinical study and to help to construct an appropriate PK/PD design that enable precise and unbiased estimation of both fixed and random PD parameters in PK/PD analysis under different doses and Hill coefficients.

Eight scenarios out of total seventy scenarios failed stochastic simulation and estimation (SSE), in spite of numerous attempts. The failed scenarios were: four scenarios of 12.5 and 25 mg at high Hill coefficients (γ = 5, 10) in single-dose study and four scenarios of 100, 200, 400 and 800 mg at low Hill coefficients (γ = 0.5) in multiple-dose study. Among the remaining sixty two scenarios which were successful, one scenario in single-dose studies had 1 failure and eight scenarios in multiple-dose studies had 1 to 3 failure(s) during 100 simulations and estimations for each scenario.

This simulation study assumes a typical human pharmacology study design where plasma concentration and effect are frequently measured in order to characterize PK and/or PD. Dense sampling strategies are used in all phases of clinical trials for changing dose, indication, or target population as well as phase 1 clinical trials in new drug development, and pharmacometric analyses are frequently done for these studies. However, not all analyses are successful and the quality of population approach depends on the quality of study design as well as the quantity of data. Therefore, this study aims to specifically investigate the estimation performance of PD parameters based on different data properties given the same quantity of data. In case of sparse sampling, it is natural that parameters are poorly estimated.

Ette et al. suggested a data supplementation (PK/PD knowledge creation or discovery) method to characterize an unknown target region of the response surface, which is very useful in finding appropriate new doses in a completed clinical study with limited information on exposure-response relationship.[13] Dutta et al. found that E_max, EC₅₀, and γ of sigmoid E_max model were poorly estimated if the maximum concentration was not attained up to 95% of E_max applying clinical data simulation. Seven escalating C_max levels with five different Hill coefficients (γ=0.2, 0.5, 1, 2 or 5) were used in simulations. The number of blood samples was ten in each scenario. However, the concentrations were derived from sigmoid E_max curve and pharmacokinetics was not considered.[14] The parameter estimation of sigmoid E_max model is an important part in dose-finding study. Dragalin et al. found that the adaptive designs with the selected dose applying Fisher information matrix were superior in finding target dose compared to the common designs with the doses allocated by equal interval.[15] Wang et al. also compared the performance of the adaptive design to the uniform design.[16] The purpose of these two studies was to show the performance of adaptive optimal designs for dose-finding studies, (corresponding to phase 2 clinical studies) and pharmacokinetics was not considered. On the other hand, the next two studies used pharmacokinetic sampling designs to investigate pharmacodynamic parameter estimation, as is done in this study. Girgis et al.[10] found that the accuracy and precision of the PD parameter estimates got worse as the number of subjects decreased from 100 to 25 with sparse sampling design. Pai et al.[11] performed a series of studies that evaluated the effect of octreotide on the insulin-like growth factor-1and the effect of remifentanil on electroencephalogram, which were successfully described using sigmoid E_max model with moderate sigmoidicity (γ=2.51) and simple E_max model (γ=1), respectively, while the previous study of theophylline-induced eosinopenia was described using a model of high sigmoidicity. The latter study investigated the change in bias and precision in relation to Hill coefficient using 100 subjects in each scenario and 4 to 5 sparse samplings based on the former study and found that the plasma concentration of two to three times EC₅₀ was needed to get reliable parameter estimates if the drug exhibited high Hill coefficient, and that the plasma sampling near EC₅₀ was important with simple E_max model. In the two studies, all the estimated ω_γ² were biased and imprecise, whereas ω_γ² of 10 scenarios, including doses of 50 (C_max=0.42 EC₅₀ units) and 100 mg (C_max=0.85 EC₅₀ units) with γ of 0.5, 50 to 200 mg with γ of 1, 100 and 200 mg with γ of 2 or 5 and 200 mg (C_max=1.7 EC₅₀ units) with γ of 10, among 31 single-dose scenarios were accurately and precisely estimated in this study with dense sampling design.

In this study, the number of subjects was set to 100 in each scenario and dense sampling design was adopted based on the results of the previous studies to investigate the estimation performance with regard to differing concentration and Hill coefficients while fixing the number of subjects and sampling points. Using virtual drugs with different Hill coefficients and ratios of concentration to EC₅₀ would facilitate general interpretation and easy application to the real drugs.

For the single-dose studies, the SSE processes failed when the doses were 12.5 (C_max=0.11 EC₅₀ units) and 50 mg (C_max=0.42 EC₅₀ units) and γ was either 5 or 10, because the simulated concentrations were located at lower flat area of the sigmoid E_max curve (Fig. 2). The accuracy and precision of θ_Emax and ω_Emax² were getting better with the increase in dose, regardless of γ, while those of θ_EC50 and ω_EC50² gradually improved with increase in both dose and γ (Fig. 3). In single-dose studies, since the baseline of effects before dosing was measurable, all the θ_E0 and ω_E0² were estimated with accuracy and precision. The estimated θ_γ showed various trends at different values of γ. The absolute values of relative bias for θ_γ were represented by a U-shaped curve when γ was 0.5, 5 or 10, while they grew with increase in dose when γ was either 1 or 2 (Fig. 3). The relative RMSE for θ_γ and the relative bias and relative RMSE for ω_γ² also represented U-shape trends regardless of γ (Fig. 3), because the bottom area of the effects were measured with lower doses and the ceiling area were measured with higher doses. Sampling enough concentrations corresponding to the whole range of the effect curve must be crucial to estimate reliable θ_γ and ω_γ² (Fig. 2).

The estimation capabilities of the PD parameters except for θ_E0 and ω_E0² with different Hill coefficients (γ) in single-dose scenarios are as follows (Fig. 3). For the γ of 0.5, all the estimates of θ_Emax and ω_Emax² were reliable when the dose was at least 50 mg or C_max was more than 0.42 EC₅₀ units. The estimates of θ_γ and ω_γ² were accurate and precise only under doses of 50 (C_max=0.42 EC₅₀ units) and 100 mg (C_max=0.84 EC₅₀ units). However, θ_EC50 and ω_EC50² were poorly estimated in all dosing scenarios. The reason seemed that the range of concentration-effect relationships was not enough to describe whole shape of sigmoid E_max model as shown in the multiple-dose scenarios. The effect-concentration profiles obtained from various doses seem to be needed to estimate accurate and precise θ_EC50 and ω_EC50². For the γ of 1, all the estimates of θ_Emax and ω_Emax² were reliable when the dose was at least 100 mg or C_max was more than 0.85 EC₅₀ units. θ_EC50 was accurately and precisely estimated when the dose was at least 100 mg (C_max=0.84 EC₅₀ units), while ω_EC50² was hardly estimated, as shown in a previous study.[11] The estimates of θ_γ and ω_γ² were accurate and precise when the doses were within the range of 50 to 400 mg (C_max=0.42 to 3.4 EC₅₀ units). This was not true, however, for 800 mg dose (C_max=6.8 EC₅₀ units) owing to lack of information on sigmoidicity in concentration-effect relationship. For the γ of 2 to 10, all the parameter estimates for E_max and EC₅₀ met the acceptance criteria of bias and precision when the doses were more than or equal to 200 mg (C_max=1.7 EC₅₀ units). This is a relatively smaller value than the result form a study which reported that the 2 to 3 EC₅₀ region is needed to get an unbiased and precise estimate of ω_Emax², when γ was 6.22 and blood sampling was sparse.[10] A dense sampling design is recommended to get a precise ωE_C502, because it was impossible with sparse sampling design.[11]

Both the θ_γ and ω_γ² met the criteria for accuracy and precision at 100 (C_max=0.84 EC₅₀ units) and 200 mg (C_max=1.7 EC₅₀ units) with γ of 2, 100 to 400 mg (C_max=0.84 to 3.4 EC₅₀ units) with γ of 5, and 200 mg (C_max=1.7 EC₅₀ units) with γ of 10. On the other hand, Girgis et al. found that all θ_γ and ω_γ² were biased and imprecise when the Hill coefficient was high (γ=6.22), 4 to 5 blood samples were gathered[10] and Pai et al. also reported similar results from spare sampling design, when γ were 1 or 2.51.[11] These results suggest that dense sampling design may be essential to get accurate and precise estimates of θ_γ and ω_γ².

These results can be applied to clinical trials and PK/PD modeling of resultant data. When a clinical study is performed to evaluate the PK/PD of an investigational product with Hill coefficients of 2, an investigator can use 200 mg (C_max=1.7 EC₅₀ units) to get the best results of parameter estimation or 100 mg (C_max=0.84 EC₅₀ units) to get useful parameter estimates with less toxicity. On the other hand, an investigator can interpret the reliability of the results from PK/PD modeling based on the results of this study.

For the multiple-dose studies (Table 3, Supplementary Fig. 1), the parameter estimation performance was not good. In this simulation setting, the accumulation index was 2.22 with the half-life of 14 hours and 12 hour dosing interval. The concentration-effect profile obtained from this setting could not include the whole range of effects and resulted in poor performance. If the dosing interval was longer than 12 hours or the drug had lower accumulation rate than 2.22, wider range of concentration-time profiles could be gathered. The worse performance of PK parameter estimation in multiple-dose settings than in single-dose settings might play a role in the poor estimates of PD parameter. On the other hand, errors in PD parameter estimates might have negatively affected PK parameter estimation, and this is hinted from different levels of bias and precision with different combinations Hill coefficients and doses, although the Hill coefficient had no relation to pharmacokinetics. If PK/PD modeling is performed with data from multiple-dose studies, the sampled concentrations may include at least 0.63 EC₅₀ and at most 3.9 EC₅₀ corresponding to the doses of 50 to 200 mg in this study.

This study has some limitations. Since the PK profiles were simulated under the assumption of dose proportionality, the relationship between dose and EC₅₀ units might change if investigational drugs have nonlinear PK properties. However, the exposure-response relationship can be easily understood and applied, if concentrations are normalized by EC₅₀. Another limitation is that only one dosing interval was used for simulation in the multiple-dose scenarios. Since a different dosing interval may yield a different steady-state fluctuation, the parameter estimation performance can be better or worse than this study according to the dosing interval. A further limitation stems from the limited range of variability parameters of this study. Since all CV% of IIV and RV were set at moderate values of 30% and 20%, respectively, prediction of exposure-response relationship of a drug with larger variability might show deviations from our findings. Simultaneous modeling of PK/PD is one of the causes of poor IIV estimation, and certain types of sequential approaches may result in better performance.[17]

This simulation study demonstrated the effect of the relative range of sampled concentrations to EC₅₀ and sigmoidicity (or sensitivity) between exposure and response on the parameter estimation performance using dense sampling design. In single-dose studies, most PD parameters of sigmoid E_max model were accurately and precisely estimated when the C_max was more than 0.85 EC₅₀ units, except for typical value and IIV of EC₅₀ which were poorly estimated at low Hill coefficients (γ=0.5 and 1) regardless of doses. The PK/PD modeling results using only steady-state data may be carefully interpreted. This study can be useful in designing a clinical study to evaluate PK/PD relationship for new drug development or drug repositioning.