# SIMEONI 2004 with IIV
# TGI model in vivo
# This model is used for estimation.
#
# Limitations:
# - among the drugs considered in the paper, only PACLITAXEL has been used for PK
#
# Dataset:
# Simulated dataset
#
# Evaluation method for this model:
# run NONMEM estimation via TEL;
# compare individual fits with figures in the paper (FIG.7)
#
# Original by Elisa Borella, Letizia Carrara, Silvia Lavezzi, Paolo Magni, Giulia Muzio, Lorenzo Pasotti, Elena Tosca
# BMS lab - Univ Pavia
# Copyright 2015
# http://repository.ddmore.eu/model/DDMODEL00000006
#

simeoni2004_data = dataObj{

DECLARED_VARIABLES{Y::observation Q1::dosingTarget}

   DATA_INPUT_VARIABLES{
      ID:{use is id}
      TIME:{use is idv}
      DV:{use is dv, variable= Y}
      AMT:{use is amt, variable=Q1}
      CMT:{use is cmt}
      MDV:{use is mdv}
      GRP:{use is catCov withCategories {Control when 0, Treated when 1}}
   }
   
   SOURCE{
   	scrfile : {
       file="simeoni_pop.csv",  
       inputFormat is nonmemFormat
       
      }
   }
}

simeoni2004_par = parObj{

   STRUCTURAL{
      LAMBDA0_POP : {value=0.3, lo=0 } # units="1/day"
      LAMBDA1_POP : {value=0.7, lo=0} #  units="g/day"
      K1_POP : {value=0.7, lo=0 } # units="1/day"
      K2_POP : {value=0.5, lo=0} # units="L/mg/day"
      W0_POP : {value=0.02, lo=0 } # units="g" 
      K10_POP : {value=20.832, fix=true} # 0.868*24 1/day
	  K12_POP : {value=0.144, fix=true} # 0.006*24 1/day
      K21_POP : {value = 2.011, fix=true} # 0.0838*24 1/day
      V1_POP : {value=0.81, fix=true} #L/kg
      PSI_POP : {value=20, fix=true}
   	}
   
   VARIABILITY{
      CV : {value=0.1, lo=0}
      OMEGA_LAMBDA0 : {value=0.02}
      OMEGA_LAMBDA1 : {value=0.3}      
      }
}

simeoni2004_pri = priorObj {
    PRIOR_VARIABLE_DEFINITION {
        logLAMBDA0_POP ~ Normal(mean=MU_LAMBDA0_POP, var=SIGMA_LAMBDA0_POP)
        logLAMBDA1_POP ~ Normal(mean=MU_LAMBDA1_POP, var=SIGMA_LAMBDA1_POP)
        logK1_POP ~ Normal(mean=MU_K1_POP, var=SIGMA_K1_POP)
        logK2_POP ~ Normal(mean=MU_K2_POP, var=SIGMA_K2_POP)
        logW0_POP ~ Normal(mean=MU_W0_POP, var=SIGMA_W0_POP)
        LAMBDA0_POP = exp(logLAMBDA0_POP)
        LAMBDA1_POP = exp(logLAMBDA1_POP)
        K1_POP = exp(logK1_POP)
        K2_POP = exp(logK2_POP)
        W0_POP = exp(logW0_POP)
        
        TAU_LAMBDA0 ~ Gamma2(shape=a_OMEGA_LAMBDA0, rate=b_OMEGA_LAMBDA0)
        OMEGA_LAMBDA0 = 1/TAU_LAMBDA0
        TAU_LAMBDA1 ~ Gamma2(shape=a_OMEGA_LAMBDA1, rate=b_OMEGA_LAMBDA1)
        OMEGA_LAMBDA1 = 1/TAU_LAMBDA1
        
        invCV2 ~ Gamma1(shape=a_CV2, scale=b_CV2)
        CV = sqrt(1/invCV2)
    }

    PRIOR_PARAMETERS {
        MU_LAMBDA0_POP = -1.4
        MU_LAMBDA1_POP = -1.3
        MU_K1_POP = -1.2
        MU_K2_POP = -3.2
        MU_W0_POP = -3.9
        SIGMA_LAMBDA0_POP = 0.8
        SIGMA_LAMBDA1_POP = 1.2
        SIGMA_K1_POP = 1.8
        SIGMA_K2_POP = 1.4
        SIGMA_W0_POP = 3.2

        a_OMEGA_LAMBDA0 = 10
        b_OMEGA_LAMBDA0 = 1      
        a_OMEGA_LAMBDA1 = 3
        b_OMEGA_LAMBDA1 = 1

        a_CV2 = 195
        b_CV2 = 1.4
        K10_POP = 20.832
        K12_POP = 0.144
        K21_POP = 2.011
        V1_POP = 0.81
        PSI_POP = 20
    }
}

simeoni2004_mdl = mdlObj{

 	IDV { T }
 	
 	VARIABILITY_LEVELS{
   	ID : { level=2, type is parameter }
   	DV : { level=1, type is observation }
   }
	
#	COVARIATES{
#		GRP withCategories {Control, Treated}
#	}
	
	STRUCTURAL_PARAMETERS{
      LAMBDA0_POP
      LAMBDA1_POP
      K1_POP
      K2_POP
      W0_POP
      K10_POP
      K12_POP
      K21_POP
      V1_POP
      PSI_POP
	}
	
	VARIABILITY_PARAMETERS{
	  OMEGA_LAMBDA0
      OMEGA_LAMBDA1
	  CV
     }
		
	INDIVIDUAL_VARIABLES{
		LAMBDA0 : {type is linear, trans is ln, pop = LAMBDA0_POP, ranEff = eta_LAMBDA0}
		LAMBDA1 : {type is linear, trans is ln, pop = LAMBDA1_POP, ranEff = eta_LAMBDA1}
		K1 = K1_POP
		K2 = K2_POP
		W0 = W0_POP
		K10 = K10_POP
		K12 = K12_POP
		K21 = K21_POP
		V1 = V1_POP
		PSI = PSI_POP
   	}
   
	
	MODEL_PREDICTION{  	
		DEQ{     	
		    C=Q1/V1
			WTOT=X1+X2+X3+X4
			Q1:{deriv=K21*Q2-(K10+K12)*Q1, init=0}
			Q2:{deriv=K12*Q1-K21*Q2, init=0}
			X1:{deriv=(LAMBDA0*X1/((1+(WTOT*LAMBDA0/LAMBDA1)^PSI)^(1/PSI)))-K2*C*X1, init=W0}
			X2:{deriv=K2*C*X1-K1*X2, init=0}
			X3:{deriv=K1*X2-K1*X3, init=0}
			X4:{deriv=K1*X3-K1*X4, init=0}	
    	}	
	}	
   
   
   	RANDOM_VARIABLE_DEFINITION (level=DV) {
	    eps_RES_W ~ Normal(mean=0, var=1)
	}
	
	RANDOM_VARIABLE_DEFINITION (level=ID) {
	    eta_LAMBDA0 ~ Normal(mean=0, var=OMEGA_LAMBDA0)
	    eta_LAMBDA1 ~ Normal(mean=0, var=OMEGA_LAMBDA1)
	}
	
   OBSERVATION {		 
 	 Y : {type is proportionalError, proportional=CV,  eps=eps_RES_W, prediction= WTOT}
  	}
 
}
   
simeoni2004_NONMEM_task = taskObj{ 
  ESTIMATE { 
        set algo is foce 
	}
}

simeoni2004_Monolix_task = taskObj{ 
  ESTIMATE { 
        set algo is saem 
	}
}
	
simeoni2004_BUGS_task = taskObj{
  ESTIMATE {
        set algo is mcmc
        TARGET_SETTINGS(target="BUGS"){
            set 
            nchains = 1,
            burnin = 1000,
            niter = 5000,
            parameters = "LAMBDA0_POP,LAMBDA1_POP,K1_POP,K2_POP,W0_POP,LAMBDA0, OMEGA_LAMBDA0, OMEGA_LAMBDA1, CV",
         #    thinupdater=1,
          #  thinsamples=5, # uno dei due
            odesolver="RK45", #oppure LSODA
        #    dic="true",
        	winbugsgui="true" 
        }
    }
} 	

simeoni2004 =mogObj{
	OBJECTS{
	simeoni2004_data : {type is dataObj}
	simeoni2004_pri : {type is priorObj}
	simeoni2004_mdl : {type is mdlObj}
	simeoni2004_BUGS_task : {type is taskObj}
	}
}
