MOD_example.mod

;; 1. Based on: 
;; 2. Description: EXAMPLE MODEL SCRIPT
;; x1. Author: DBi


$PROBLEM    Estimate individual scaling factor for Emax

$INPUT      ID TIME AMT	RATE MDV LINDV DV EVID CMT FLAG BW BMI HT AGE SEX CBG PHENOTYPE OCC GEN_ACT GEN_ACT_MIN GEN_MM_ACT GEN_ACT_FACTOR OBE ADRENAL ADR_RANK EXCL

$DATA       XXX.csv IGNORE=@ 
$SUBROUTINE ADVAN13 TOL=12
$MODEL      NCOMP=4 COMP=(DEPOT) COMP=(ACTH) COMP=(CORT) COMP=(CORTP) 
$PK 
"FIRST
" USE PRDATA, ONLY: MXSTP01
" MXSTP01=2147483647 
 
;------------------------ ABSORPTION ------------------
;KA=THETA(25)*EXP(ETA(9))

;MTT=THETA(26)*(DOSECOV/5)**THETA(36)*EXP(VAR)
;NN=THETA(34)*EXP(ETA(10))
;KTR=(NN+1)/MTT

;F1=0
;TVBIO=THETA(27)
;LOGITBIO=LOG(TVBIO/(1-TVBIO))
;BIO=EXP(LOGITBIO+ETA(11))/(1+EXP(LOGITBIO+ETA(11)))    

;LNFAC=LOG(2.5066)+(NN+0.5)*LOG(NN)-NN  
;NFAC =SQRT(2*3.1415)*NN**(NN+0.5)*EXP(-NN) 

;----------------- ACTH SURGE FUNCTIONS ---------------
SA1=THETA(1)*EXP(ETA(1))
SW1=THETA(2)
Pt=THETA(3)

SA2=THETA(4)
SW2=THETA(20)
Pt2=THETA(5)

KOUT=THETA(6)*EXP(ETA(2))
BASE=THETA(9)*EXP(ETA(7))
KIN  =BASE*KOUT 

;---------------- INHIBITION OF ACTH ------------------
IMAX=THETA(7)
IC50=THETA(8)
GAMMAI=THETA(22)
IMAX_BASE=THETA(23)

;---------------- CORT SECRETION (PD) -----------------
IF(GEN_ACT_FACTOR.EQ.1) THEN
LOGITIEMAX=THETA(38)
ENDIF
IF(GEN_ACT_FACTOR.EQ.2) THEN
LOGITIEMAX=THETA(39)
ENDIF
IF(GEN_ACT_FACTOR.EQ.3) THEN
LOGITIEMAX=THETA(40)
ENDIF
IF(GEN_ACT_FACTOR.EQ.4) THEN
LOGITIEMAX=THETA(41)
ENDIF
IF(GEN_ACT_FACTOR.EQ.5) THEN
LOGITIEMAX=THETA(42)
ENDIF 
IF(GEN_ACT_FACTOR.EQ.6) THEN
LOGITIEMAX=THETA(43)+ETA(6)
ENDIF
IF(GEN_ACT_FACTOR.EQ.7) THEN
LOGITIEMAX=THETA(44)+ETA(3)
ENDIF
IF(GEN_ACT_FACTOR.EQ.8) THEN
LOGITIEMAX=THETA(45)+ETA(4)
ENDIF
IF(GEN_ACT_FACTOR.EQ.9) THEN
LOGITIEMAX=THETA(46)
ENDIF 

IF(AGE.LE.13) THEN
IEMAX=EXP(LOGITIEMAX)/(1+EXP(LOGITIEMAX))*THETA(47)
ELSE
IEMAX=EXP(LOGITIEMAX)/(1+EXP(LOGITIEMAX))
ENDIF

IF(GEN_ACT_FACTOR.EQ.0) IEMAX=1
EMAX =THETA(10)*IEMAX
EC50  =THETA(11)*EXP(ETA(8))
GAMMAE=THETA(12)

;-------------- CORT DISPOSITION ----------------------
;CL = ICL
CL = THETA(28)*(BW/70)**0.75*EXP(ETA(12))
;V2= IV2
V2 = THETA(29)*(BW/70)*EXP(ETA(13))
;Q = IQ
Q = THETA(30)*(BW/70)**0.75*EXP(ETA(14))
;V3= IV3
V3 = THETA(31)*(BW/70)*EXP(ETA(15))

k20=CL/V2
k23=Q/V2
k32=Q/V3

;--------------- BINDING MODEL ------------------------
;CBGmol=CBG*1000/52 ; Converting CBG from ug/mL to nmol/L. Mw=52000g/mol
CBG2=CBG            ;CBG given in dataset for study 2.
IF(CBG.EQ.-99) THEN
CBG2=22.4 ; CBG is sampled from the distribution for study 1.
ENDIF
CBGmol=CBG2*1000/52

BMAX=CBGmol*EXP(LOG(THETA(13))+ETA(5)) ; Bmax=CBG concentration * number of binding sites at CBG (fixed to 1)
AMAX=BMAX*V2

;KD= IKD
KD = THETA(32)*EXP(ETA(16))
K=KD*V2

;NS= INS ; Nonspecific binding parameter
NS = THETA(33)*EXP(ETA(17))

; ------------- INITIALISE CMTs -----------------------
A_0(2)=BASE
IBASE=THETA(35)
A_0(3)=IBASE

;-------------- DEX Inihibition -----------------------
IF(FLAG.EQ.20) THEN
IDEX=1
IDEX_BASE=1
ENDIF

IF(FLAG.NE.20)THEN
IDEX=1-(THETA(21)/10000)
IDEX_BASE=1-(THETA(24)/10000)
ENDIF

;-------------- DAY - repeat surges -------------------
DAY=-1
IF(TIME.GE.9.AND.TIME.LT.33) THEN
DAY=0
ENDIF
IF(TIME.GE.33.AND.TIME.LT.57) THEN 
DAY=1
ENDIF
IF(TIME.GE.57.AND.TIME.LT.81) THEN
DAY=2
ENDIF 

;--------------DIFFERENTIAL EQUATIONS------------------
$DES 
SUR1=SA1/((ABS(T-(Pt+(24*DAY)))/SW1)**4+1)
SUR2=SA2/((ABS(T-(Pt2+(24*DAY)))/SW2)**4+1)

AU = (A(3)-K*(1+NS)-AMAX+SQRT(((A(3)-K*(1+NS)-AMAX)**2)+4*K*A(3)*(1+NS)))/(2*(1+NS)) 
FL=1-((IMAX/10000)*((AU/V2)**GAMMAI)/((AU/V2)**GAMMAI+IC50**GAMMAI)) ; Feedback loop 

;DADT(1)=EXP(LOG(BIO*DOSETRR+.00001)+LOG(KTR)+NN*LOG(KTR*TAD+.00001)-KTR*TAD-LNFAC)-KA*A(1) 
DADT(2)=KIN+((SUR1+SUR2)*IDEX*FL)-KOUT*A(2)

EACTH=EMAX*A(2)**GAMMAE/(EC50**GAMMAE+A(2)**GAMMAE)
DADT(3)=EACTH-(k20+k23)*AU+k32*A(4)
DADT(4)=-k32*A(4)+k23*AU

$ERROR
AU2 = (A(3)-K*(1+NS)-AMAX+SQRT(((A(3)-K*(1+NS)-AMAX)**2)+4*K*A(3)*(1+NS)))/(2*(1+NS)) 

IF(CMT.EQ.2) THEN
IPRED=A(2)
W=SQRT(THETA(15)**2+(THETA(14)**2/IPRED**2))
ENDIF
IF(CMT.EQ.3) THEN
IPRED=A(3)/V2
W=SQRT(THETA(17)**2+(THETA(16)**2/IPRED**2))
ENDIF

IF(IPRED.GT.0) THEN 
IPRED=LOG(IPRED)
ELSE
IPRED=LOG(IPRED+0.01)
ENDIF

IRES=DV-IPRED
IWRES=IRES/W
Y=IPRED+W*EPS(1)

$THETA
(0, 357) FIX ; SA1
(1.03) FIX ; SW
(0, 14.9,24) FIX ; Pt
(0, 54.3) FIX ; SA2
(0, 21,24) FIX ; Pt2
(0, 0.753) FIX ; Kout
(9800) FIX ; IMAX
(0, 4.51) FIX ; IC50
(0, 1.41) FIX ; BASE
(0, 4010) FIX ; EMAX
(0, 5.49) FIX ; EC50
(3.3) FIX ; GAMMAeff
(1) FIX ; NCBG
(0) FIX ; ADD_ERR_ACTH
(0, 0.26) FIX ; PROP_ERR_ACTH
(0) FIX ; ADD_ERR_CORT
(0, 0.146) FIX ; PROP_ERR_CORT
(0) FIX ; ADD_ERR_UNBCORT
(0) FIX ; PROP_ERR_UNBCORT
(0, 2.38) FIX ; SW2
(10000) FIX ; IDEX
(0, 7.09) FIX ; GAMMAinh
(0) FIX ; IMAX_BASE
(0) FIX ; IDEX_BASE
(0, 28.5) FIX ; KA [nmol/h]
(0, 0.865) FIX ; MTT
(0, 0.35,1) FIX ; F1 [-]
(0, 108) FIX ; 1. CL [L/h]
(0.001, 2.22) FIX ; 2. V2 [L]
(0, 93.9) FIX ; 3. Q [L/h]
(0.01, 63.2) FIX ; 4. V3 [L]
(9.71) FIX ; 9. KD [nmol/L],fixed from binding model
(4.15) FIX ; 10. NS [-],fixed from binding model
(0, 2.15) FIX ; NN
(0, 15.7) FIX ; Cort base
(0.179) FIX ; expMTTdose
(0) FIX ; SA1BW  
(-2.23) ; NULL
(-2.06) ; In2G
(-1.38) ; I77T
(-1.86) ; I172N
(0.822) FIX ; R124C
(-0.166) FIX ; P30L
(1.3) FIX ; P453S
(1.5) FIX ; V281L 
(0.0358) FIX ; P482S 
(0.316) ; AGE-IEMAX

$OMEGA
 0.636 FIX  ; IIV_SA-SW
 0.288 FIX  ; IIV_Kout
 0 FIX  ; IIV_P453S
 0.787 ; IIV_V281L
 0.049 FIX  ; IIV_NCBG
 0 FIX  ; IIV_P30L
$OMEGA BLOCK(2)
 0.0579 FIX ;  IIV_BASE
 0.0618 0.072 ;  IIV_EC50 
$OMEGA
 0 FIX  ; IIV_KA
 0 FIX  ; IIV_NN
 0 FIX  ; IIV_F 
 0.0133 FIX  ; IIV_CL
 0 FIX  ; IIV_V2
 0 FIX  ; IIV_Q
 0.0151 FIX  ; IIV_V3
 0 FIX  ; IIV_KD
 0 FIX  ; IIV_NS   

$SIGMA 1 FIX  ; SIGMA    
;$SIMULATION (12345) (54321) ONLYSIM NSUB=1000

$ESTIMATION MAXEVAL=9999 PRINT=1 METH=1 INTER NOABORT SIGDIG=3
$COVARIANCE FPOSDEF=1 MATRIX=R PRINT=E UNCONDITIONAL

$TABLE ID TIME MDV DV EVID CMT EC50 IC50 PRED IPRED CWRES IEMAX SA1 BASE KOUT CL V2 V3 Q BMAX ETA1 ETA2 ETA3 ETA4 ETA5 ETA6 ETA7 ETA8 ETA12 ETA15 LINDV NOPRINT NOAPPEND ONEHEADER FILE=sdtabXXX