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

#Abbreviations: F = flux, Q = pool size, dQ = differential of Q with respect to time, 
#	C = concentration, K = mass action rate constant, Rum = rumen, SI = small intestine, 
#	Fec = fecal output, RD = ruminally degraded, TD = total tract digested, Abs = absorbed, 
#	In = input, Fa,b = the flux from pool a to pool b
#
#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)

#Digestion model function - Rum plus SI
#The model is defined as a function that can be called with input arguments of:
#	t, and parms.  Can be called by hand as below or by an optimizer solving for
#	parameter estimates or attempting to minimize or maximize any objective function
 
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, in this case 19.6% CP at 16% N and 53% ruminally insoluble
Kdeg = 0.075;		#degradation rater per hour
Kpas = 0.04;		#passage rate per hour
KSIAbs = 0.75;		#SI digestibilty of escaped nutrient


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

plot(out[,1:2])