Commentary: Determination of the 95% effective dose of remimazolam tosylate in anesthesia induction inhibits endotracheal intubation response in senile patients

A Commentary on Determination of the 95% effective dose of remimazolam tosylate in anesthesia induction inhibits endotracheal intubation response in senile patients by Liangchao Qu, Mei Liu, Ru Ouyang, Tianyuan Li, Dingde Long, Yao Jiang, Chengyu Wang and Liqin Cheng (2023). Front Pharmacol. 14:1136003. doi: 10.3389/fphar.2023.1136003 Introduction Medication use in special populations—such as older adults, children, frail patients, and obese individuals—often presents clinical challenges due to insufficient evidence-based guidance. Consequently, research targeting these groups holds significant value for informing pharmacotherapeutic decisions. This imperative extends to emergency surgical patients, whose distinct physiological responses may yield clinically actionable evidence not generalizable from routine populations(Pouraghaei et al. 2013). A review revealed that preoperative frailty is significantly correlated with adverse perioperative outcomes(Kikuchi et al. 2021). Consequently, it is highly important for anesthesiologists to pay more attention to frail patients. In frail patients undergoing general anesthesia, the dose of an anesthetic drug necessary to suppress the response to tracheal intubation might be lower than that in non-frail patients. Therefore, applying the recommended dose for non-frail patients to frail patients might lead to harm to patients, such as hemodynamic fluctuations. Remimazolam, a novel drug, is recommended for anesthesia induction because of its lesser impact on circulation(Hu et al. 2022). Studies have demonstrated that the 50% effective dose (ED50) and 95% effective dose (ED95) of remimazolam for sedation during gastroscopy in elderly patients are 0.153 mg/kg and 0.162 mg/kg(Sun et al. 2022), respectively . This could offer a dose recommendation for remimazolam to suppress the response to tracheal intubation in non-frail elderly patients. However, the ED95 obtained through probit regression based on the ED50 is not reliable, and researchers have proposed biased coin design (BCD). Qu et al. used a BCD and isotonic regression to determine that the 95% effective dose of remimazolam used for anesthesia induction inhibits the endotracheal intubation response in senile patients(Qu et al. 2023). This is an interesting study, and the results of the precise ED95 are sufficiently attractive for clinical anesthesiologists. Nevertheless, we observe that the results presented in the article are perplexing and that the analysis of the data is incomplete. The data of adjusted response rates were not reproducible. Qu et al. assert that isotonic regression is employed to calculate ED95, and the pool-adjacent-violators algorithm (PAVA) is employed to adjust response rates. PAVA replaces the data of adjacent violations with a weighted mean to guarantee the nondecreasing nature of the data, which can be acquired by utilizing the gpava function within the isotone package. However, the adjusted response rates of non-fail patients were not in accordance with the PAVA results (Table 1). Particularly in the Frail group, the response rates and doses were monotonically increasing. Thus, the results of PAVA could not alter the data. In reality, simple rules based on PAVA do not have to utilize a computer to obtain an answer. For the aforementioned reasons, the outcomes calculated by authors are highly astonishing. Point and interval estimations demand further assessment. A linear interpolation was employed to estimate the dose with a target probability of Γ = 0.95, which essentially constitutes a straight line fit with explicit formulas for calculation(Pace et al. 2007). Hence, the ED95 calculated by the Qu et al. on the basis of dubious PAVA-adjusted response rates is also subject to skepticism. The authors assert that the 95% confidence interval (CI) for ED95 was calculated via a bias-corrected percentile method(Stylianou M et al. 2003). Software has been written in GAUSS and R(Pace et al. 2007), but the code seems not to be publicly accessible, and researchers have to program it on their own. Discrepancies in code may result in marginally different outcomes. In this algorithm, the interval estimate of the dose is computed on the basis of the inverse function of the empirical cumulative distribution; therefore, the range does not exceed that of the trial doses (Table 2). An R package (https://rdrr.io/cran/ed50/) for ED50 provides a similar calculation, and it is also explicitly indicated that the CIs are bounded as the bootstrap parametric resample function generates additional doses strictly within the dose sample space. We observe that the results of Qu et al. extend far beyond the range of trial doses, a situation that does not seem to align with the bias-corrected percentile method. An algorithm known as centered isotonic regression improves the PAVA by achieving the monotonically increasing of the adjusted response rates and doses(Oron and Flournoy 2017). The CI with this algorithm would lie outside the range of the trial dose, as it offers an algorithm distinct from the bias-corrected percentile method. Consequently, the authors are obligated to verify the correctness of the code. In addition, Figure 3B provided by the author is misplaced because it is a copy of Figure 3A. Detailed reporting is needed for the statistical analyses of secondary clinical outcomes. Qu et al. reported that DBP was significantly greater in the non-frail group than in the frail group at 1 and 2 min postintubation (p = 0.006 and 0.010, respectively). As depicted in their statistical analysis, two-factor repeated measures analysis of variance (ANOVA) was used to compare the differences between the two groups. The authors fail to report whether the assumptions employed for repeated-measures ANOVA are valid, such as normality and homogeneity of variance across groups and levels. More crucially, the authors fail to report relevant statistics or effect sizes. If normality and homogeneity of variance are not fulfilled, it is suggested that more robust statistical methods, such as generalized estimating equations, be employed. In addition, from Figure 4B in the paper, the baseline heart rate (HR) of the two groups was significantly different (approximately 15 bpm), and it is not evident whether the authors corrected it as a covariate. Qu et al. compared the demographic characteristics of the Non-frail group and the frail group to demonstrate the balance of the baseline data of the two cohorts. This is indispensable for the comparison of ED95 between the two cohorts. However, it appears that the authors intend to explore the influences of remimazolam on BP and HR in both cohorts. There is no evidence suggesting that the drug dose has no effect on BP or HR; rather, it is more probable that the drug dose is a confounding factor influencing BP and HR in different cohorts. The proportions of drug doses utilized in the two cohorts were dissimilar (Figure 1), so we propose that the dose should also be incorporated as a covariate in the analysis to clarify whether the variations in BP and HR are attributed to exposure factors (frail) or confounding factors (drug). Given that the authors failed to report detailed information regarding the statistical analyses, the reliability of their results cannot be trusted. Conclusion Dose finding studies typically utilize standardized methodological frameworks and often generate expected results—characteristics that could foster excessive confidence in their conclusions. Although computational discrepancies are relatively rare in dosing trials, they remain a potential concern. We emphasize the importance of critical appraisal when interpreting research outcomes; clinical application of conclusions based on biased evidence poses substantial risks to patient safety. Prior to translating such evidence into practice, thorough validation of the reliability of study conclusions is essential, particularly for findings lacking comprehensive evaluation. The study conducted by Qu et al., which employed BCD to determine the ED95 of remimazolam in senile patients, is of interest. Nevertheless, the incomplete reporting of data processing has led to doubts regarding the accuracy of the results. As readers, we look forward to more detailed data reports or errata from the authors. Table 1 Adjusted response rates of different algorithms. Dose Observed reponse rates Calculated by Qu et al. gpava# cirPAVA## Non-frail Group 0.30 0.929 0.870 0.828 0.828 0.32 0.727 0.920 0.828 0.885 0.34 Frail Group 0.26 0.667 0.630 0.667 0.667 0.28 0.875 0.940 0.875 0.875 0.30 Note: # gpava is a function in the R package of 'isotone'; ## cirPAVA is a function in the R package of 'cir'. Table 2 Calculation results of ED95 and 95% CI for different codes. Note: * The codes are shown in Appendix I. Group Range of dose ED95 and 95%CI Calculated by Qu et al. Calculated by us* Calculated with CIR Non-frail 0.30-0.34 0.331 (0.272-0.472) 0.334 (0.294-0.340) 0.331 (0.318-0.456) Frail 0.26-0.30 0.297 (0.231-0.451) (0.263-0.300) 0.292 (0.287-0.406) Figure 1 Comparison of the dose distributions between the two groups. 1. options(digits = 10) 2. data0<-read.csv("name.csv",header = T) 3. IR<-function(data0){ 4. freq0<-as.data.frame(table(data0[,1])) 5. positive<-c(1:nrow(freq0)) 6. p<-c(1:nrow(freq0)) 7. freq1<-cbind(freq0,positive,p) 8. 9. for (i in 1:nrow(freq1)) { 10. freq1$positive[i]<-sum(data0[,1]==freq1$Var1[i]&data0[,2]==1) 11. } 12. freq1$p<-freq1$positive/freq1$Freq 13. c0<-as.numeric(as.character(freq1$Var1)) 14. p0<-as.numeric(as.character(freq1$p)) 15. library(isotone) 16. 17. p1<-gpava(y=p0,z=c0) 18. } 19. 20. p1<-IR(data0) 21. c0<-unlist(p1[2]) 22. p0<-unlist(p1[3]) 23. ap<-unlist(p1[1]) 24. 25. finddose<-function(Γ,z,x){ 26. 27. if(Γ max(x)){ 33. km<-length(x) 34. p_max<-x[km]-x[km-1] 35. if(p_max==0){μ<-Γ/x[km]*z[km]}else{ 36. μ<-(Γ-x[km-1])/p_max*(z[km]-z[km-1])+z[km-1]} 37. 38. }else if(Γ %in% x){ 39. sameX<- match(Γ,x) 40. μ<-mean(z[sameX]) 41. 42. }else{ 43. k<-findInterval(Γ,x) 44. dp<-x[k+1]-x[k] 45. if(dp==0){μ3<-0.5*(z[k]+z[k+1])}else{ 46. μ<-(Γ-x[k])/dp*(z[k+1]-z[k])+z[k]} 47. }} 49. u3<-finddose(0.95,c0,ap) 50. names(u3)<-"thet" 51. 52. boots<-function(Γ,data){ 53. df0<-c(1:2000) 54. for(i in 1:2000){ 55. 56. dataS<-function(data){ 57. rowS<- sample(1:nrow(data),nrow(data),replace = TRUE) 58. data[rowS,]} 59. 60. datai<-dataS(data) 61. p1_i<-IR(datai) 62. df0[i]<-finddose(Γ,unlist(p1_i[2]),unlist(p1_i[1])) 63. } 64. 65. mean(df0) 66. } 68. u_mean<-boots(0.95,data0) 69. 70. findp<-function(dose,z,x){ 71. h<-findInterval(dose,z) 72. if(x[h+1]==x[h]){p<-x[h]}else{ 73. p<-(dose-z[h])/(z[h+1]-z[h])*(x[h+1]-x[h])+x[h]} 74. } 76. p_mean<-findp(u_mean,c0,ap) 77. names(p_mean)<-"Probability" 78. 79. doseCI<-function(p){ 80. z0<-qnorm(p) 81. za<-qnorm(0.975) 82. a1<-pnorm(2*z0-za) 83. a2<-pnorm(2*z0+za) 84. u1<-finddose(a1,c0,ap) 85. u2<-finddose(a2,c0,ap) 86. c(u1,u2)} 87. 88. CI_95<-doseCI(p_mean) 89. names(CI_95)<-c("CI_L","CI_U") 90. 91. list(data0,c0,p0,ap,u3,CI_95) Appendix Ⅱ A Brief Introduction to Statistical Methods Probit regression Probit regression is a statistical method used to model the relationship between one or more predictor variables and a binary outcome (e.g., success/failure). In dose-finding studies, the model is employed to fit the cumulative distribution function of drug response. This allows for the estimation of doses corresponding to specific response probabilities, such as ED₅₀ and ED₉₀. It is important to note that this model rely on the key assumption that the probability density function follows a normal distribution. Biased coin design Biased coin design is a response-adaptive randomization technique used primarily in clinical trials to balance treatment group sizes while incorporating accumulating outcome data to preferentially assign patients to the seemingly superior treatment. In dose-finding studies targeting metrics such as ED90, biased coin design concentrates experimental doses near the target probability threshold (e.g., 0.9), thereby enhancing the precision of parameter estimation. Isotonic regression Isotonic regression is a nonparametric statistical technique that fits a monotonic function to observed data while minimizing the sum of squared errors. Unlike parametric models (e.g., probit regression), it imposes no predefined functional form but strictly enforces that the estimated response function is either entirely non-decreasing or non-increasing. In dose-finding studies, isotonic regression corrects non-monotonic dose-response patterns caused by random variability, ensuring response rates adhere to biological plausibility. The adjusted response rates are subsequently used to fit the empirical cumulative distribution function. Pool-adjacent-violators algorithm Pool-adjacent-violators algorithm is the fundamental computational procedure for solving the isotonic regression problem. It efficiently finds the best monotonic least-squares fit to a sequence of observed data points while preserving their original order. Bias-corrected percentile method Bias-corrected percentile method is an enhanced bootstrap technique for constructing confidence intervals that adjusts for both bias and skewness in the sampling distribution, offering improved accuracy over the standard percentile method. In dose-finding studies, isotonic regression typically derives point estimates (e.g., ED90) via linear interpolation, while the bias-corrected percentile method computes interval estimates; bootstrap resampling implement the robustness of interval estimates. Centered isotonic regression Centered isotonic regression is a modified isotonic regression that estimates a monotonic function centered around a target parameter, minimizing the sum of squared errors in transformed space. In dose-finding studies, the non-decreasing property of standard isotonic regression can create mathematical challenges when inverting the function to derive dose estimates. Centered isotonic regression addresses this limitation by generating a strictly increasing function, which better aligns with the biologically imperative monotonicity of drug response cumulative distribution functions. This approach simultaneously circumvents mathematical complexities inherent in dose inversion procedures. Repeated measures analysis of variance Repeated Measures analysis of variance analyzes within-subject changes across ≥ 3 timepoints or conditions while accounting for intra-individual correlations. Critical assumptions: A. Normality: Residuals at each timepoint should be normally distributed. B. Sphericity: Variances of differences between all pairs of repeated measurements must be equal. C. Tested via Mauchly's test. D. Homoscedasticity: Equal error variances across experimental groups. E. Other assumptions. Generalized estimating equation Generalized estimating equation is a semi-parametric method for analyzing correlated data (e.g., longitudinal/repeated measures). This model demonstrates considerable robustness in statistical inference.