# WHIG Model: Weight-HbA1c-Insulin-Glucose Disease Progression Model
# Reference: Choy et al. CPT Pharmacometrics Syst Pharmacol 2016;5:11-19

library(shiny)
library(bslib)
library(mrgsolve)
library(dplyr)
library(plotly)
library(DT)

model_code <- '
$PARAM @annotated
BLWT     : 104   : Baseline weight (kg)
BLFSI    : 19.2  : Baseline FSI (mIU/mL)
BLFPG    : 7.8   : Baseline FPG (mmol/L)
BLHbA1c  : 6.7   : Baseline HbA1c (pct)
halfWGT  : 96.9  : Half-life of weight compartment (days)
EFDE     : 4.08  : Diet and exercise effect (pct)
EFP      : 2.28  : Placebo effect (pct)
EFUP     : 2.99  : Counter-effect on weight (pct/year)
IS0      : 0.25  : Baseline insulin sensitivity
ScaleEFS : 0.0514 : Scaling factor for weight effect on IS
B0       : 0.61  : Baseline beta-cell function
RB       : 0.209 : Rate of beta-cell decline per year
EFBmax   : 0.171 : Maximal treatment effect on B
SEFBI    : 23.69 : Shape for logistic increase
SEFBD    : 8.05  : Shape for logistic decrease
EFB50    : 190   : Time at half decline of EFB (days)
MTT      : 38.9  : Mean transit time (days)
KinHbA1c : 0.0129 : HbA1c production rate
PPG      : 0.0709 : Postprandial glucose contribution
ScalePPG : 0.963  : PPG scaling factor during treatment
DIET_ON  : 1     : Diet/exercise intervention

$CMT @annotated
WGT      : Body weight (kg)
HBA1     : HbA1c transit 1
HBA2     : HbA1c transit 2
HBA3     : HbA1c transit 3

$MAIN
WGT_0 = BLWT;
double KoutHbA1c = 3.0 / MTT;
double HbA1c_ss_rate = PPG + KinHbA1c * BLFPG;
double HbA1c_ss = HbA1c_ss_rate / KoutHbA1c;
HBA1_0 = HbA1c_ss; HBA2_0 = HbA1c_ss; HBA3_0 = HbA1c_ss;

$ODE
double t_years = SOLVERTIME / 365.0;
double tTRT = 42.0;
double EF_diet_ex = DIET_ON * EFDE / 100.0;
double EF_placebo = (SOLVERTIME >= tTRT) ? (EFP / 100.0) : 0.0;
double EF_up = (EFUP / 100.0) * t_years;
double EF_total = EF_diet_ex + EF_placebo;
double EFW = (1.0 + EF_up) * (1.0 - EF_total);
double KoutWGT = log(2.0) / halfWGT;
double KinWGT = KoutWGT * BLWT;
dxdt_WGT = EFW * KinWGT - KoutWGT * WGT;

double DWGT = BLWT - WGT;
double EFS = 1.0 + ScaleEFS * DWGT;
double IS = IS0 * EFS;

double B = 1.0 / (1.0 + exp(-(log(B0/(1.0-B0)) - RB * t_years)));
double EFBI = 1.0; double EFBD = 1.0;
if (SOLVERTIME > 0 && DIET_ON > 0) {
  EFBI = EFBmax / (1.0 + pow(tTRT / SOLVERTIME, SEFBI));
  EFBD = (1.0 + EFBI) / (1.0 + pow(SOLVERTIME / EFB50, SEFBD));
}
double EFB = 1.0 + EFBI * EFBD - EFBI;

double B_eff = EFB * B;
double alpha = B_eff * 7.8;
double beta = 35.1;
double a = IS * alpha;
double b_coef = -3.5 * IS * alpha;
double c_coef = -beta;
double discriminant = fmax(b_coef * b_coef - 4.0 * a * c_coef, 0.0);
double FPG = (-b_coef + sqrt(discriminant)) / (2.0 * a);
FPG = fmax(FPG, 3.5);
double FSI = alpha * (FPG - 3.5);
FSI = fmax(FSI, 0.1);

double KoutHbA1c_rate = 3.0 / MTT;
double PPG_eff = (SOLVERTIME > 0 && DIET_ON > 0) ? (PPG * ScalePPG) : PPG;
double KinHbA1c_rate = PPG_eff + KinHbA1c * FPG;
dxdt_HBA1 = KinHbA1c_rate - KoutHbA1c_rate * HBA1;
dxdt_HBA2 = KoutHbA1c_rate * HBA1 - KoutHbA1c_rate * HBA2;
dxdt_HBA3 = KoutHbA1c_rate * HBA2 - KoutHbA1c_rate * HBA3;

$TABLE
double HbA1c = HBA1 + HBA2 + HBA3;
double IS_pct = IS * 100.0;
double B_pct = (EFB * B) * 100.0;
double WGTout = WGT;
double WGT_change = BLWT - WGT;
double WGT_change_pct = (WGT_change / BLWT) * 100.0;

$CAPTURE @annotated
WGTout       : Body weight (kg)
WGT_change   : Weight change from baseline (kg)
WGT_change_pct : Weight change (pct)
IS_pct       : Insulin sensitivity (pct of normal)
B_pct        : Beta-cell function (pct of normal)
FSI          : Fasting serum insulin (mIU/mL)
FPG          : Fasting plasma glucose (mmol/L)
HbA1c        : Glycated hemoglobin (pct)
'

mod <- mcode("whig_model", model_code)

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; }
  .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 = "T2DM Disease Progression (WHIG Model)",
  theme = app_theme,
  sidebar = sidebar(
    title = "Simulation Settings", width = 340,
    h6("Patient Characteristics"),
    sliderInput("blwt", "Baseline Weight (kg)", min=70, max=160, value=104, step=1),
    hr(), h6("Baseline Disease State"),
    sliderInput("is0", "Insulin Sensitivity (% normal)", min=10, max=50, value=25, step=1),
    sliderInput("b0", "β-cell Function (% normal)", min=30, max=90, value=61, step=1),
    sliderInput("blfpg", "Fasting Glucose (mmol/L)", min=5.0, max=14.0, value=7.8, step=0.1),
    sliderInput("blhba1c", "HbA1c (%)", min=5.5, max=9.5, value=6.7, step=0.1),
    hr(), h6("Intervention"),
    checkboxInput("diet_on", "Diet & Exercise Program", value=TRUE),
    sliderInput("diet_effect", "Diet/Exercise Intensity (%)", min=0, max=10, value=4, step=0.5),
    hr(),
    sliderInput("sim_weeks", "Duration (weeks)", min=12, max=260, value=67, step=4),
    actionButton("run_sim", "Run Simulation", class="btn-primary btn-block")
  ),
  navset_card_tab(
    title = "WHIG Model Simulator",
    full_screen = TRUE,
    nav_panel("Simulation",
      layout_columns(col_widths=c(3,3,3,3),
        div(class="metric-card", uiOutput("metric_weight")),
        div(class="metric-card", uiOutput("metric_insulin")),
        div(class="metric-card", uiOutput("metric_glucose")),
        div(class="metric-card", uiOutput("metric_hba1c"))
      ),
      tabsetPanel(
        tabPanel("Weight & IS",
          layout_columns(col_widths=c(6,6),
            plotlyOutput("plot_weight", height="420px"),
            plotlyOutput("plot_is", height="420px"))),
        tabPanel("FSI & FPG",
          layout_columns(col_widths=c(6,6),
            plotlyOutput("plot_fsi", height="420px"),
            plotlyOutput("plot_fpg", height="420px"))),
        tabPanel("HbA1c", plotlyOutput("plot_hba1c", height="500px")),
        tabPanel("β-cell Function", plotlyOutput("plot_bcell", height="500px")),
        tabPanel("All Biomarkers", plotlyOutput("plot_all", height="500px")),
        tabPanel("Data", DT::dataTableOutput("data_table"))
      )
    ),
    nav_panel("Compare Scenarios",
      layout_columns(col_widths=c(3,9),
        card(card_header("Scenario Settings"),
          sliderInput("sc_blwt", "Weight (kg)", min=70, max=160, value=104),
          sliderInput("sc_is0", "Insulin Sensitivity (%)", min=10, max=50, value=25),
          sliderInput("sc_b0", "β-cell Function (%)", min=30, max=90, value=61),
          checkboxGroupInput("scenarios", "Compare:",
            choices=c("No intervention"="none","Mild diet (2%)"="mild","Moderate diet (4%)"="mod",
                      "Intensive diet (6%)"="intense","Very intensive (8%)"="very"),
            selected=c("none","mod","intense")),
          sliderInput("sc_weeks", "Duration (weeks)", min=26, max=260, value=104),
          actionButton("run_scenarios", "Compare", class="btn-primary btn-block")
        ),
        div(
          tabsetPanel(
            tabPanel("Weight", plotlyOutput("sc_weight", height="500px")),
            tabPanel("Insulin Sensitivity", plotlyOutput("sc_is", height="500px")),
            tabPanel("HbA1c", plotlyOutput("sc_hba1c", height="500px")),
            tabPanel("Fasting Glucose", plotlyOutput("sc_fpg", height="500px"))
          )
        )
      )
    ),
    nav_panel("Model Information", markdown("
## WHIG Model Overview

The Weight-HbA1c-Insulin-Glucose (WHIG) model describes disease progression in Type 2 Diabetes Mellitus.

### Key Components
- **Weight:** Turnover model with diet/exercise effects
- **Insulin Sensitivity:** Inversely proportional to weight change (~1.5% per kg lost)
- **β-cell Function:** Progressive decline with temporary treatment effect
- **FSI-FPG:** Homeostatic feedback (linked turnover)
- **HbA1c:** 3 transit compartments driven by FPG and PPG
    ")),
    nav_panel("References", div(class="ref-box",
      tags$h5("📚 Key References"),
      tags$ol(tags$li("Choy S, Kjellsson MC, Karlsson MO, de Winter W. CPT Pharmacometrics Syst Pharmacol. 2016;5(1):11-19. ",
                      tags$a(href="https://doi.org/10.1002/psp4.12051", target="_blank", "DOI"))),
      tags$h5("💊 Therapeutic Context"),
      tags$ul(
        tags$li(tags$strong("Disease:"), " Type 2 Diabetes Mellitus"),
        tags$li(tags$strong("Model:"), " Semi-mechanistic disease progression"),
        tags$li(tags$strong("Biomarkers:"), " Weight, HbA1c, FSI, FPG"),
        tags$li(tags$strong("Intervention:"), " Diet & exercise program")
      )
    ))
  ),
  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) {
  sim_counter <- reactiveVal(0)
  observeEvent(input$run_sim, { sim_counter(sim_counter() + 1) })

  sim_data <- reactive({
    sim_counter()
    mod %>% param(BLWT=input$blwt, BLFPG=input$blfpg, BLHbA1c=input$blhba1c,
      IS0=input$is0/100, B0=input$b0/100, EFDE=input$diet_effect,
      DIET_ON=ifelse(input$diet_on,1,0)) %>%
      mrgsim(end=input$sim_weeks*7, delta=1) %>% as.data.frame() %>%
      mutate(week=time/7, month=time/30.44)
  })

  output$metric_weight <- renderUI({
    d <- sim_data(); final <- tail(d,1); bl <- d$WGTout[1]; ch <- final$WGTout-bl
    div(div(class="metric-value", sprintf("%.1f kg",final$WGTout)),
        div(class="metric-label", "Final Weight"),
        div(style="font-size:11px;color:#7f8c8d;", sprintf("%+.1f kg",ch)))
  })
  output$metric_insulin <- renderUI({
    d <- sim_data(); final <- tail(d,1)
    div(div(class="metric-value", sprintf("%.1f%%",final$IS_pct)), div(class="metric-label","Insulin Sensitivity"))
  })
  output$metric_glucose <- renderUI({
    d <- sim_data(); final <- tail(d,1)
    div(div(class="metric-value", sprintf("%.1f",final$FPG)), div(class="metric-label","Fasting Glucose (mmol/L)"))
  })
  output$metric_hba1c <- renderUI({
    d <- sim_data(); final <- tail(d,1)
    div(div(class="metric-value", sprintf("%.1f%%",final$HbA1c)), div(class="metric-label","HbA1c"))
  })

  output$plot_weight <- renderPlotly({
    d <- sim_data()
    plot_ly(d, x=~week, y=~WGTout, type="scatter", mode="lines", line=list(color="#f5576c",width=2)) %>%
      layout(title="Body Weight", xaxis=list(title="Time (weeks)"), yaxis=list(title="Weight (kg)"), hovermode="x unified")
  })
  output$plot_is <- renderPlotly({
    d <- sim_data()
    plot_ly(d, x=~week, y=~IS_pct, type="scatter", mode="lines", line=list(color="#4facfe",width=2)) %>%
      layout(title="Insulin Sensitivity", xaxis=list(title="Time (weeks)"), yaxis=list(title="IS (% of normal)"), hovermode="x unified")
  })
  output$plot_fsi <- renderPlotly({
    d <- sim_data()
    plot_ly(d, x=~week, y=~FSI, type="scatter", mode="lines", line=list(color="#764ba2",width=2)) %>%
      layout(title="Fasting Serum Insulin", xaxis=list(title="Time (weeks)"), yaxis=list(title="FSI (mIU/mL)"), hovermode="x unified")
  })
  output$plot_fpg <- renderPlotly({
    d <- sim_data()
    plot_ly(d, x=~week, y=~FPG, type="scatter", mode="lines", line=list(color="#43e97b",width=2)) %>%
      layout(title="Fasting Plasma Glucose", xaxis=list(title="Time (weeks)"), yaxis=list(title="FPG (mmol/L)"), hovermode="x unified")
  })
  output$plot_hba1c <- renderPlotly({
    d <- sim_data()
    plot_ly(d, x=~week, y=~HbA1c, type="scatter", mode="lines", line=list(color="#fa709a",width=3)) %>%
      add_trace(y=rep(6.5,nrow(d)), mode="lines", line=list(dash="dash",color="orange"), name="Diagnostic (6.5%)") %>%
      add_trace(y=rep(7.0,nrow(d)), mode="lines", line=list(dash="dash",color="red"), name="Target (<7%)") %>%
      layout(title="HbA1c", xaxis=list(title="Time (weeks)"), yaxis=list(title="HbA1c (%)"), hovermode="x unified")
  })
  output$plot_bcell <- renderPlotly({
    d <- sim_data()
    plot_ly(d, x=~week, y=~B_pct, type="scatter", mode="lines", line=list(color="#667eea",width=2)) %>%
      layout(title="β-cell Function", xaxis=list(title="Time (weeks)"), yaxis=list(title="β-cell function (% of normal)"), hovermode="x unified")
  })
  output$plot_all <- renderPlotly({
    d <- sim_data() %>% mutate(WGT_n=WGTout/WGTout[1]*100, IS_n=IS_pct/IS_pct[1]*100, FPG_n=FPG/FPG[1]*100, HbA1c_n=HbA1c/HbA1c[1]*100)
    plot_ly(d, x=~week) %>%
      add_trace(y=~WGT_n, name="Weight", mode="lines", line=list(color="#f5576c",width=2)) %>%
      add_trace(y=~IS_n, name="IS", mode="lines", line=list(color="#4facfe",width=2)) %>%
      add_trace(y=~FPG_n, name="FPG", mode="lines", line=list(color="#43e97b",width=2)) %>%
      add_trace(y=~HbA1c_n, name="HbA1c", mode="lines", line=list(color="#fa709a",width=2)) %>%
      layout(title="All Biomarkers (% of Baseline)", xaxis=list(title="Time (weeks)"), yaxis=list(title="% of Baseline"), hovermode="x unified")
  })
  output$data_table <- DT::renderDataTable({
    sim_data() %>% filter(time%%7==0) %>% select(week, WGTout, WGT_change, IS_pct, B_pct, FSI, FPG, HbA1c) %>%
      mutate(across(where(is.numeric), ~round(.,2)))
  }, options=list(pageLength=15, scrollX=TRUE))

  # Scenarios
  sc_counter <- reactiveVal(0)
  observeEvent(input$run_scenarios, { sc_counter(sc_counter() + 1) })

  scenario_data <- reactive({
    sc_counter(); shiny::req(length(input$scenarios)>0)
    scenario_map <- list(none=list(name="No intervention",EFDE=0,DIET_ON=0),mild=list(name="Mild diet (2%)",EFDE=2,DIET_ON=1),
      mod=list(name="Moderate diet (4%)",EFDE=4,DIET_ON=1),intense=list(name="Intensive diet (6%)",EFDE=6,DIET_ON=1),
      very=list(name="Very intensive (8%)",EFDE=8,DIET_ON=1))
    all_data <- lapply(input$scenarios, function(sc) {
      sp <- scenario_map[[sc]]
      mod %>% param(BLWT=input$sc_blwt, IS0=input$sc_is0/100, B0=input$sc_b0/100, EFDE=sp$EFDE, DIET_ON=sp$DIET_ON) %>%
        mrgsim(end=input$sc_weeks*7, delta=7) %>% as.data.frame() %>% mutate(week=time/7, scenario=sp$name)
    })
    bind_rows(all_data)
  })

  output$sc_weight <- renderPlotly({
    d <- scenario_data()
    plot_ly(d, x=~week, y=~WGTout, color=~scenario, type="scatter", mode="lines", line=list(width=2)) %>%
      layout(title="Weight by Intervention", xaxis=list(title="Time (weeks)"), yaxis=list(title="Weight (kg)"), hovermode="x unified")
  })
  output$sc_is <- renderPlotly({
    d <- scenario_data()
    plot_ly(d, x=~week, y=~IS_pct, color=~scenario, type="scatter", mode="lines", line=list(width=2)) %>%
      layout(title="Insulin Sensitivity by Intervention", xaxis=list(title="Time (weeks)"), yaxis=list(title="IS (% normal)"), hovermode="x unified")
  })
  output$sc_hba1c <- renderPlotly({
    d <- scenario_data()
    plot_ly(d, x=~week, y=~HbA1c, color=~scenario, type="scatter", mode="lines", line=list(width=2)) %>%
      layout(title="HbA1c by Intervention", xaxis=list(title="Time (weeks)"), yaxis=list(title="HbA1c (%)"), hovermode="x unified")
  })
  output$sc_fpg <- renderPlotly({
    d <- scenario_data()
    plot_ly(d, x=~week, y=~FPG, color=~scenario, type="scatter", mode="lines", line=list(width=2)) %>%
      layout(title="Fasting Glucose by Intervention", xaxis=list(title="Time (weeks)"), yaxis=list(title="FPG (mmol/L)"), hovermode="x unified")
  })
}

shinyApp(ui = ui, server = server)