############################################################################ #Example Rumen Passage and SI Digestion Model for NANP Modeling Workshop #This model contains CP, NDF, and Microbial CP pools #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) #Gut model - Rum plus SI Gut <- function(t, parms) { #Derivative section derivs <- function(t, state, parms) { with(as.list(c(state, parms)), { #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 #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),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)) }) } state <- c(QCpRum=0.120*6.25, QChoRum=4.5, QCpMiRum=0.186*6.25) #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.04; #per hour KCpAbs = 0.75; KChoAbs = 0.45 KRdp <- 0.5 pars <- c(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])