SIM_example.mod

;; 2. Description: EXAMPLE SIMULATION SCRIPT

;; x1. Author: DBi

$PROBLEM    Simulate ACTH-CORT using Individual parameters (with/without HC)

$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 ISA1 IBASE IEC50 IBMAX IIEMAX IKOUT ICL IV3 ETC1	
            ETC2 ETC3 ETC4 ETC5 ETC6 ETC7 ETC8 ETC12 ETC15 DOSEN DOSE 

$DATA       XXX.csv IGNORE=@ 
$SUBROUTINE ADVAN13 TOL=12
$MODEL      NCOMP=6 COMP=(DEPOT) COMP=(ACTH) COMP=(CORT) COMP=(CORTP) COMP=(AUC_CORT) COMP=(AUC_ACTH)
$ABBR	    COMRES=6
$PK 
"FIRST
" USE PRDATA, ONLY: MXSTP01
" MXSTP01=2147483647 
  
;------------------------ IR HC ABSORPTION (TCAM) ------------------
IF(EVID.EQ.3) THEN
COM(1)=0
COM(2)=0
COM(3)=0
COM(4)=0
COM(5)=0
COM(6)=0
ENDIF
IF(NEWIND.LE.1) THEN
COM(1) = TIME  
COM(2) = 0
COM(3) = TIME
COM(4) = 0
COM(5) = TIME
COM(6) = 0
ENDIF
IF(EVID.EQ.1.AND.DOSEN.EQ.1) COM(1) = TIME 
TEMP=COM(1)
IF(EVID.EQ.1.AND.DOSEN.EQ.2) COM(3) = TIME 
TEMP2=COM(3)
IF(EVID.EQ.1.AND.DOSEN.EQ.3) COM(5) = TIME 
TEMP3=COM(5)
IF(EVID.EQ.1.AND.DOSEN.EQ.1) COM(2) = AMT 
IF(EVID.EQ.1.AND.DOSEN.EQ.2) COM(4) = AMT 
IF(EVID.EQ.1.AND.DOSEN.EQ.3) COM(6) = AMT 
IF(FLAG.EQ.10) THEN
COM(2)=0
COM(4)=0
COM(6)=0
ENDIF
DOSTC = COM(2)  
DOSTC2 = COM(4)
DOSTC3 = COM(6)

KA=THETA(25)*EXP(ETA(9))

VAR=0 
IF(FLAG.EQ.1) THEN
VAR=ETA(18)
ENDIF
IF(FLAG.EQ.2) THEN
VAR=ETA(19)
ENDIF
IF(FLAG.EQ.3) THEN
VAR=ETA(20)
ENDIF
IF(FLAG.EQ.4) THEN
VAR=ETA(21)
ENDIF
IF(FLAG.EQ.6) THEN
VAR=ETA(22)
ENDIF

IF(DOSE.NE.0.AND.FLAG.NE.10)THEN
DOSECOV=DOSE
ENDIF

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))
SA1=ISA1*EXP(ETA(1)*SQRT(ETC1))
SW1=THETA(2)
Pt=THETA(3)

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

;---------------- INHIBITION OF ACTH ------------------
;KOUT=THETA(6)*EXP(ETA(2))
KOUT=IKOUT*EXP(ETA(2)*SQRT(ETC2)) 
IMAX=THETA(7)
IC50=THETA(8)
GAMMAI=THETA(22)
IMAX_BASE=THETA(23)

;BASE=THETA(9)*EXP(ETA(7))
BASE=IBASE*EXP(ETA(7)*SQRT(ETC7)) 

;---------------- CORT SECRETION (PD) -----------------
; IEMAX=Scaling factor for EMAX, IIEMAX=Individual IEMAX from modelling 

IF(IIEMAX.NE.1) LOGITIEMAX=LOG(IIEMAX/(1-IIEMAX))

IF(GEN_ACT_FACTOR.EQ.8) LOGITIEMAX=LOG(IIEMAX/(1-IIEMAX))+(ETA(4)*SQRT(ETC4))

IF(IIEMAX.NE.1) IEMAX=EXP(LOGITIEMAX)/(1+EXP(LOGITIEMAX)) 
IF(IIEMAX.EQ.1) IEMAX=IIEMAX

EMAX =THETA(10)*IEMAX
;EC50  =THETA(11)*EXP(ETA(8))
EC50=IEC50*EXP(ETA(8)*SQRT(ETC8)) 
GAMMAE=THETA(12)

KIN  =BASE*KOUT 

;-------------- CORT DISPOSITION ----------------------
CL = ICL*EXP(ETA(12)*SQRT(ETC12))
;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*EXP(ETA(15)*SQRT(ETC15))
;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)
BMAX=IBMAX*EXP(ETA(5)*SQRT(ETC5))
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_C=THETA(35)*EXP(ETA(23))
A_0(3)=IBASE_C

;-------------- 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)

TSTAR = TIME - TEMP
IF(TSTAR.LE.0) TSTAR=0.0001
IF(DOSTC.GT.0) THEN
INP = EXP(LOG(DOSTC*BIO+.00001) + LOG(KTR) + NN*LOG(KTR*TSTAR+.00001) - KTR*TSTAR - LNFAC) 
ELSE
INP = 0
ENDIF  
TSTAR2 = TIME - TEMP2
IF(TSTAR2.LE.0) TSTAR2=0.0001
IF(DOSTC2.GT.0) THEN
INP2 = EXP(LOG(DOSTC2*BIO+.00001) + LOG(KTR) + NN*LOG(KTR*TSTAR2+.00001) - KTR*TSTAR2 - LNFAC) 
ELSE
INP2 = 0
ENDIF 
TSTAR3 = TIME - TEMP3
IF(TSTAR3.LE.0) TSTAR3=0.0001
IF(DOSTC3.GT.0) THEN
INP3 = EXP(LOG(DOSTC3*BIO+.00001) + LOG(KTR) + NN*LOG(KTR*TSTAR3+.00001) - KTR*TSTAR3 - LNFAC)  
ELSE
INP3 = 0
ENDIF
TINP=INP+INP2+INP3

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)=TINP-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)+KA*A(1)
DADT(4)=-k32*A(4)+k23*AU

DADT(5)=A(3)
DADT(6)=A(2)

AUC_CORT=A(5)
AUC_ACTH=A(6)


$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) FIX ; PROP_ERR_ACTH
(0) FIX ; ADD_ERR_CORT
(0) FIX ; PROP_ERR_CORT
(0) FIX ; ADD_ERR_UNBCORT
(0) FIX ; PROP_ERR_UNBCORT
(0, 2.38) FIX ; SW2
(0) 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

$OMEGA
 1 FIX  ; IIV_SA-SW
 1 FIX  ; IIV_Kout
 0 FIX ; IIV_P453S
 1 FIX ; IIV_V281L
 1 FIX  ; IIV_NCBG
 0 FIX  ; IIV_P30L
 1 FIX ;  IIV_BASE
 1 FIX ;  IIV_EC50 
 0 FIX  ; IIV_KA
 0.166 FIX  ; IIV_NN
 0.215 FIX  ; IIV_F 
 1 FIX  ; IIV_CL
 0 FIX  ; IIV_V2
 0 FIX  ; IIV_Q
 1 FIX  ; IIV_V3
 0 FIX  ; IIV_KD
 0 FIX  ; IIV_NS   
$OMEGA BLOCK(1) 0.0792 FIX  ; IOV_MTT
$OMEGA BLOCK (1) SAME ; VAR 2
$OMEGA BLOCK (1) SAME ; VAR 5
$OMEGA BLOCK (1) SAME ; VAR 10
$OMEGA BLOCK (1) SAME ; VAR 20 
$OMEGA
 0 FIX  ; IIV_Base 

$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 GEN_ACT_FACTOR EC50 IC50 PRED IPRED IEMAX AUC_CORT AUC_ACTH NOPRINT NOAPPEND ONEHEADER FILE=simtabXXX