############################################################################
#Example Rumen Passage and SI Digestion Model for NANP Modeling Workshop
#This model contains CP, NDF, and Microbial CP pools and isotopic
#tracer pools for CP. 
#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, 300, by = 5)

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

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

		#Isotope Infusion
		if(t > TStartIsot & t < TStopIsot) {
			FCpIn_1 <- FCp_1
		} else {
			FCpIn_1 <- 0.0
		}

		#Rumen Fluxes
		FCpIn <- fDMI * cCp
		FCpDeg <- QCpRum * KdCp
		FCpSi <- QCpRum * Kp

		FChoIn <- fDMI * cCho
		FChoDeg <- QChoRum * KdCho
		FChoSi <- QChoRum * Kp

		FCpMi <- 0.25 * FChoDeg / (1 + KRdp*FChoDeg/FCpDeg)
		FCpMiSi <- QCpMiRum * Kp

		#Rumen Differentials
		dQCpRum <- FCpIn - FCpDeg - FCpSi
		dQChoRum <- FChoIn - FChoDeg - FChoSi
		dQCpMiRum <- FCpMi - FCpMiSi

		#Isootope enrichment, fluxes, and differentials
		#Always multiply the cold fluxes by enrichment to calculate the label fluxes
		eCpRum_1 <- QCpRum_1 / QCpRum;	#fractional enrichment
		eCpMiRum_1 <- QCpMiRum_1 / QCpMiRum

		dQCpRum_1 <- FCpIn_1 - FCpDeg*eCpRum_1 - FCpSi*eCpRum_1
		dQCpMiRum_1 <- FCpMi * eCpRum_1 - FCpMiSi * eCpMiRum_1

		#IntestInal Fluxes
		FAaAbs <- (FCpSi + FCpMiSi) * KCpAbs
		FCpFec <- FCpSi + FCpMiSi - FAaAbs 

		FGlcAbs <- FChoSi * KChoAbs
		FChoFec <- FChoSi - FGlcAbs 

		#Digestibility Calculations
		RDCp <- FCpDeg/FCpIn			#Ruminal
		RDCho <- FChoDeg/FChoIn

		TDCp <- (FCpIn - FCpFec) / FCpIn		#Total Tract
		TDCho <- (FChoIn - FChoFec) / FChoIn

		return(list(c(dQCpRum,dQChoRum,dQCpMiRum,dQCpRum_1,dQCpMiRum_1),
			FCpIn=FCpIn,FCpDeg=FCpDeg,FCpSi=FCpSi,FAaAbs=FAaAbs,RDCp=RDCp,
			TDCp=TDCp,FChoIn=FChoIn,FChoDeg=FChoDeg,FChoSi=FChoSi,
			FGlcAbs=FGlcAbs,RDCho=RDCho,TDCho=TDCho,FCpMi=FCpMi,
			FCpMiSi=FCpMiSi,eCpRum_1=eCpRum_1,eCpMiRum_1=eCpMiRum_1))
	   })
	}
	state <- c(QCpRum=0.120*6.25, QChoRum=4.5, QCpMiRum=0.186*6.25,QCpRum_1=0.0,
		QCpMiRum_1=0.0) #Observed rumInal particulate pool sizes 
	return(ode(y=state, times=t, func=derivs, parms=pars, method="rk4"))
}	#End of the Gut model function

#Model Inputs
fDMI = 22/24;		#kg DM/h
cCp = 0.196*(0.53+.114)	#kg/kg DM
cCho = 0.50
KdCp = 0.075;		#per hour
KdCho = 0.055
Kp = 0.05;			#per hour
KCpAbs = 0.75;
KChoAbs = 0.45
KRdp <- 0.1

TStartIsot <- 100
TStopIsot <- 196
FCp_1 <- 0.002;		# 2 g of isotope/hour

pars <- c(TStartIsot,TStopIsot,FCp_1,fDMI,cCp,KdCp,Kp,cCho, KdCho, KRdp, 
	KCpAbs, KChoAbs); #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])

Time <- out[,1]
Kp5 <- out[,21]

Kp = 0.08;			#per hour
pars <- c(TStartIsot,TStopIsot,FCp_1,fDMI,cCp,KdCp,Kp,cCho, KdCho, KRdp, 
	KCpAbs, KChoAbs); #Define the list of model inputs
out <- Gut(time, pars);		#Call the model passing time and pars vectors

Kp8 <- out[,21]

df <- data.frame(Time, Kp5, Kp8)
ggplot(df, aes(y=Kp5,x=Time)) + geom_point(linetype="solid")