## ====================================================================
## Mammary Model Neal and Thornley, 1983
## includes changes made by Baldwin, 1995
## includes fixed substrate blood concentrations (caa cgl cAc cFa)
## 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
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
VAcTmV = 0.007, #vmax acetate denovo milk fat synthesis relative to cac
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)
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
#EFFECTS OF SUBSTRATE AVAILABILITY - Lactose
M = ULm / PCLm
Kminh = (Mm - M)/(Mm - M + Kr) #milk retention effects
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)
AcTmV = (VAcTmV * Cs * Kminh) / (1 + KAcTmV / cAc + K1AcTm / cGl)
AcTmV1 = AcTmV * 0.041667
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, 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
#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[,13:15], type = "l", lwd = 3, col = c(13,14,15), ylim=c(0,5),
xlab="Time,d", ylab="Amount,kg")
legend("topright", c("dMlkLm", "dMlkTm", "dMlkPm"), fill=c(13,14,15), 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[,10:12], type = "l", lwd = 3, col = c(10,11,12), ylim=c(0,60),
xlab="Time,d", ylab="Amount,kg")
legend("topright", c("Mave", "M", "dmilk"), fill=c(10,11,12), cex = .8)
matplot(x=df$time, y=df[,17:18], type = "l", lwd = 3, col = c(17,18), ylim=c(0,5),
xlab="Time,d", ylab="Amount,kg")
legend("topright", c("pTm", "pPm"), fill=c(17,18), 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)