#This script allows to calculate the parameters of a mathematical model for homeostasis of glucose and insulin applied to healthy rats
#It needs glucose and insulin values measured after an oral glucose, the values have to load arranged in a txt file with 3 columns separated by space
#The data should have the following order and the following headings: t G I (where t is the time sample was taken, G is glycemia in mg/dl, I is insulinol in pmol/l)
print("introduce the name of file with data including the extension, example dato.txt")
file<- scan(file = "", what = "",nmax=1)
data<-read.table(file,header=TRUE)
print(data) 

print("introduce the value of time for the highest value of glycaemia (tMg)") #tMg=time for the highest value of glycaemia
tMg<-as.numeric(scan(file="",what="",nmax=1))

print("introduce the value of the highest value of glycaemia (GMg)") #GMg= the highest value of glycaemia
GMg<-as.numeric(scan(file="",what="",nmax=1))

print("introduce the value of time for the lowest value of glycaemia (tmg)") #tmg=time for the lowest value of glycaemia
tmg<-as.numeric(scan(file="",what="",nmax=1))

print("introduce the value of the oral dose of glucose (D0)")	#D0 is the oral dose of glucose
D0<-as.numeric(scan(file="",what="",nmax=1))

print("introduce the value of time of the highest value of insulinemia (tMi)") #tMi=time of the highest value of insulinemia
tMi<-as.numeric(scan(file="",what="",nmax=1))

print("introduce the value of the highest value of insulinemia (IMi)") #IMi= the highest value of insulinemia
IMi<-as.numeric(scan(file="",what="",nmax=1))

Ga<-data[1,2] 

Ia<-data[1,3]

#ka and ke estimation
tabla1<-data.frame(subset(data,t>=tMg&t<=tmg,select=c(t,G)))
tabla1<-cbind(tabla1,ln=log(tabla1$G))
print(tabla1)
fit1<-lsfit(tabla1$t,tabla1$ln)
print(fit1$coefficients)
tabla2<-data.frame(subset(data,t<tMg,select=c(t,G)))
tabla2<-cbind(tabla2,lnGr=log(exp(fit1$coefficients[1]+fit1$coefficients[2]*tabla2$t)-tabla2$G))
print(tabla2)
fit2<-lsfit(tabla2$t,tabla2$lnGr)
ka<--(fit2$coefficients[2])
ke<--(fit1$coefficients[2])

#k0, Vd and k1 estimation
tabla3<-data.frame(subset(data,t<=15,select=c(t,G))) 
print(tabla3)
ajustek0<-nls(G~Ga+((k0*D0)/ka)*(1-exp(-ka*t)),data=tabla3,start=list(k0=0.01),trace=T)
k0<-coef(ajustek0)
tabla4<-data.frame(subset(data,t<=15,select=c(t,I)))
print(tabla4)
ajustek1<-nls(I~Ia+(((k0*D0)/ka)+Ga)*k1*t+((k0*D0*k1)/ka^2)*(exp(-ka*t-1)),data=tabla4,start=list(k1=0.01),trace=T)
k1<-coef(ajustek1)
Vd<-ka/k0

#k6 and Ipi estimation
tabla5<-data.frame(subset(data,t>=tMg,select=c(t,G))) 
print(tabla5)
ajustetpi<-nls(G~Ga+(GMg-Ga)/(1+exp((t-tpi)/B)),data=tabla5,start=list(tpi=120, B=10),trace=T)
tpi<-coef(ajustetpi)[1]
B<-coef(ajustetpi)[2]
print(tpi)
tabla6<-data.frame(subset(data,t>=tMi&t<=tmg,select=c(t,I)))
print(tabla6)
ajustek6<-nls(I~Ia+(IMi-Ia)*exp(-k6*(t-tMi)),data=tabla6,start=list(k6=0.0002),trace=T)
k6<-coef(ajustek6)
print(k6)
Ipi<-Ia+(IMi-Ia)*exp(-k6*(tpi-tMi))

#k2, k3, k4 estimation
dG1<-((Ga-GMg)*exp((-10)/B))/(B*(1+exp(-10/B))^2)
dG2<-((Ga-GMg)*exp((5)/B))/(B*(1+exp(5/B))^2)
I1<-Ia+(IMi-Ia)*exp(-k6*(tpi-10-tMi))
I2<-Ia+(IMi-Ia)*exp(-k6*(tpi+5-tMi))
k2<-(dG1-dG2)/(I2-I1)
k3<--dG1-k2*I1
k4<-k3/(Ipi-Ia)
print("Values of parameters")
writeLines(paste("ka=",ka))
writeLines(paste("ke=",ke))
writeLines(paste("k0=",k0))
writeLines(paste("k1=",k1))
writeLines(paste("k2=",-k2))
writeLines(paste("k3=",k3))
writeLines(paste("k4=",k4))
writeLines(paste("k6=",k6))
writeLines(paste("Ipi=",Ipi))
writeLines(paste("Vd=",Vd))
writeLines(paste("Ga=",Ga))
writeLines(paste("Ia=",Ia))
writeLines(paste("D0=",D0))