$PROB  ETAS WERE ADDED TO ALL THE PK PR AND RET0, RBC0, KSAI, KD, MCH ERR2 FIXED TO 0
$INPUT      ID AMT TIME CMT DV DRUG ROUT RATE 
$DATA Simulated_EPO_dataset.csv IGNORE=I IGNORE=(DRUG=1)

$SUBROUTINES  ADVAN6 TOL=6 
$EST  MAXEVALS=9999 SIGDIGITS=2 PRINT=1 METHOD=1 INTERACTION NOABORT MSFO=EPO_FINAL_simulated_dataset.msf
$COV PRINT=E

;*** PK PARAMETERS THETAS ***
$THETA (0,2.4,10)   ; D1
$THETA (0,0.045,5)  ; KA 
$THETA (0,4.9,10)   ; VC
$THETA (0,0.3,0.6)  ; CL
$THETA (0,0.33,1)   ; F1
$THETA 0 FIXED      ; ALAG1
$THETA (0,8.2,20)   ; C0
$THETA (0,0.06,1)   ; K23
$THETA (0,0.1,1)    ; K32

;*** PD PARAMETERS ***
$THETA (0,54,100) ; RET0  
$THETA 0 FIXED    ; TD3   
$THETA (0,45,200) ; TP1 
$THETA 0 FIXED    ; TP2  
$THETA 0 FIXED    ; TP3 DUMMY  
$THETA (0,3,40)   ; SC50  
$THETA (0,0.04,1) ; SMAX
$THETA (50,700,2000) ; KD
$THETA (0,0.8)    ; KSAI
$THETA (0,4,20)   ; RBO
$THETA (0,34,100) ; MCH

$OMEGA 0.5  ; ETAD1
$OMEGA 0.1  ; ETAKA
$OMEGA 0.09 ; ETAVC
$OMEGA 0.5  ; ETACL
$OMEGA 0.1  ; ETAC0
$OMEGA 0.05 ; ETARET0
$OMEGA 0 FIXED ; TP1
$OMEGA 0 FIXED ; SC50
$OMEGA 0.8  ; ETASMAX
$OMEGA 0 FIXED ; ETAKD
$OMEGA 0.2  ; ETAKSAI
$OMEGA 0.01 ;ETARBC0    
$OMEGA 0.01 ;ETAMCH 

$SIGMA 0.04; ERRSD1 
$SIGMA 0 FIXED; ERRSD2
$SIGMA 0.04; ERRSD3 
$SIGMA 0.01; ERRSD4


$MODEL  NCOMP=9
COM  = (A1)              ;COMPARTMENT #1 EPO SC
COM  = (A2)              ;COMPARTMENT #2 EPO SERUM
COM  = (A3)              ;COMPARTMENT #3 EPO TISSUE
COM  = (A4)              ;COMPARTMENT #4 PRECURSOR 1
COM  = (A5)              ;COMPARTMENT #5 PRECURSOR 2
COM  = (A6)              ;COMPARTMENT #6 PRECURSOR 3
COM  = (A7)              ;COMPARTMENT #7 RET       
COM = (A8)               ;COMPARTMENT #8 RBC
COM = (A9)               ;COMPARTMENT #9 HGB
$PK
;*** PK PARAMETERS ***  
D1 = THETA(1)*EXP(ETA(1))         ;ZERO ORDER INFUSION DURATION
KA = THETA(2)*EXP(ETA(2))         ;FIRST ORDER ABSORPTION RATE   
VC = THETA(3)*EXP(ETA(3))         ;VOLUME OF DISTRIBUTION OF SERUM COMPARTMENT
CL = THETA(4)*EXP(ETA(4))         ; LINEAR CLEARANCE
F1 = THETA(5)                     ; BIOAVAILABILITY 
ALAG1 = THETA(6)                  ; SHOULD FIXED TO ZERO
C0 = THETA(7)*EXP(ETA(5))         ; EPO BASELINE
K23 = THETA(8)                    ;TISSUE DISTRIBUTION CONSTANT
K32 = THETA(9)                    ;TISSUE DISTRIBUTION CONSTANT
   
;*** PD PARAMETERS ***
RR0   =  THETA(10)*EXP(ETA(6))   ; RET0                      
TD3    =  THETA(11)              ; TR
TP1   =  THETA(12)*EXP(ETA(7))   ; TP1
TD2   =  THETA(13)               ; TP2
TD3   =  THETA(14)               ; TP3 DUMMY
SC50  =  THETA(15)*EXP(ETA(8))   ; SC50
SMAX  =  THETA(16)*EXP(ETA(9))   ; SMAX
KD    = THETA(17)*EXP(ETA(10))   ; KD 
KSAI  = THETA(18)*EXP(ETA(11))   ; PROPORTIONAL COEFFICIENT   
RB0 =   THETA(19)*EXP(ETA(12))   ; RBC BASELINE 
MCH =   THETA(20)*EXP(ETA(13))   ; COEFFICIENT FOR HGB
TR =TP1
TP3=TP1
TP2=TP1
TB=(RB0-RR0/1000)*TR*1000/RR0     ;RBC LIFESPAN
;*** INITIAL CONDITION FOR PK ***
A_0(1) = 0                        ;INITIAL CONDITION FOR SC COMPARTMENT 
A_0(2) = C0*VC                    ;INITIAL CONDITION FOR SERUM COMPARTMENT
A_0(3) = C0*VC*K23/K32            ;INITIAL CONDITION FOR TISSUE COMPARTMET
   
;*** INITIAL CONDITION FOR PD ***
RTOT0 = RR0*TP2/TR
RC0 = RTOT0*C0/(KD+C0)
SC0 = SMAX*RC0/(SC50+RC0)
A_0(4) = RR0*TP1/(TR*SC0)*(1/2**5)
A_0(5) = RR0*TP2/TR
A_0(6) = RR0*TP3/TR
A_0(7) = RR0
A_0(8) = RB0-RR0/1000
 
$DES
;*** EQUATIONS FOR PK ***
RTOT = A(5)
KEPO = CL*C0+RTOT0*KSAI*VC*C0/(KD+C0)  ; ZERO ORDER PRODUCTION OF ENDOGENOUS EPO    
 
DADT(1) =     -KA*A(1)                                                                           ;EPO SC
DADT(2) = KEPO+KA*A(1)-CL*A(2)/VC-K23*A(2)+K32*A(3)-RTOT*KSAI*VC*(A(2)/VC)/(KD+A(2)/VC)          ;EPO SERUM
DADT(3) =                         K23*A(2)-K32*A(3)                                              ;EPO TISSUE

;*** EQUATIONS FOR PD ***
C2 = A(2)/VC
RC = RTOT*C2/(KD+C2)
SC = SMAX*RC/(SC50+RC)
KIN0  =  RR0/TR*(1/2**5)
DADT(4) =  KIN0-1/TP1*SC*A(4)
DADT(5) = 2**5*(1/TP1)*SC*A(4)-(1/TP2)*A(5)
DADT(6) =                      (1/TP2)*A(5)-(1/TP3)*A(6)
DADT(7) =                                   (1/TP3)*A(6)-A(7)/TR
DADT(8) =                                                A(7)/(TR*1000)-A(8)/TB
 
$ERROR  
IF (CMT.EQ.2) THEN
IPRED=A(2)/VC
Y = IPRED*(1 + EPS(1))+EPS(2)
ELSE
ENDIF
      
IF (CMT.EQ.7) THEN
IPRED=A(7)
Y = IPRED*(1+EPS(3))
ELSE
ENDIF

IF (CMT.EQ.8) THEN                   ;RBC = A(8)+A(7)/1000
IPRED=A(8)+A(7)/1000
Y = IPRED*(1+EPS(4)) 
ELSE
ENDIF
      
IF (CMT.EQ.9) THEN                   ;HGB= MCH*RBC/10
IPRED=MCH*(A(8)+A(7)/1000)/10
Y = IPRED*(1+EPS(4)) 
ELSE
ENDIF