#Example problem for NANP Modeling Workshop
#Time is in hours, mass in kg, and fluxes in kg/h

library(deSolve)			#Load the deSolve package)

#A list of times to output predictions 
time <- seq(0, 96, by = 4)

#Gut model - Rum plus SI
Gut <- function(t, parms) { 

	#Derivative section
	derivs <- function(t, state, parms) {
	   with(as.list(c(state, parms)), {

		#Constants and Fluxes
		Fin <- DMI * Cnut
		FRumDeg <- Qrum * Kdeg
		FRumSI <- Qrum * Kpas
		RDnut <- FRumDeg/Fin

		#Differentials
		dQrum <- Fin - FRumDeg - FRumSI

		#Intestinal Components
		FSIAbs <- FRumSI * KSIAbs
		FSIFec <- FRumSI - FSIAbs 

		TDnut <- (Fin - FSIFec) / Fin


		return(list(dQrum,Fin=Fin,FRumDeg=FRumDeg,FRumSI=FRumSI,RDnut=RDnut,FSIAbs=FSIAbs,TDnut=TDnut))
	   })
	}
	state <- c(Qrum=5.0)
	return(ode(y=state, times=t, func=derivs, parms=pars, method="rk4"))
}	#End of the Gut model function


#Model Inputs
DMI0 <- 1.0
Cnut0 <- 0.3
Kdeg0 = 0.032
Kpas0 = 0.025
KSIAbs0 = 0.95

tempres <- NULL
results <- NULL
loopnum <- NULL
for(i in 1:100) {
	run <- i

	#Random components
	DMI <- DMI0 + rnorm(1, mean = 0, sd = DMI0*.15)
	Cnut <- Cnut0 + rnorm(1, mean = 0, sd = Cnut0*.15)
	Kdeg <- Kdeg0 + rnorm(1, mean = 0, sd = Kdeg0*.15)
	Kpas <- Kpas0 + rnorm(1, mean = 0, sd = Kpas0*.15)
	KSIAbs <- KSIAbs0 + rnorm(1, mean = 0, sd = KSIAbs0*.15)

	pars <- c(DMI, Cnut, Kdeg, Kpas, KSIAbs)
	out <- Gut(time, pars)
	loopnum <- rbind(loopnum, i)
	tempres <- rbind(tempres, out[nrow(out),])	
}
results <- cbind(loopnum, tempres)
summary(results)