############################################################################
#Example Rumen Passage and SI Digestion Model for NANP Modeling Workshop
#This model contains a single CP pool. 
#Written by Mark D. Hanigan, Virginia Tech
#############################################################################

#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, 100, by = 4)

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

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

		#Rumen Fluxes
		Fin <- DMI * Cnut
		FRumDeg <- Qrum * Kdeg
		FRumSI <- Qrum * Kpas

		#Rumen Differentials
		dQrum <- Fin - FRumDeg - FRumSI

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

		#Digestibility Calculations
		RDnut <- FRumDeg/Fin			#Ruminal
		TDnut <- (Fin - FSIFec) / Fin		#Total Tract

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

#Model Inputs
DMI = 22/24;		#kg DM/h
Cnut = 0.196*0.16*.53	#kg/kg DM
Kdeg = 0.075;		#per hour
Kpas = 0.04;		#per hour
KSIAbs = 0.75;


pars <- c(DMI,Cnut,Kdeg,Kpas,KSIAbs); #Define the list of model inputs
out <- Gut(time, pars);		#Call the model passing time and pars vectors
 					#Model output is collected in the out matrix
summary(out);			#Summarize output
plot(out)

plot(out[,1:2])