;; 1. Based on: run20 RW04_07
;; 2. Description: Final PKPD model for isoniazid response in zebrafish larvae infected with M. marinum (E11)
;; x1. Author: R.C. van Wijk (r.c.van.wijk@lacdr.leidenuniv.nl)
;; Author: R.C. van Wijk (r.c.van.wijk@lacdr.leidenuniv.nl)
;; Dataset: VUmc zebrafish larvae constant treated to isoniazid (INH) from 3-5 dpf at 0 and 0.25x-10MIC (MIC = 15 mg/L or 109 uM) 

$PROBLEM PKPD
$INPUT ID TIME AMT DV EVID MDV CMT BQL AQL DOSE AGE AGE_H XEXP BLOOD FLAG
$DATA Real_PKPD_isoniazid_zebrafish.csv IGNORE=@ IGNORE=(BQL.EQ.1) IGNORE=(AQL.EQ.1) 

; time in hours of exposure
; amt in uM (concentration)
; DV in pmole / larvae or pmole / uL (uM) / log10(bac burden)
; Vd in uL
; AGE_H age in hpf 

$SUBROUTINE ADVAN13 TOL=9

$MODEL 
COMP ;CMT 1 dosing
COMP ;CMT 2 isoniazid in larva
COMP ;CMT 3 bacterial burden measured in log10(fluorescence pixels)

$PK

;-----PK-----

TVKE = THETA(1)        ;first order elimination rate
TVV2 = THETA(2)        ;distribution volume in uL/larva
TVKA = THETA(3) / 1e2  ;absorption at start of experiment
KA_GI = THETA(4)       ;effect of GI opening between 3 and 4 dpf on absorption
KA_HPF = THETA(5)      ;effect of age in hours post fertilization on absorption

;age effect on KA
KA = TVKA * KA_HPF**(AGE_H / 101)  ;median age is 101 hpf

;GI effect on KA
IF(AGE_H.LE.90) KIN = KA * DOSE 
IF(AGE_H.GT.90) KIN = KA * DOSE * (1 + KA_GI)

KE = TVKE 
V2 = TVV2

;-----PD-----

KG = THETA(6)                        ;exponential growth rate
INOC = THETA(7) * 1e1                ;inoculum estimated
SLP = (THETA(8) / 1e2) * EXP(ETA(2)) ;slope (linear drug effect)

A_0(3) = INOC * EXP(ETA(1))          ;inoculum set as compartment initial value
	
$DES 
DADT(1) = 0                ;dosing compartment
DADT(2) = KIN - KE * A(2)  ;larva compartment isoniazid

;drug effect

;Distribution volume is estimated based on 5 dpf blood samples, but is expected to be lower with younger larvae of 3 and 4 dpf. 
;It is here assumed that the distribution volume scales with total volume.
IF(AGE.EQ.3) THEN
V = V2 * 253/300 ;taking total volume into account (Guo et al, Biomed. Opt. Express 8, 2611 (2017))
ENDIF
IF(AGE.EQ.4) THEN
V = V2 * 263/300 ;taking total volume into account (Guo et al, Biomed. Opt. Express 8, 2611 (2017))
ENDIF
IF(AGE.EQ.5) THEN
V = V2  ;taking total volume into account (Guo et al, Biomed. Opt. Express 8, 2611 (2017))
ENDIF

C = A(2)/V

EFF = 0 ;for C is zero or NA
IF(C.GT.0) THEN
EFF = SLP * C
ENDIF

DADT(3) = KG * A(3) * (1 - EFF) ;bacterial burden (log10)

$ERROR

;-----PK-----

IF(FLAG.EQ.0.AND.BLOOD.EQ.0) THEN
IPRED = A(2)
Y = IPRED * (1 + EPS(1)) + EPS(2) ; comb error
W = SQRT(IPRED**2*SIGMA(1,1)**2    + SIGMA(2,2)**2)
IRES = DV - IPRED
ENDIF 

IF(FLAG.EQ.0.AND.BLOOD.EQ.1) THEN
IPRED = A(2) / V2
Y = IPRED * (1 + EPS(3)) + EPS(4) ; comb error
W = SQRT(IPRED**2*SIGMA(3,3)**2    + SIGMA(4,4)**2)
IRES = DV - IPRED
ENDIF

;-----PD-----

IF(FLAG.EQ.1) THEN
IPRED = LOG10(A(3)+0.00001)
Y = IPRED + EPS(5)
W = EPS(5)
IRES = DV - IPRED
ENDIF 

IF(W.EQ.0)W=1 
IWRES = IRES/W


$THETA
(0, 0.58)   ; 1 KE 
(0, 0.325)  ; 2 V2
(0, 0.355)  ; 3 KA
(0, 0.178)  ; 4 KA_GI
(1, 7.45)   ; 5 KA_HPF
(0, 0.093)  ; 6 KG 
(0, 1.6)   ; 7 *1e1 INOC
(0, 0.991)  ; 8 /1e2 SLP

$OMEGA 1.64  ;IIV INOC 
$OMEGA 0.227 ;IIV SLP 

$SIGMA 0.0611 ; 1 prop homogenate INH
$SIGMA 0.59   ; 2 add homogenate INH
$SIGMA 0.485  ; 3 prop blood INH
$SIGMA 0 FIX  ; 4 add blood INH
$SIGMA 0.125  ; 5 prop bacterial burden

$ESTIMATION METHOD=1 INTER MAXEVAL=8000 NOABORT PRINT=10 NSIG=3 SIGL=9 POSTHOC
$COVARIANCE PRINT=E

$TABLE ID TIME KIN KA TVKA KE V2 KA_GI KA_HPF KG INOC SLP C V ETA(1) ETA(2) DV EVID MDV CMT BQL AQL DOSE AGE AGE_H BLOOD XEXP FLAG IPRED IRES IWRES CWRES NOPRINT ONEHEADER FILE=tab1