;model definition
<MODEL>	
[LONGITUDINAL]
input={V,kep,kon,KD,kin,R0,kel, k12,k21, Vm, Km}

PK:
; bolus administrations
depot(target=Lfull, p=1/V)
depot(target=LtotQE, p=1/V)
depot(target=LR0, p=1/V)
depot(target=LIB, p=1/V)
depot(target=LtotW, p=1/V)
depot(target=LMM, p=1/V)
depot(target=LR0IB, p=1/V)

EQUATION:
odeType=stiff

; parameter transformations
kout = kin/R0
koff=KD*kon
;Vm = kep*R0  
;Km = KD

;===== full ==============================
t_0    = 0 
Lfull_0    = 0 
Rfull_0    = R0 
Pfull_0    = 0 
Afull_0    = 0

ddt_Rfull    = kin-kout*Rfull-kon*Lfull*Rfull+koff*Pfull
ddt_Pfull    = kon*Lfull*Rfull-koff*Pfull-kep*Pfull
ddt_Lfull    = -kel*Lfull-kon*Lfull*Rfull+koff*Pfull -k12*Lfull + k21*Afull/V
ddt_Afull    = k12*Lfull*V - k21*Afull

L_full = Lfull
P_full = Pfull
R_full = Rfull
Ltot_full = Lfull+Pfull
Rtot_full = Rfull+ Pfull

;======= QE ============================
LtotQE_0 = 0 
RtotQE_0 = R0 
AQE_0 = 0

LQE  = 1/2*((LtotQE-RtotQE-KD)+sqrt((LtotQE-RtotQE-KD)^2+4*KD*LtotQE))
ddt_RtotQE = kin - (kep-kout)*(LtotQE-LQE) - kout*RtotQE
ddt_LtotQE = -kep*LtotQE - (kel-kep)*LQE -k12*LQE + k21*AQE/V
ddt_AQE    = k12*LQE*V - k21*AQE                                          

L_QEQSS = LQE
R_QEQSS = RtotQE - (LtotQE - LQE)  
P_QEQSS  = LtotQE - LQE
Ltot_QEQSS = LtotQE
Rtot_QEQSS = RtotQE

;========= constant Rtot =========================
LR0_0 = 0
PR0_0 = 0
AR0_0 = 0

ddt_LR0 = -(kon*R0 + kel)*LR0 + (kon*LR0 + koff)*PR0 -k12*LR0 + k21*AR0/V
ddt_AR0 = k12*LR0*V - k21*AR0
ddt_PR0 = kon*R0*LR0 - (kon*LR0 + koff + kep)*PR0

L_constRtot = LR0
R_constRtot = R0 - PR0
P_constRtot = PR0
Ltot_constRtot = LR0 + PR0
Rtot_constRtot = R0


;============ irreversible binding IB ===================
LIB_0    = 0 
RIB_0    = R0 
PIB_0    = 0 
AIB_0    = 0

ddt_RIB    = kin-kout*RIB-kon*LIB*RIB
ddt_PIB    = kon*LIB*RIB-kep*PIB
ddt_LIB    = -kel*LIB-kon*LIB*RIB -k12*LIB + k21*AIB/V
ddt_AIB    = k12*LIB*V - k21*AIB

L_IB = max(LIB,1e-6)
R_IB = RIB
P_IB = PIB
Ltot_IB = max(LIB,1e-6) + PIB
Rtot_IB = RIB + PIB

;==========const Rtot +  IB===========================
LR0IB_0 = 0
PR0IB_0 = 0
AR0IB_0 = 0

ddt_LR0IB = -(kon*R0 + kel)*LR0IB + (kon*LR0IB )*PR0IB -k12*LR0IB + k21*AR0IB/V
ddt_AR0IB = k12*LR0IB*V - k21*AR0IB
ddt_PR0IB = kon*R0*LR0IB - (kon*LR0IB + kep)*PR0IB

L_constRtotIB = max(LR0IB,1e-6)
R_constRtotIB = R0 - PR0IB
P_constRtotIB = PR0IB
Ltot_constRtotIB = max(LR0IB,1e-6) + PR0IB
Rtot_constRtotIB = R0

;========= Wagner ==========================
LtotW_0=0
AW_0 = 0

LW = 1/2*((LtotW-R0-KD)+sqrt((LtotW-R0-KD)^2+4*KD*LtotW))
ddt_LtotW = -kep*LtotW - (kel - kep)*LW -k12*LW + k21*AW/V
ddt_AW    = k12*LW*V - k21*AW

L_Wagner = LW
R_Wagner = R0 - (LtotW - LW)
P_Wagner = LtotW - LW
Ltot_Wagner = LtotW
Rtot_Wagner = R0

;======= MM ===============================
LMM_0 = 0
AMM_0 = 0

ddt_LMM = -kel*LMM - Vm*LMM/(Km+LMM) -k12*LMM + k21*AMM/V
ddt_AMM    = k12*LMM*V - k21*AMM

L_MM = max(LMM, 1e-6)


<PARAMETER>

kel=0.0019
V = 0.9
kon = 0.34
KD = 0.00076
kep = 0.0007
kin = 0.22
R0 = 400
k12 = 0.03
k21 = 0.0087
Vm = 0.6
Km = 0.2


 <DESIGN>
[ADMINISTRATION]
;adm100  = {time=0, amount=100}
adm500  = {time=0, amount=500}
;adm1000  = {time=0, amount=1000}

<OUTPUT>
list={L_full,L_QEQSS,L_constRtot,L_IB,L_Wagner,L_constRtotIB,L_MM}
;list={Ltot_full,Ltot_QEQSS,Ltot_constRtot,Ltot_IB,Ltot_Wagner,Ltot_constRtotIB}
;list={R_full,R_QEQSS,R_constRtot,R_IB,R_Wagner,R_constRtotIB}
;list={Rtot_full,Rtot_QEQSS,Rtot_constRtot,Rtot_IB,Rtot_Wagner,Rtot_constRtotIB}
;list={P_full,P_QEQSS,P_constRtot,P_IB,P_Wagner,P_constRtotIB}
grid=0.001:0.1:1200

<RESULTS>
[GRAPHICS]
p1 = {y={L_full,L_QEQSS,L_constRtot,L_IB,L_Wagner,L_constRtotIB,L_MM}, ylabel='Free ligand concentration L', xlabel='Time (hr)'}
;p2 = {y={Ltot_full,Ltot_QEQSS,Ltot_constRtot,Ltot_IB,Ltot_Wagner,Ltot_constRtotIB}, ylabel='Total ligand concentration Ltot', xlabel='Time (hr)'}
;p3 = {y={R_full,R_QEQSS,R_constRtot,R_IB,R_Wagner,R_constRtotIB}, ylabel='Free receptor concentration R', xlabel='Time (hr)'}
;p4 = {y={Rtot_full,Rtot_QEQSS,Rtot_constRtot,Rtot_IB,Rtot_Wagner,Rtot_constRtotIB}, ylabel='Total receptor concentration Rtot', xlabel='Time (hr)'}
;p5 = {y={P_full,P_QEQSS,P_constRtot,P_IB,P_Wagner,P_constRtotIB}, ylabel='Complex concentration P', xlabel='Time (hr)'}