library(shiny) library(bslib) library(mrgsolve) library(dplyr) library(ggplot2) aggregate_model_code <- ' $PARAM @annotated KM : 1.92 : Michaelis-Menten constant (mg/L) VMAX1 : 110 : Maximum metabolic rate at 1 h (mg/h) CLL : 3.60 : Linear clearance (L/h) EXPONENT : 0.31 : Exponent of lean body weight on CL and Vmax V2 : 59 : Central volume (L) V3 : 79 : Peripheral volume (L) Q : 27 : Intercompartmental clearance (L/h) KA : 1.0 : Absorption rate constant (h^-1) LOGITF : 1.32 : Logit oral bioavailability parameter T50 : 2.41 : Time to 50% inhibition of Vmax (h) VMAXINH : 0.82 : Maximum proportion reduction of Vmax PM_VMAX : 0.59 : Proportion of Vmax for poor metabolizers EINH_CL : 0.50 : Enzyme inhibitor effect on clearance EIND_CL : 2.00 : Enzyme inducer effect on clearance WT : 67 : Total body weight (kg) LBW : 54 : Lean body weight (kg) NSEX : 0 : Sex flag (0 male, 1 female) PM : 0 : CYP2C19 poor metabolizer flag EINH : 0 : Enzyme inhibitor coadministration flag EIND : 0 : Enzyme inducer coadministration flag $CMT @annotated DEPOT : Oral depot compartment (mg) CENT : Central compartment (mg) PERIPH : Peripheral compartment (mg) $MAIN double main_f1 = exp(LOGITF) / (1.0 + exp(LOGITF)); F_DEPOT = main_f1; $ODE double ode_sex_ref = NSEX > 0.5 ? 42.0 : 54.0; double ode_vmax1 = VMAX1 * pow(PM_VMAX, PM) * pow(LBW / ode_sex_ref, EXPONENT); double ode_cll = CLL * pow(LBW / ode_sex_ref, EXPONENT); double ode_v2 = V2 * pow(WT / 70.0, 1.0); double ode_v3 = V3 * pow(WT / 70.0, 1.0); double ode_q = Q * pow(WT / 70.0, 1.0); double ode_tadj = SOLVERTIME < 1.0 ? 1.0 : SOLVERTIME; double ode_vmax_now = ode_vmax1 * (1.0 - VMAXINH * ((ode_tadj - 1.0) / ((ode_tadj - 1.0) + (T50 - 1.0)))); double ode_enz = pow(EINH_CL, EINH) * pow(EIND_CL, EIND); double ode_cp = CENT / ode_v2; double ode_clt = (ode_cll + ode_vmax_now / (KM + ode_cp)) * ode_enz; dxdt_DEPOT = -KA * DEPOT; dxdt_CENT = KA * DEPOT - (ode_clt / ode_v2) * CENT - (ode_q / ode_v2) * CENT + (ode_q / ode_v3) * PERIPH; dxdt_PERIPH = (ode_q / ode_v2) * CENT - (ode_q / ode_v3) * PERIPH; $TABLE double tab_sex_ref = NSEX > 0.5 ? 42.0 : 54.0; double tab_vmax1 = VMAX1 * pow(PM_VMAX, PM) * pow(LBW / tab_sex_ref, EXPONENT); double tab_cll = CLL * pow(LBW / tab_sex_ref, EXPONENT); double tab_v2 = V2 * pow(WT / 70.0, 1.0); double tab_tadj = TIME < 1.0 ? 1.0 : TIME; double tab_vmax_now = tab_vmax1 * (1.0 - VMAXINH * ((tab_tadj - 1.0) / ((tab_tadj - 1.0) + (T50 - 1.0)))); double tab_enz = pow(EINH_CL, EINH) * pow(EIND_CL, EIND); double CP = CENT / tab_v2; double CLLIN = tab_cll * tab_enz; double CLNL = (tab_vmax_now / (KM + CP)) * tab_enz; double CLTOT = CLLIN + CLNL; double F_ORAL = exp(LOGITF) / (1.0 + exp(LOGITF)); $CAPTURE CP CLTOT CLLIN CLNL F_ORAL ' m1_model_code <- ' $PARAM @annotated CL : 5.2 : Linear clearance (L/h) V : 92 : Central volume (L) KA : 1.1 : Absorption rate constant (h^-1) F1 : 0.63 : Oral bioavailability (fraction) RIF : 0 : Rifampicin coadministration flag CHOL : 0 : Grade 3 cholestasis flag $CMT @annotated DEPOT : Oral depot compartment (mg) CENT : Central compartment (mg) $MAIN F_DEPOT = F1; $ODE double ode_cl = CL * (1.0 + 3.0 * RIF) * (1.0 - 0.52 * CHOL); double ode_v = V; dxdt_DEPOT = -KA * DEPOT; dxdt_CENT = KA * DEPOT - (ode_cl / ode_v) * CENT; $TABLE double CLi = CL * (1.0 + 3.0 * RIF) * (1.0 - 0.52 * CHOL); double CP = CENT / V; double F_ORAL = F1; $CAPTURE CP CLi F_ORAL ' m2_model_code <- ' $PARAM @annotated CL : 6.95 : Linear clearance (L/h) V : 200 : Central volume (L) KA : 1.1 : Absorption rate constant (h^-1) F1 : 0.895 : Oral bioavailability (fraction) AGE : 61 : Age (years) PM : 0 : CYP2C19 poor metabolizer flag ALP : 104 : Alkaline phosphatase (U/L) AGE_CL : 0.012 : Age effect on clearance PM_CL : 0.37 : Poor metabolizer effect on clearance ALP_CL : 0.0016 : ALP effect on clearance $CMT @annotated DEPOT : Oral depot compartment (mg) CENT : Central compartment (mg) $MAIN F_DEPOT = F1; $ODE double ode_cl = CL * (1.0 - AGE_CL * (AGE - 61.0)) * (1.0 - PM_CL * PM) * (1.0 - ALP_CL * (ALP - 104.0)); if(ode_cl < 0.2) ode_cl = 0.2; double ode_v = V; dxdt_DEPOT = -KA * DEPOT; dxdt_CENT = KA * DEPOT - (ode_cl / ode_v) * CENT; $TABLE double CLi = CL * (1.0 - AGE_CL * (AGE - 61.0)) * (1.0 - PM_CL * PM) * (1.0 - ALP_CL * (ALP - 104.0)); if(CLi < 0.2) CLi = 0.2; double CP = CENT / V; double F_ORAL = F1; $CAPTURE CP CLi F_ORAL ' mod_aggregate <- mcode("voriconazole_aggregate_312", aggregate_model_code, quiet = TRUE) mod_m1 <- mcode("voriconazole_m1_312", m1_model_code, quiet = TRUE) mod_m2 <- mcode("voriconazole_m2_312", m2_model_code, quiet = TRUE) model_choices <- c( "Aggregate model (McDougall hybrid, best paper performance)" = "aggregate", "Model 1, Pascual et al. 2012" = "m1", "Model 2, Wang et al. 2014" = "m2" ) model_label <- function(model_choice) { switch( model_choice, aggregate = "Aggregate model", m1 = "Model 1 (Pascual 2012)", m2 = "Model 2 (Wang 2014)", "Aggregate model" ) } make_event_df <- function(route, use_loading, loading_dose, maintenance_dose, dose_unit, wt, interval, n_days, infusion_h) { total_hours <- n_days * 24 dose_to_mg <- function(x) { if (dose_unit == "mg/kg") x * wt else x } build_rows <- function(times, amt_mg, cmt, route_type) { if (length(times) == 0) return(NULL) data.frame( ID = 1, time = times, amt = amt_mg, cmt = cmt, evid = 1, rate = if (route_type == "iv") amt_mg / infusion_h else 0, stringsAsFactors = FALSE ) } oral_cmt <- 1 iv_cmt <- 2 active_cmt <- if (route == "oral") oral_cmt else iv_cmt load_times <- numeric(0) if (isTRUE(use_loading)) { load_times <- seq(0, total_hours - 1e-8, by = interval) load_times <- head(load_times, 2) } maint_start <- if (isTRUE(use_loading) && length(load_times) > 0) max(load_times) + interval else 0 maint_times <- seq(maint_start, total_hours - 1e-8, by = interval) load_df <- build_rows(load_times, dose_to_mg(loading_dose), active_cmt, route) maint_df <- build_rows(maint_times, dose_to_mg(maintenance_dose), active_cmt, route) bind_rows(load_df, maint_df) |> arrange(time) } get_model_object <- function(model_choice) { switch( model_choice, aggregate = mod_aggregate, m1 = mod_m1, m2 = mod_m2, mod_aggregate ) } get_model_params <- function(model_choice, wt, lbw, sex, age, alp, pm, einh, eind, rif, chol) { switch( model_choice, aggregate = list( WT = wt, LBW = lbw, NSEX = if (sex == "Female") 1 else 0, PM = pm, EINH = einh, EIND = eind ), m1 = list( RIF = rif, CHOL = chol ), m2 = list( AGE = age, PM = pm, ALP = alp ), list() ) } simulate_voriconazole <- function(model_choice = "aggregate", route = "oral", use_loading = TRUE, loading_dose = 400, maintenance_dose = 200, dose_unit = "mg", wt = 67, lbw = 54, sex = "Male", age = 60, alp = 104, pm = 0, einh = 0, eind = 0, rif = 0, chol = 0, interval = 12, n_days = 7, infusion_h = 2, delta = 0.25) { event_df <- make_event_df( route = route, use_loading = use_loading, loading_dose = loading_dose, maintenance_dose = maintenance_dose, dose_unit = dose_unit, wt = wt, interval = interval, n_days = n_days, infusion_h = infusion_h ) model <- get_model_object(model_choice) params <- get_model_params(model_choice, wt, lbw, sex, age, alp, pm, einh, eind, rif, chol) out <- model %>% param(params) %>% data_set(event_df) %>% mrgsim(end = n_days * 24, delta = delta) %>% as.data.frame() out$time_h <- out$time out$model_choice <- model_choice out$model_name <- model_label(model_choice) out$route <- route out } compute_metrics <- function(df, interval_h) { last_end <- max(df$time_h, na.rm = TRUE) last_start <- max(0, last_end - interval_h) last_df <- df |> filter(time_h >= last_start, time_h <= last_end) if (nrow(last_df) < 2) { return(list(cmax = NA_real_, ctrough = NA_real_, auc = NA_real_, thalf = NA_real_)) } auc <- sum(diff(last_df$time_h) * (head(last_df$CP, -1) + tail(last_df$CP, -1)) / 2) positive_df <- last_df |> filter(CP > 0) if (nrow(positive_df) > 10) { positive_df <- tail(positive_df, 10) } thalf <- NA_real_ if (nrow(positive_df) >= 3) { fit <- try(lm(log(CP) ~ time_h, data = positive_df), silent = TRUE) if (!inherits(fit, "try-error")) { slope <- coef(fit)[[2]] if (is.finite(slope) && slope < 0) { thalf <- log(2) / abs(slope) } } } list( cmax = max(last_df$CP, na.rm = TRUE), ctrough = min(last_df$CP[last_df$CP > 0], na.rm = TRUE), auc = auc, thalf = thalf ) } app_theme <- bs_theme( version = 5, bootswatch = "flatly", primary = "#8b5cf6" ) |> bs_add_rules(" .metric-card { background: #f8f9fa; border-radius: 8px; padding: 15px; margin: 5px; text-align: center; border: 1px solid #dee2e6; } .metric-value { font-size: 24px; font-weight: bold; color: #2c3e50; } .metric-label { font-size: 12px; color: #7f8c8d; } .metric-success .metric-value { color: #10b981; } .metric-warning .metric-value { color: #f59e0b; } .metric-primary .metric-value { color: #8b5cf6; } .metric-info .metric-value { color: #0dcaf0; } .small-note { font-size: 12px; color: #6c757d; line-height: 1.5; } .ref-box { background: #f0f4ff; border-left: 4px solid #8b5cf6; padding: 12px 16px; border-radius: 4px; margin-top: 10px; font-size: 13px; } .ref-box a { color: #8b5cf6; } ") ui <- page_sidebar( title = "Voriconazole MIPD Simulator", theme = app_theme, sidebar = sidebar( title = "Simulation Settings", width = 360, accordion( open = c("Model", "Dosing"), accordion_panel( "Model", selectInput("model_choice", "Population model", choices = model_choices, selected = "aggregate"), radioButtons("route", "Route", choices = c("Oral" = "oral", "IV infusion" = "iv"), selected = "oral"), checkboxInput("log_scale", "Log scale", value = FALSE), div(class = "small-note", "The aggregate model generally performed best in the paper. M1 and M2 are included for comparison with the published standard NLME models." ) ), accordion_panel( "Dosing", checkboxInput("use_loading", "Use 2 loading doses", value = TRUE), selectInput("dose_unit", "Dose units", choices = c("mg" = "mg", "mg/kg" = "mg/kg"), selected = "mg"), numericInput("loading_dose", "Loading dose", value = 400, min = 1, step = 10), numericInput("maintenance_dose", "Maintenance dose", value = 200, min = 1, step = 10), sliderInput("interval", "Dosing interval (h)", min = 8, max = 24, value = 12, step = 1), numericInput("n_days", "Simulation duration (days)", value = 7, min = 2, max = 30, step = 1), numericInput("infusion_h", "IV infusion duration (h)", value = 2, min = 0.25, max = 6, step = 0.25), div(class = "small-note", "Paper regimens: oral 2 × 400 mg q12h, then 200 mg q12h. IV 2 × 6 mg/kg q12h, then 4 mg/kg q12h over 2 h." ) ), accordion_panel( "Patient characteristics", sliderInput("wt", "Total body weight (kg)", min = 35, max = 120, value = 67, step = 1), sliderInput("lbw", "Lean body weight (kg)", min = 25, max = 90, value = 54, step = 1), radioButtons("sex", "Sex", choices = c("Male", "Female"), selected = "Male", inline = TRUE), sliderInput("age", "Age (years)", min = 18, max = 90, value = 60, step = 1), numericInput("alp", "Alkaline phosphatase (U/L)", value = 104, min = 10, step = 5), checkboxInput("pm", "CYP2C19 poor metabolizer", value = FALSE), checkboxInput("einh", "Enzyme inhibitor coadministration", value = FALSE), checkboxInput("eind", "Enzyme inducer coadministration", value = FALSE), checkboxInput("rif", "Rifampicin coadministration", value = FALSE), checkboxInput("chol", "Grade 3 cholestasis", value = FALSE), div(class = "small-note", "Aggregate model uses weight, lean body weight, sex, CYP2C19 PM, and enzyme inhibitor/inducer status. M1 uses rifampicin and cholestasis. M2 uses age, ALP, and CYP2C19 PM." ) ) ) ), layout_column_wrap( width = 1/4, fill = FALSE, div(class = "metric-card metric-success", div(class = "metric-value", textOutput("cmax")), div(class = "metric-label", HTML("Cmax, last interval (mg/L)")) ), div(class = "metric-card metric-warning", div(class = "metric-value", textOutput("ctrough")), div(class = "metric-label", HTML("Ctrough, last interval (mg/L)")) ), div(class = "metric-card metric-primary", div(class = "metric-value", textOutput("auc")), div(class = "metric-label", HTML("AUC0-τ (mg·h/L)")) ), div(class = "metric-card metric-info", div(class = "metric-value", textOutput("thalf")), div(class = "metric-label", HTML("Apparent t½ (h)")) ) ), card( full_screen = TRUE, plotOutput("pkPlot", height = "500px") ), navset_card_underline( nav_panel( "Model Information", markdown("## Voriconazole MIPD simulator **Source paper:** Kluwe F, Michelet R, Huisinga W, Zeitlinger M, Mikus G, Kloft C. *Towards model-informed precision dosing of voriconazole: Challenging published voriconazole nonlinear mixed-effects models with real-world clinical data.* Clinical Pharmacokinetics. ### Structural models implemented **Aggregate model (best paper performance)** - Two-compartment model with first-order absorption - Parallel linear and nonlinear elimination - Time-dependent reduction in Vmax to reflect autoinhibition - Covariates: weight, lean body weight, sex, CYP2C19 poor metabolizer status, enzyme inhibitor or inducer coadministration **Model 1, Pascual et al. 2012** - One-compartment model with first-order absorption and linear elimination - CL = 5.2 L/h, V = 92 L, Ka = 1.1 h⁻¹, F = 0.63 - Covariates: rifampicin coadministration and grade 3 cholestasis **Model 2, Wang et al. 2014** - One-compartment model with first-order absorption and linear elimination - CL = 6.95 L/h, V = 200 L, Ka = 1.1 h⁻¹, F = 0.895 - Covariates: age, alkaline phosphatase, CYP2C19 poor metabolizer status ### Clinical targets used in the paper - Frequently used voriconazole trough target range: **1 to 6 mg/L** - Paper example: with the aggregate model, standard oral maintenance 200 mg q12h after loading may underperform, while higher maintenance doses can improve target attainment ### Default regimen in this simulator - Oral: **2 × 400 mg q12h**, then **200 mg q12h** - IV infusion: **2 × 6 mg/kg q12h**, then **4 mg/kg q12h over 2 h** ### Notes - This app is for research and educational use only - It does not implement full Bayesian MAP updating from measured TDM samples - The aggregate model is the most appropriate starting point for exploring MIPD concepts from this paper ") ), nav_panel( "References", div(class = "ref-box", tags$h5("Key references"), tags$ol( tags$li("Kluwe F et al. Towards model-informed precision dosing of voriconazole: Challenging published voriconazole nonlinear mixed-effects models with real-world clinical data. Clin Pharmacokinet."), tags$li("McDougall DAJ, Martin J, Playford EG, Green B. Determination of a suitable voriconazole pharmacokinetic model for personalised dosing. J Pharmacokinet Pharmacodyn. 2016."), tags$li("Pascual A et al. Challenging recommended oral and intravenous voriconazole doses for improved efficacy and safety: population pharmacokinetics-based analysis of adult patients with invasive fungal infections. Clin Infect Dis. 2012."), tags$li("Wang T et al. Identification of factors influencing the pharmacokinetics of voriconazole and optimization of dosage regimens based on Monte Carlo simulation in patients with invasive fungal infections. J Antimicrob Chemother. 2014."), tags$li("Voriconazole therapeutic drug monitoring literature supporting target trough concentrations around 1 to 6 mg/L.") ), tags$h5("Drug context"), tags$ul( tags$li(tags$strong("Drug class:"), " Triazole antifungal"), tags$li(tags$strong("Major metabolism:"), " CYP2C19, CYP3A4, CYP2C9"), tags$li(tags$strong("Routes:"), " Oral and IV infusion"), tags$li(tags$strong("Clinical use:"), " Treatment and prophylaxis of invasive fungal infections"), tags$li(tags$strong("Why MIPD matters:"), " High interindividual variability, nonlinear PK, and narrow therapeutic range") ) ) ) ), div(style = "text-align: center; padding: 20px; margin-top: 30px; border-top: 1px solid #e9ecef; color: #6c757d; font-size: 12px;", "Powered by ", tags$a(href = "https://www.pkpdbuilder.com", target = "_blank", style = "color: #8b5cf6; font-weight: 500;", "PKPDBuilder.com"), " • Built by Sunny ☀️ (Husain Attarwala's AI Assistant)", br(), tags$span(style = "font-size: 10px;", "For research and educational purposes only. Not for clinical decision-making.") ) ) server <- function(input, output, session) { observeEvent(input$route, { if (identical(input$route, "oral")) { updateSelectInput(session, "dose_unit", selected = "mg") updateNumericInput(session, "loading_dose", value = 400) updateNumericInput(session, "maintenance_dose", value = 200) updateNumericInput(session, "infusion_h", value = 2) } else { updateSelectInput(session, "dose_unit", selected = "mg/kg") updateNumericInput(session, "loading_dose", value = 6) updateNumericInput(session, "maintenance_dose", value = 4) updateNumericInput(session, "infusion_h", value = 2) } }, ignoreInit = TRUE) sim_data <- reactive({ shiny::req(input$loading_dose, input$maintenance_dose, input$interval, input$n_days) simulate_voriconazole( model_choice = input$model_choice, route = input$route, use_loading = input$use_loading, loading_dose = input$loading_dose, maintenance_dose = input$maintenance_dose, dose_unit = input$dose_unit, wt = input$wt, lbw = input$lbw, sex = input$sex, age = input$age, alp = input$alp, pm = if (isTRUE(input$pm)) 1 else 0, einh = if (isTRUE(input$einh)) 1 else 0, eind = if (isTRUE(input$eind)) 1 else 0, rif = if (isTRUE(input$rif)) 1 else 0, chol = if (isTRUE(input$chol)) 1 else 0, interval = input$interval, n_days = input$n_days, infusion_h = input$infusion_h ) }) metrics <- reactive({ compute_metrics(sim_data(), input$interval) }) output$cmax <- renderText({ sprintf("%.2f", metrics()$cmax) }) output$ctrough <- renderText({ sprintf("%.2f", metrics()$ctrough) }) output$auc <- renderText({ sprintf("%.1f", metrics()$auc) }) output$thalf <- renderText({ if (is.na(metrics()$thalf)) return("--") sprintf("%.1f", metrics()$thalf) }) output$pkPlot <- renderPlot({ df <- sim_data() plot_df <- df if (isTRUE(input$log_scale)) { plot_df <- plot_df |> filter(CP > 0.01) } dose_label <- paste0( if (isTRUE(input$use_loading)) paste0("LD ", input$loading_dose, " ", input$dose_unit, " ×2, then ") else "", input$maintenance_dose, " ", input$dose_unit, " q", input$interval, "h" ) p <- ggplot(plot_df, aes(x = time_h, y = CP)) + annotate("rect", xmin = -Inf, xmax = Inf, ymin = 1, ymax = 6, fill = "#10b981", alpha = 0.10) + geom_hline(yintercept = c(1, 6), linetype = "dashed", color = "#10b981", alpha = 0.75) + annotate("text", x = max(plot_df$time_h) * 0.02, y = 6.35, label = "Target trough range: 1–6 mg/L", hjust = 0, size = 3.5, color = "#10b981") + geom_line(color = "#8b5cf6", linewidth = 0.9) + labs( x = "Time (hours)", y = "Voriconazole concentration (mg/L)", title = paste0("Voriconazole simulator — ", model_label(input$model_choice)), subtitle = paste0(tools::toTitleCase(input$route), " | ", dose_label, " | WT ", input$wt, " kg") ) + theme_minimal(base_size = 14) + theme( plot.title = element_text(face = "bold", size = 15), plot.subtitle = element_text(color = "#6c757d", size = 11) ) if (isTRUE(input$log_scale)) { p <- p + scale_y_log10() } p }) } shinyApp(ui = ui, server = server)