INPUT:
parameter={V_b, V_l, ke_C, NB0, kout_NB, LF_B, alpha_S, B_max, eps_N, K_Dn, NL0, K_NnB, K_NnM, S_NL,MB0, kout_MB, r_max, d_M, AML0, S_AML, kout_AM,Emax_clo,CL50,K_Nm,I_k_AML, w_R, LF_P, phi, P_C,ecl, burst, decay, psi, iota_M, gam_M, M_50}

PK:
depot(type=2, target=S)
depot(type=3, target=P)
depot(type=4, target=R)
depot(type=5, target=CL)


EQUATION:
odeType = stiff
t0=0

; Parameter transformation
Rin_NB = NB0 * kout_NB
MoM0 = AML0/3
Rin_MB = MB0 * kout_MB

k_Mbl=kout_MB*2 ;no idea
k_Mlb=k_Mbl*MB0/MoM0*V_b/V_l

AM_max = AML0*3 ;data ~baseline (PAKP1 treated) / Koogushi: Baseline * 2 /Misharin2018 baseline *2-3 
alpha_AM = kout_AM / (1-AML0/AM_max)
K_rm=K_NnB 

alpha_R=w_R*alpha_S
B=S+I1+I2+R


;====== PD part of the model
; Initial values
CL_0=0
NB_0 = NB0
NL_0 = NL0
MB_0 = MB0
ML0_0 = MoM0
ML1_0 = 0
ML2_0 = 0
ML3_0 = 0
ML4_0 = 0
ML5_0 = 0
AML_0 = AML0

S_0 = 0
I1_0 = 0
I2_0 = 0
R_0 = 0
P_0 = 0


                               
                               
; ODEs
ddt_S = alpha_S*S*(1-B/B_max) - ((eps_N*NL*S)/(1+B/K_Dn)) - phi*(1-AML^gam_M*iota_M/(AML^gam_M+M_50^gam_M))*P/(1+P/P_C)*S
ddt_R = alpha_R*R*(1-B/B_max) - ((eps_N*NL*R)/(1+B/K_Dn))
ddt_I1 = phi*(1-AML^gam_M*iota_M/(AML^gam_M+M_50^gam_M))*S*P/(1+P/P_C) - ((eps_N*NL*I1)/(1+B/K_Dn)) - ecl * I1
ddt_I2 = ecl * I1 - ((eps_N*NL*I2)/(1+B/K_Dn)) - ecl * I2

ddt_P =  - decay*P - phi*(1-AML^gam_M*iota_M/(AML^gam_M+M_50^gam_M))*(S+I1+I2)*P/(1+P/P_C) - psi*P*AML + burst*ecl*I2

ddt_NB = Rin_NB - kout_NB*NB - r_max * NB * B/(B+K_NnB) *(1-I_k_AML*(1-AML/(AML+K_NnM)))
ddt_NL =  r_max * NB * B/(B+K_NnB) * (V_b/V_l) *(1-I_k_AML*(1-AML/(AML+K_NnM)))

ddt_MB = Rin_MB - (kout_MB + k_Mbl*(1+(r_max-k_Mbl)/k_Mbl*(B/(B+K_rm)))) * MB + k_Mlb * ML0 * (V_l/V_b)  ;blood monocytes
ddt_ML0= k_Mbl * MB*(1+(r_max-k_Mbl)/k_Mbl*(B/(B+K_rm))) * V_b/V_l - k_Mlb * ML0 - d_M*ML0*(B/(B+K_Nm))*(1-AML/AM_max)  ;lung monocytes 
ddt_ML1=  d_M*(B/(B+K_Nm))*(1-AML/AM_max)  * (ML0 - ML1) ;lung monocytes transit1
ddt_ML2=  d_M*(B/(B+K_Nm))*(1-AML/AM_max)  * (ML1 - ML2) ;lung monocytes transit2
ddt_ML3=  d_M*(B/(B+K_Nm))*(1-AML/AM_max)  * (ML2 - ML3) ;lung monocytes transit3
ddt_ML4=  d_M*(B/(B+K_Nm))*(1-AML/AM_max)  * (ML3 - ML4) ;lung monocytes transit4
ddt_ML5=  d_M*(B/(B+K_Nm))*(1-AML/AM_max)  * (ML4 - ML5) ;lung monocytes transit5
ddt_AML =  alpha_AM*AML * (1-AML/AM_max) - kout_AM * AML * (1+Emax_clo*CL/(CL+CL50))+d_M*ML5*(B/(B+K_Nm))*(1-AML/AM_max) ; macrophages

ddt_CL = -ke_C * CL ;Clodronate amound in the lung


log_B=log10(B/(V_l*1000*LF_B)+1e-6)
log_P=log10(P/(V_l*1000*LF_P)+1e-6)

N_B = NB+1e-6
N_BAL = NL * S_NL
AML_BAL = AML * S_AML



OUTPUT: 
output = {AML_BAL,log_B,MB,NB,N_BAL,log_P}
table = {NL,ML0,AML,CL,S,R,I1,I2}