### Model for comparison with IVGTT Minimal model

### Solve first with fixglu=1 to steady-state, and use this state as initial condition for your simulation such as in Fig. 2 

# IVGTT minimal model approximation
X'=Y-m*X
Y'=alpha*(Yinf-Y)

###

# comment out if you don't have IVGTT glucose data - you can use the exponential G description as surrogate
table GGG  gluXPP_IVGTT.dat

par fixglu=0, glu0=4.9127222
G=(1-fixglu)*GGG(t)+fixglu*18*glu0
#G=(1-fixglu)*18*(glu0+(17.2-glu0)*exp(-t/30))+fixglu*18*glu0
aux glu=G





% [1..500]
h[j]'=pplus*I*phi([j])-pminus*h[j]-fplus*h[j]*Heav(G-[j])
%
F'=fplus*(sum(0,G-1)of(shift(h1,i'))) -m*F

aux actRRP=(sum(0,G-1)of(shift(h1,i')))



### 43 changes ug/pancreas/4L blood to pM [= 1e6/(4*5808) ]

par SecrBase=0.58
secr= m*F + SecrBase
aux SR=secr


I'=Mob - r*I - pplus*I +pminus*(sum(0,499)of(shift(h1,i')))
I(0)=311

Mob=M0+c*(G^nM)/(KmM^nM+G^nM)
aux m=mob

phi(ggg)=nH*max(0,ggg-90)^(nH-1)*(KmH-90)^nH/((KmH-90)^nH+max(0,ggg-90)^nH)^2
par nH=3, KmH=130


!ppp=phi(200)


# Fitted to Grodsky 1972 Fig 10 (KmM, nM), and adjusted to reproduce C-peptide (c, p+, p-)
par tau=15, KmM=180, nM=4, c=200, M0=14, r=0.08


par pminus=0.01, pplus=0.003, k=0, fplus=6.2, m=0.62


aux RRP=(sum(0,499)of(shift(h1,i')))

@ meth=5dp, dt=1, total=100, bound=500000, maxstor=500000
@ yp=f, xp=t
@ xlo=0, xhi=360, ylo=0, yhi=10, 
@ BUT=quit:fq, BUT=export:go

### Add Cpeptide for control
### 43 changes ug/pancreas/4L blood to pM [= 1e6/(4*5808) ]

par k21=0.0559, k12=0.0492, k01=0.06, fudge=0.5

CP1'= fudge* 43*secr +k12*CP2-k21*CP1-k01*CP1
CP2'= k21*CP1-k12*CP2



############ simplifying stuff

hG=sum(0,499)of((if(abs(i'-G+1)<0.5)then(1)else(0))*shift(h1,i'))
aux hG=hG
#aux test=sum(0,499)of((if(abs(i'-G+1)<0.5)then(1)else(0)))


PhiG=max(0,G-90)^nH /( (KmH-90)^nH + max(0,G-90)^nH )

# IVGTT minimal model approximation
#X'=Y-m*X
#Y'=alpha*(Yinf-Y)

alpha = r+pplus*PhiG*fplus/(pminus+fplus)
aux alpha=alpha
Yinf  = fplus*pplus*PhiG*Mob/(pminus+fplus)/alpha
aux Yinf=Yinf

X(0)=92.7
!X0=(pplus/(pminus+f))*311*(max(0,18*17.2 - 90)^nH /( (KmH - 90)^nH + max(0,17.2*18 - 90)^nH ))

DD'=alpha*(mob/alpha-DD)
aux HH=pplus*I*PhiG/(pminus+fplus)
aux RX=pplus*(pminus+(1-PhiG)*fplus)/pminus/(pminus+fplus)*I

done