############################################################################ #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 ############################################################################# #Abbreviations: F = flux, Q = pool size, dQ = differential of Q with respect to time, e = enrichment, # 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 # A _1 at the end of a variable denotes a labelled nutrient where the numeral can be used to # denote the position of the label, # #Time is in hours, mass in kg, and fluxes in kg/h require(deSolve); #Load the deSolve package) require(ggplot2) #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 #Isotope 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]) #Collect CP enrichment and time for plotting at 5 and 8% passage rates Time <- out[,1] eCP <- out[,21] df5 <- data.frame(Time, eCP) df5$Rate <- 5 #Run the model again at a passage rate of 8%/h 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 #Collect CP enrichment eCP <- out[,21] df8 <- data.frame(Time, eCP) df8$Rate <- 8 #Merge the data df <- rbind(df5, df8) df$Rate <- as.factor(df$Rate) #Convert rates to factors so they plot as a discrete category #Plot both enrichment curves on a common plot ggplot(df, aes(x=Time, y = eCP, colour = Rate)) + geom_line() #Read some observed data and plot with the model predictions #Set the working directory to the place on your computer containing the example data #Note that the normal \ denoting subdirectories in Windows must be reversed to / for R. setwd("L:/My Documents/") Obs2 <- read.csv("Exercise #1 data file.csv") #Plot enrichment curves on a common plot p <- ggplot() p <- p + geom_line(data=df, aes(x=Time, y = eCP, colour = Rate)) + geom_point(data=Obs2,aes(x=Time, y=Obs_eCP, color = "Obs")) p