## ==================================================================== ## Mammary Model Neal and Thornley, 1983 ## includes changes made by Baldwin, 1995 ## substrate blood concentrations can change (cAc) ## includes milk lactose (Lm) ## milk yield est from Lm ## includes milk fat (Tm, FaTmV=1.8, AcTmV=16.8) ## includes milk protein (Pm, AaPmV=8.6) ## ==================================================================== #load packages library(deSolve) library(FME) library(ggplot2) library(gridExtra) # SECTION 1. Initial conditions and parameters times = seq(0, 305, by = 0.1) #Specify parameters parms = c(Kh = 0.0102, #hormone decay rate for LHor /d (KLhor) Cu = 1000, #number of udder cells (ucells) Vm = 1, #max enzyme synthetic capacity /cell/d (Vusyn) KH = 0.2, #enzymatic response to hormone kg hormone/m3 (Kusyn) Ks = 0.1, #enzyme degradation rate constant /d (Kudeg) KR = 3, #milk secretion rate mm constant:persistency kg milk (KmlkI) Ksm = 0.2, #enzyme degrade rate constant indiced by milk retention /d (KudegM) Mh = 27, #half response pt for degrade due to Umilk kg (KmDeg) Kr = 0.048, #milk averaging constant /d (TaveM) Km = 0.00506, #milk secretion constant kg/uenz/d (VFaTmV, VAcTmV, VAaPmV, VGmLmV) KM = 4.43, #milk removal constant kg=4.43 (Kmilk=2.91/umlkcr) Mm = 30, #milk capacity of animal kg (Mlkmax) theta = 10, #parameter of equation (theta) Rc = 40, #milk removal by suckling calf kg/d (umlkcr=1) PCLm = 0.048, #percent lactose in milk VGlLmV = 0.0039, #maximal velocity glucose to lactose synthesis in mammary KGlLmV = 0.003, #Km glucose to lactose synthesis in mammary cGl = 0.0036, #blood concentration of glucose = 63mEq/dl KAaLmV = 0.002, #Km amino acids to lactose sytnesis in mammary cAa = 0.004, #blood concentration of total amino acids 0.0025-0.005 KAaPmV = 0.0021, #Km amino acids to milk protein sythesis in mammary VAaPmV = 0.0025, #vmax amino acids tomilk protein synthesis in mammary Kmilk = 2.91, #udder milk retention feedback effect kg/kg cFa = 0.00025, #blood concentration of fatty acids mole/L #cAc = 0.0017, #blood concentration of acetate mole/L vAc = 392.4, #body volume of distribution for Ac pool EBW*0.654 L KAaGlV = 0.01, #Km amino acids to gluconeogenesis VAaGlV = 0.228, #vmax amino acids to gluconeogenesis KAcTsF = 0.0018, #Km acetate to denovo body fat synthesis 0.0018 K1AcTs = 0.001, #vmax acetate to denovo body fat synthesis KFaTmV = 0.0005, #Km fat incorporated into milk fat relative to cfa VFaTmV = 0.001, #vmax fat incorporated into milk fat K1FaTm = 0.0015, #Km fat incorporated into milk fat relative to cgl KAcTmV = 0.0018, #Km acetate denovo milk fat synthesis in mammary rel to cac 0.0018 VAcTmV = 0.005, #vmax acetate denovo milk fat synthesis relative to cac 0.007 K1AcTm = 0.001) #Km acetate denovo milk fat synthesis relative to cgl #Specify initial values for state Variables qinit = c(H = 1, #hormone effector of udder enzymes kg/m3 (Lhor) #M = 0, #quantity of milk in cow kg (Umilk) ULm = 0.1, #udder milk lactose kg UTm = 0.01, #udder milk fat kg UPm = 0.01, #udder milk protein kg TMlkLm = 0, #total milk lactose kg TMlkTm = 0, #tot milk fat kg TMlkPm = 0, #tot milk protein kg Cs = 520, #number of udder enzymes (uenz) Ac = 0.005, #initial acetate pool size (EBW*0.0.0118) mole VAcTs = 42, #moles Ac to body fat synthesis? Mave = 0) #time avg or Umilk kg (Umave) #tvmlk = 0, #total milk yield kg # SECTION 2. Derivatives eqns = function(t, q, parms) { with(as.list(c(parms, q)), { #Eq1:lactation hormones decrease as lactation progresses in kg/m3 dH = -Kh * H #Eq5:udder enzyme growth influenced by lactation hormone Ps = Vm * Cu * H / (KH + H) #Eq6:udder enzyme death influenced by milk retention Ls = Cs * (Ks + Ksm * ((Mave / Mh)^theta / (1.0 + (Mave / Mh)^theta))) #Eq7:number of udder enzymes dCs = Ps - Ls M = ULm / PCLm Kminh = (Mm - M)/(Mm - M + Kr) #milk retention effects #ACETATE SUBSTRATE CONCENTRATIONS IN BLOOD cAc = Ac / vAc Ahor = cGl/0.003 absAc = 65.9 #20-70 (65.9 moles) dVAcTs = 12*Ahor-0.35*VAcTs AcTsF = VAcTs / (1+KAcTsF/cAc +K1AcTs/(Ahor*cGl)) #4-20 (16 moles) AaGlV = VAaGlV*121/(1+KAaGlV/cAa) AcTmV = (VAcTmV * Cs * Kminh) / (1 + KAcTmV / cAc + K1AcTm / cGl) #2-28 (16.8 moles) AcTmV1 = AcTmV * 0.041667 AaAc = AaGlV*0.6 #3.2 - 4.7 (2.4 moles) AcCd = 35.5 #24-34 (35.5 moles) dAc = absAc+AaAc-AcCd-AcTsF-AcTmV #MILK LACTOSE GlLmV = VGlLmV * Cs * Kminh / (1 + KGlLmV / cGl + KAaLmV / cAa) dMlkLm = ULm * Kmilk dULm = GlLmV * 0.5 * 0.342 - dMlkLm dmilk = dMlkLm / PCLm TMlkLm = dMlkLm tvmlk = TMlkLm / PCLm # MILK FAT FaTmV = (VFaTmV * Cs * Kminh) / (1 + KFaTmV / cFa + K1FaTm / cGl) FaTmV1 = FaTmV * 0.333 dMlkTm = UTm * Kmilk dUTm = AcTmV1 + FaTmV1 * 0.806 - dMlkTm TMlkTm = dMlkTm pTm = 100*dMlkTm / dmilk # MILK PROTEIN AaPmV = (VAaPmV * Cs * Kminh) / (1 + KAaPmV / cAa) dMlkPm = UPm * Kmilk dUPm = AaPmV * 0.110 - dMlkPm TMlkPm = dMlkPm pPm = 100*dMlkPm / dmilk #eq8:rate of secretion of milk kg/d is proportional to number of udder enzymes Cs Pm = Km * Cs * Kminh #eq9:rate of removal of milk in kg/d #dmilk = (M/(KM + M))*Rc #Eq14:rate of change in secretion of milk in kg/d is proportional to uenz and inhibited by retained milk #as approaches max capacity assuming no milk pulsing function m(t) in paper #dM = Pm - dmilk #Eq16:length of time M is high influences rate of secretion of milk in kg over 1/kr d dMave = Kr * (M - Mave) #tvmlk = dmilk res = c(dH, dULm, dUTm, dUPm, dMlkLm, dMlkTm, dMlkPm, dCs, dAc, dVAcTs, dMave) list(res) }) } # SECTION 3. Results integrate numerically out1 = rk(qinit, times, eqns, parms, method = "rk4") df = as.data.frame(out1) #Add non-integrated variables to dataframe df$M = df$ULm/(parms["PCLm"]) df$dmilk = df$ULm*parms["Kmilk"]/(parms["PCLm"]) df$dMlkLm = df$ULm*parms["Kmilk"] df$dMlkTm = df$UTm*parms["Kmilk"] df$dMlkPm = df$UPm*parms["Kmilk"] df$tvmlk = df$TMlkLm/(parms["PCLm"]) df$pTm = df$dMlkTm*100/df$dmilk df$pPm = df$dMlkPm*100/df$dmilk df$cAc = df$Ac/parms["vAc"] df$AaAc = parms["VAaGlV"]*121/(1+parms["KAaGlV"]/parms["cAa"])*0.6 df$AcTsF = df$VAcTs / (1+parms["KAcTsF"]/df$cAc + parms["K1AcTs"]/((parms["cGl"]/0.003)*parms["cGl"])) df$AcTmV = (parms["VAcTmV"] * df$Cs * (parms["Mm"] - df$M)/(parms["Mm"] - df$M + parms["Kr"])) / (1 + parms["KAcTmV"] / df$cAc + parms["K1AcTm"] / parms["cGl"]) #see what you have in the first rows head(df, 70) #see what you have in the last rows tail(df, 70) #1 plots colnames(df) matplot(x=df$time, y=df[,15:17], type = "l", lwd = 3, col = c(15,16,17), ylim=c(0,5), xlab="Time,d", ylab="Amount,kg") legend("topright", c("dMlkLm", "dMlkTm", "dMlkPm"), fill=c(15,16,17), cex = .8) matplot(x=df$time, y=df[,2:5], type = "l", lwd = 3, col = c(2,3,4,5), ylim=c(0,1), xlab="Time,d", ylab="Amount,kg") legend("topright", c("H", "ULm", "UTm", "UPm"), fill=c(2,3,4,5), cex = .8) matplot(x=df$time, y=df[,12:14], type = "l", lwd = 3, col = c(12,13,14), ylim=c(0,60), xlab="Time,d", ylab="Amount,kg") legend("topright", c("Mave", "M", "dmilk"), fill=c(12,13,14), cex = .8) matplot(x=df$time, y=df[,19:20], type = "l", lwd = 3, col = c(19,20), ylim=c(0,5), xlab="Time,d", ylab="Amount,kg") legend("topright", c("pTm", "pPm"), fill=c(19,20), cex = .8) matplot(x=df$time, y=df[,9], type = "l", lwd = 3, col = c(9), ylim=c(0,8000), xlab="Time,d", ylab="Amount") legend("topright", c("Cs"), fill=c(9), cex = .8) matplot(x=df$time, y=df[,21], type = "l", lwd = 3, col = c(21), ylim=c(0,5), xlab="Time,d", ylab="Amount") legend("topright", c("cAc"), fill=c(21), cex = .8) matplot(x=df$time, y=df[,22:24], type = "l", lwd = 3, col = c(22,23,24), ylim=c(0,50), xlab="Time,d", ylab="Amount") legend("topright", c("AaAc", "AcTsF", "AcTmV"), fill=c(22, 23, 24), cex = .8)