
#clear workspace
rm(list = ls())

library(deSolve)

#set the working directory to where you downloaded the scripts
this.dir <- dirname(parent.frame(2)$ofile)
setwd(this.dir)


#Get model File
source("..\\Model\\model-new-annotated.R")

#get paramFile
chem = "Carbaryl"
gender = "M"
age = "25_R_Y" 

param_file_name = paste("..",chem,paste("params",gender,paste(age,".R","sep" = ""),sep="_"),sep = "/")
source(param_file_name)

#set dosing parameters only keep one for the scenario file


#inhalation
params[["INHDOSE"]] <-7.1
params[["INHTLEN"]] <-1
params[["INHDAYS"]] <-7


#set simulation Parameters

#simulation
params[["TSTART"]] <-0
params[["TOTDAYS"]] <-180

initial_params <- within(as.list(params),{
  CINH <- INHDOSE/24450.0*1000
  QALV <- (TV-DS)*RESPR
  RDRMEXP <- RDRMEXPC/MW
  EXPAREA <- EXPAREAC*BSA*10000 #cm^2
  # PARAMETER SCALING
  #TISSUE VOLUME CHECK
  VOLTOTALC <- VOLBLOODC + VOLFATC+ VOLBRAINC + VOLLIVERC +VOLRESTC
  
  # SCALED TISSUE VOLUMES (L)
  VOLBLOOD <- VOLBLOODC *(0.85/VOLTOTALC)* BODYWT    # L; BLOOD
  VOLFAT <- VOLFATC *(0.85/VOLTOTALC)* BODYWT      # L; FATOSE TISSUE
  VOLBRAIN <- VOLBRAINC *(0.85/VOLTOTALC)* BODYWT    # L; BRAIN TISSUE
  VOLLIVER <- VOLLIVERC *(0.85/VOLTOTALC)* BODYWT    # L; LIVER TISSUE
  VOLREST <-  VOLRESTC*(0.85/VOLTOTALC)*BODYWT       # L; Rest of the body
  VOLPLS <- VOLBLOOD*(1-HCT)                         # L; Plasma
  VOLRBC <- VOLBLOOD*HCT                             # L; RBC
  
  #BLOOD FLOW CHECK
  FRTOTALC <- FRFATC + FRBRAINC  + FRLIVC + FRRESTC
  
  # SCALED BLOOD FLOWS (L / HR)
  CARDOUTP <- CARDOUTPC * (1-HCT)               # L/HR; CARDIAC OUTPUT # use different type of QC
  QFAT <- FRFATC *(1/FRTOTALC)* CARDOUTP      # L/HR; FATOSE BLOOD FLOW
  QBRAIN <- FRBRAINC *(1/FRTOTALC)* CARDOUTP    # L/HR; BRAIN BLOOD FLOW
  QLIV <- FRLIVC * (1/FRTOTALC)*CARDOUTP        # L/HR; LIVER ARTERIAL BLOOD FLOW
  QREST <- FRRESTC* (1/FRTOTALC)*CARDOUTP       # L/HR; Rest of the body
  
  #permiability constants
  
  PAF <- PAFC* (VOLFAT ^0.75)	                # fat
  PAR <- PARC* (VOLREST ^0.75)	                # remaining
  PAL <- PALC* (VOLLIVER ^0.75)                	# liver
  PAB <- PABC* (VOLBRAIN ^0.75)                	# brain
  
  #carbaryl metabolism
  
  VKM1  <- VKM1C  * VOLLIVER					          # carbaryl CLint in liver
  VKM2  <- VKM2C  * (VOLBLOOD*(1-HCT))          # carbaryl CLint in plasma, VPLS= (VBLD*(1-HCT))
  UCLN  <- UCLNC  * BODYWT^(-0.25)	          	# carbaryl rate of excretion from plasma to urine
  
  #Synthesis Rate of ACHe add
  
  KsAChB <- (BACh*VOLBRAIN/TRACh)*KdAChB        #Brain (umole/hr)
  KsAChPLS <- (PLSACh*VOLPLS/TRACh)*KdAChPLS    #Plasma  
  KsAChRBC <- (RBCACh*VOLRBC/TRACh)*KdAChRBC    #RBC  
  
  #Synthesis Rate of BChe Add
  KsBChPLS <- (PLSBCh*VOLPLS/TRBCh)*KdBChPLS    #Plasma
  
  #Simulation Paramters
  TSTOP <- TOTDAYS*24
  
})


ABACh0 <- initial_params[["BACh"]]*initial_params[["VOLBRAIN"]]/initial_params[["TRACh"]]
APLSACh0 <- initial_params[["PLSACh"]]*initial_params[["VOLPLS"]]/initial_params[["TRACh"]]
APLSBCh0 <- initial_params[["PLSBCh"]]*initial_params[["VOLPLS"]]/initial_params[["TRBCh"]]
ARBCACh0 <- initial_params[["RBCACh"]]*initial_params[["VOLRBC"]]/initial_params[["TRACh"]]



#function for dosing
bw <- initial_params[["BODYWT"]]
mw <- initial_params[["MW"]]

#ORAL
bdose <- initial_params[["BDOSE"]]
breps <- initial_params[["BREPS"]]
blen <- initial_params[["BLEN"]]

totbreps <- breps*blen
#Drinking Water
ddose <- initial_params[["DDOSE"]]
vdw <- initial_params[["VDW"]]
dreps <- initial_params[["DREPS"]]

#inhalation
inhdose <- initial_params[["INHDOSE"]]
inhtlen <- initial_params[["INHTLEN"]]
inhdays <- initial_params[["INHDAYS"]]

#dermal
drmdlen  <- initial_params[["DRMDLEN"]]

#simulation
tstart <- initial_params[["TSTART"]]
numdays <- initial_params[["TOTDAYS"]]
tstop <- initial_params[["TSTOP"]]


#dermal
drmdlen  <- initial_params[["DRMDLEN"]] # Time for 100% dermal absorption
drmelen <- initial_params[["DRMELEN"]]  # Total Length of Dermal Exposure (Initial to termination) 

# If time to 100% absortion is greater than total length of exposure, switch exsposure off at total length of exposure.
if (drmdlen>drmelen){
  totdrmlen <- drmelen
}else{
  totdrmlen <- drmdlen
}

#if bolus oral dose is administered
if (bdose > 0){
  # var to change
  state_Var = c("ODOSE","TOTODOSE")
  # Value  of change 
  change_val1= (bdose*bw*1000/mw)/totbreps
  change_val2= change_val1
  change_arr = c(change_val1,change_val2)
  # operation of event
  operation = c("add","add")
  # times of event
  if (breps==1){
    #only one bolus dose per day
    event_times <- head(seq(tstart,tstop,24),-1)
  }else{
    #multiple bolus doses per day
    event_times <- unlist(lapply(X = 1:numdays, FUN = function(x){head(seq(0,blen,1/breps),-1)+(24*(x-1))}))
  }
  
  eventDat <- data.frame(
    
    var = rep(x = state_Var,each = length(event_times)),
    time = rep(event_times,length(state_Var)),
    value = rep(x = change_arr,each = length(event_times)),
    method = rep(x = operation,each = length(event_times))
    
  )
  
  # if drinking water dose is administered
}else if (ddose >0){
  # var to change
  state_Var = c("DDOSE","TOTDDOSE")
  # Value  of change 
  change_val1= (ddose*vdw*1000/mw)/dreps
  change_val2= change_val1
  change_arr = c(change_val1,change_val2)
  # operation of event
  operation = c("add","add")
  # times of event
  event_times = unlist(lapply(X = 1:numdays,function(x){head(seq(0,24,by = 24/dreps),-1)+24*(x-1)}))
  
  eventDat <- data.frame(
    
    var = rep(x = state_Var,each = length(event_times)),
    time = rep(event_times,length(state_Var)),
    value = rep(x = change_arr,each = length(event_times)),
    method = rep(x = operation,each = length(event_times))
    
  )
  
  # if inhalation dose is administered  
}else if (inhdose >0){
  # var to change
  state_var1 = "INHSWTCH"
  state_var2 ="INHSWTCH"
  # Value  of change 
  change_val1= 1
  change_val2= 0
  # operation of event
  operation1 = "rep"
  operation2 = "rep"
  # times of event
  
  #days on which dosing can occue
  event_days = unlist(lapply(X=1:numdays,function(x){lapply(1:inhdays,function(y){(x-1)*7+y})}))
  
  event_times1 = unlist(lapply(event_days,function(x){0+24*(x-1)}))
  event_times2 = unlist(lapply(event_days,function(x){inhtlen+24*(x-1)}))
  eventDat <- data.frame(
    
    var = c(rep(x = state_var1,each = length(event_times1)),rep(x = state_var2,each = length(event_times2))),
    time = c(event_times1,event_times2),
    value = c(rep(x = change_val1,each = length(event_times1)),rep(x = change_val2,each = length(event_times2))),
    method = c(rep(x = operation1,each = length(event_times1)),rep(x = operation2,each = length(event_times2)))
    
  )
  
  # if dermal exposure (not used in pyrethroid model)
}else if (totdrmlen >0){
  # var to change
  state_var1 = "DRMSWTCH"
  state_var2 ="DRMSWTCH"
  # Value  of change 
  change_val1= 1
  change_val2= 0
  # operation of event
  operation1 = "rep"
  operation2 = "rep"
  # times of event
  event_times1 = unlist(lapply(1:numdays,function(x){0+24*(x-1)}))
  event_times2 = unlist(lapply(1:numdays,function(x){totdrmlen+24*(x-1)}))
  eventDat <- data.frame(
    
    var = c(rep(x = state_var1,each = length(event_times1)),rep(x = state_var2,each = length(event_times2))),
    time = c(event_times1,event_times2),
    value = c(rep(x = change_val1,each = length(event_times1)),rep(x = change_val2,each = length(event_times2))),
    method = c(rep(x = operation1,each = length(event_times1)),rep(x = operation2,each = length(event_times2)))
    
  )
  
}else{
  print("No exposure selected")
}


state <- c(INHSWTCH=0,DRMSWTCH=0,ODOSE=0,TOTODOSE=0,DDOSE = 0, TOTDDOSE = 0,AINH = 0,AEXH=0,ADRM=0,ANABSC=0,AABSC = 0,AABSN=0,AVL=0, AL=0,LAM=0, AVF=0, AF=0, AVR=0, AR=0,ABI=0,ABACh=ABACh0,ABAChi=0,
           AVB=0, AB=0, PLSAM=0,APLS=0,APLSACh=APLSACh0,APLSAChi=0,APLSBCh=APLSBCh0,APLSBChi=0, ARBC = 0,ARBCACh=ARBCACh0,
           ARBCAChi=0,ABLDI = 0,AAN=0,AUN=0,AUCBRAIN=0, AUCLIVER=0, AUCFATOSE=0, AUCBLOOD=0,AUCRBC=0)


#total simulation output times
times <- seq(tstart,tstop,by=0.01)


modelOutput<- ode(y = state, times = times,method = "lsodes",
                  func = genericCarbarylModel, parms = initial_params,
                  events = list(data = eventDat))

result <- as.data.frame(modelOutput)
print("RBC")
print(100- min(result$pctRBCACH))
print("BRAIN")
print(100-min(result$pctBACh))

calcDermDose <- drmdlen*initial_params[["RDRMEXPC"]]
print("Dermal Dose")
print(calcDermDose)