# Script generated by the UniPV WinBUGS converter
# 'WinBUGSConverter' copyright (c) BMS - University of Pavia (2016)
# Converter Version: 1.0.0
# Source: PharmML 0.8.1
# Target: WINBUGS 1.4.3
# Run ID: run_id_0
# Model: Generated from MDL. MOG ID: warfarin_PK_ODE_mog
# File: UseCase3
# Dated: Mon Aug 01 23:09:59 CEST 2016

model {
  for (ind_subj in 1:N_subj) {
    for (ind_t in 1:N_t[ind_subj]) {
      WTinterp[ind_subj,ind_t] <- interp.function(grid[ind_subj,ind_t], grid_WT[ind_subj,1:N_t_WT[ind_subj]], WT[ind_subj,1:N_t_WT[ind_subj]])
      logtWTinterp[ind_subj,ind_t] <- (log((WTinterp[ind_subj,ind_t] / 70)))
      CP_obs_mean[ind_subj,ind_t] <- (CC[ind_subj,ind_t] + 0)
      CP_obs_prec[ind_subj,ind_t] <- 1 / ((pow(((RUV_ADD + (RUV_PROP * CC[ind_subj,ind_t]))),2))*1)
      CP_obs[ind_subj,ind_t] ~ dnorm(CP_obs_mean[ind_subj,ind_t], CP_obs_prec[ind_subj,ind_t])
      CP_obs_res[ind_subj,ind_t] <- CP_obs[ind_subj,ind_t]-CP_obs_mean[ind_subj,ind_t]
      PCA_obs_mean[ind_subj,ind_t] <- (PCA[ind_subj,ind_t] + 0)
      PCA_obs[ind_subj,ind_t] ~ dnorm(PCA_obs_mean[ind_subj,ind_t], PCA_obs_prec)
      PCA_obs_res[ind_subj,ind_t] <- PCA_obs[ind_subj,ind_t]-PCA_obs_mean[ind_subj,ind_t]
      log(CL[ind_subj,ind_t]) <- log(POP_CL) + BETA_CL_WT * logtWTinterp[ind_subj,ind_t] + ETA_CL[ind_subj]
      log(V[ind_subj,ind_t]) <- log(POP_V) + BETA_V_WT * logtWTinterp[ind_subj,ind_t] + ETA_V[ind_subj]
    }
    ETA_KA[ind_subj] ~ dnorm(0, ETA_KA_prec)
    ETA_TLAG[ind_subj] ~ dnorm(0, ETA_TLAG_prec)
    ETA_PCA0[ind_subj] ~ dnorm(0, ETA_PCA0_prec)
    ETA_EMAX[ind_subj] ~ dnorm(0, ETA_EMAX_prec)
    ETA_C50[ind_subj] ~ dnorm(0, ETA_C50_prec)
    ETA_TEQ[ind_subj] ~ dnorm(0, ETA_TEQ_prec)
    log(KA[ind_subj])<- log(POP_KA) + ETA_KA[ind_subj]
    log(TLAG[ind_subj])<- log(POP_TLAG) + ETA_TLAG[ind_subj]
    log(PCA0[ind_subj])<- log(POP_PCA0) + ETA_PCA0[ind_subj]
    EMAX[ind_subj] <- POP_EMAX + ETA_EMAX[ind_subj]
    log(C50[ind_subj])<- log(POP_C50) + ETA_C50[ind_subj]
    log(TEQ[ind_subj])<- log(POP_TEQ) + ETA_TEQ[ind_subj]
    ETA1[ind_subj,1:2] ~ dmnorm(mu1[], T1[,])
    ETA_CL[ind_subj] <- ETA1[ind_subj,1]
    ETA_V[ind_subj] <- ETA1[ind_subj,2]
    der99wb_unipv[ind_subj,1:N_t[ind_subj], 1:3]<- function.model1(grid[ind_subj,1:N_t[ind_subj]], amt[ind_subj,1:N_t[ind_subj]],  rate[ind_subj,1:N_t[ind_subj]], ii[ind_subj,1:N_t[ind_subj]], evid[ind_subj,1:N_t[ind_subj]], cmt[ind_subj,1:N_t[ind_subj]], addl[ind_subj,1:N_t[ind_subj]], ss[ind_subj,1:N_t[ind_subj]], theta[ind_subj,])
    theta[ind_subj,1] <- n_cov_cont
    theta[ind_subj,2] <- max_m_WT
    theta[ind_subj,3] <- N_t_WT[ind_subj]
    theta[ind_subj,3+2*max_m_WT +1] <- ETA_KA[ind_subj]
    theta[ind_subj,3+2*max_m_WT +2] <- POP_KA
    theta[ind_subj,3+2*max_m_WT +3] <- ETA_CL[ind_subj]
    theta[ind_subj,3+2*max_m_WT +4] <- POP_CL
    theta[ind_subj,3+2*max_m_WT +5] <- BETA_CL_WT
    theta[ind_subj,3+2*max_m_WT +6] <- ETA_V[ind_subj]
    theta[ind_subj,3+2*max_m_WT +7] <- POP_V
    theta[ind_subj,3+2*max_m_WT +8] <- BETA_V_WT
    theta[ind_subj,3+2*max_m_WT +9] <- ETA_PCA0[ind_subj]
    theta[ind_subj,3+2*max_m_WT +10] <- POP_PCA0
    theta[ind_subj,3+2*max_m_WT +11] <- ETA_EMAX[ind_subj]
    theta[ind_subj,3+2*max_m_WT +12] <- POP_EMAX
    theta[ind_subj,3+2*max_m_WT +13] <- ETA_C50[ind_subj]
    theta[ind_subj,3+2*max_m_WT +14] <- POP_C50
    theta[ind_subj,3+2*max_m_WT +15] <- ETA_TEQ[ind_subj]
    theta[ind_subj,3+2*max_m_WT +16] <- POP_TEQ
    for (i in 1 : max_m_WT){ 
	theta[ind_subj,i+3] <- grid_WT[ind_subj,i]
	}
	for (i in 1 : max_m_WT){ 
	theta[ind_subj,i+3+max_m_WT] <- WT[ind_subj,i]
	}
      theta[ind_subj,3+2*max_m_WT+17] <- PCA0[ind_subj] #IC PCA
      theta[ind_subj,3+2*max_m_WT+18] <- 0.0 #IC Ad1
      theta[ind_subj,3+2*max_m_WT+19] <- 0.0 #IC CENTRAL
      theta[ind_subj,3+2*max_m_WT+20] <- 1 #F PCA
      theta[ind_subj,3+2*max_m_WT+21] <- 1 #F Ad1
      theta[ind_subj,3+2*max_m_WT+22] <- 1 #F CENTRAL
      theta[ind_subj,3+2*max_m_WT+23] <- 0 #TLAG PCA
      theta[ind_subj,3+2*max_m_WT+24] <- TLAG[ind_subj] #TLAG Ad1
      theta[ind_subj,3+2*max_m_WT+25] <- 0 #TLAG CENTRAL
      for (ind_q in 1:N_t[ind_subj]) {
        PCA[ind_subj,ind_q] <- der99wb_unipv[ind_subj,ind_q,1] + PCA0[ind_subj]
        Ad1[ind_subj,ind_q] <- der99wb_unipv[ind_subj,ind_q,2] + 0.0
        CENTRAL[ind_subj,ind_q] <- der99wb_unipv[ind_subj,ind_q,3] + 0.0
        PD[ind_subj,ind_q] <- (1 - ((EMAX[ind_subj] * CC[ind_subj,ind_q]) / (C50[ind_subj] + CC[ind_subj,ind_q])))
        DPCA[ind_subj,ind_q] <- PCA[ind_subj,ind_q]
        CC[ind_subj,ind_q] <- (CENTRAL[ind_subj,ind_q] / V[ind_subj,ind_q])
        RPCA[ind_subj,ind_q] <- (PCA0[ind_subj] * KPCA[ind_subj,ind_q])
        KPCA[ind_subj,ind_q] <- ((log(2)) / TEQ[ind_subj])
      }
    }
    PCA_obs_prec <- 1 / ((pow((RUV_FX),2))*1)
    ETA_KA_prec <- (1/(pow(PPV_KA,2)))
    ETA_TLAG_prec <- (1/(pow(PPV_TLAG,2)))
    ETA_PCA0_prec <- (1/(pow(PPV_PCA0,2)))
    ETA_EMAX_prec <- (1/(pow(PPV_EMAX,2)))
    ETA_C50_prec <- (1/(pow(PPV_C50,2)))
    ETA_TEQ_prec <- (1/(pow(PPV_TEQ,2)))
    T1[1:2,1:2] <- inverse(OMEGA1[,])
    mu1[1] <- 0.0
    mu1[2] <- 0.0
    OMEGA1[1,1] <- pow(PPV_CL,2)
    OMEGA1[1,2] <- CORR_CL_V*PPV_CL*PPV_V
    OMEGA1[2,1] <- CORR_CL_V*PPV_CL*PPV_V
    OMEGA1[2,2] <- pow(PPV_V,2)
    MU_CL <- (log(0.1))
    MU_V <- (log(8))
    MU_KA <- (log(0.362))
    MU_TLAG <- (log(1))
    MU_PPV_CL <- (log(0.1))
    MU_PPV_V <- (log(0.1))
    MU_PPV_KA <- (log(0.1))
    MU_CORR_CL_V <- (log(0.01))
    MU_RUV_PROP <- (log(0.1))
    MU_RUV_ADD <- (log(0.1))
    VAR_CL <- pow((log(pow(0.1,2))),2)
    VAR_V <- pow((log(pow(8,2))),2)
    VAR_KA <- pow((log(pow(0.362,2))),2)
    VAR_TLAG <- pow((log(pow(1,2))),2)
    VAR_PPV_CL <- pow((log(pow(0.1,2))),2)
    VAR_PPV_V <- pow((log(pow(0.1,2))),2)
    VAR_PPV_KA <- pow((log(pow(0.1,2))),2)
    VAR_CORR_CL_V <- pow((log(pow(0.01,2))),2)
    VAR_RUV_PROP <- pow((log(pow(0.1,2))),2)
    VAR_RUV_ADD <- pow((log(pow(0.1,2))),2)
    MU_PCA0 <- (log(96.7))
    MU_C50 <- (log(1.2))
    MU_TEQ <- (log(12.9))
    MU_PPV_PCA0 <- (log(0.1))
    MU_PPV_C50 <- (log(0.1))
    MU_PPV_TEQ <- (log(0.1))
    MU_RUV_FX <- (log(4))
    VAR_PCA0 <- pow((log(pow(96.7,2))),2)
    VAR_C50 <- pow((log(pow(1.2,2))),2)
    VAR_TEQ <- pow((log(pow(12.9,2))),2)
    VAR_PPV_PCA0 <- pow((log(pow(0.1,2))),2)
    VAR_PPV_C50 <- pow((log(pow(0.1,2))),2)
    VAR_PPV_TEQ <- pow((log(pow(0.1,2))),2)
    VAR_RUV_FX <- pow((log(pow(4,2))),2)
    BETA_CL_WT <- 0.75
    BETA_V_WT <- 1
    POP_EMAX <- 1
    PPV_TLAG <- 0.1
    PPV_EMAX <- 0.1
    POP_CL <- 0.1
    POP_V <- 8
    POP_KA <- 0.36
    POP_TLAG <- 1
    PPV_CL <- 0.1
    PPV_V <- 0.14
    PPV_KA <- 0.94
    CORR_CL_V <- 0.29
    RUV_PROP <- 0.11
    RUV_ADD <- 0.001
    POP_PCA0 <- 96
    POP_C50 <- 1.2
    POP_TEQ <- 12
    PPV_PCA0 <- 0.1
    PPV_C50 <- 0.1
    PPV_TEQ <- 0.1
    RUV_FX <- 4
    # PRIORS
}
