## ====================================================================
## Mammary Model Neal and Thornley, 1983
## includes changes made by Baldwin, 1995
## includes substrate (caa cgl)
## includes milk lactose (Lm)
## milk yield est from Lm
##
##
## ====================================================================
#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 = 63 mg/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
Kmilk = 2.91) #udder milk retention feedback effect
#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
TMlkLm = 0, #total milk lactose 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
#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, dMlkLm, 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$tvmlk = df$TMlkLm/(parms["PCLm"])
#see what you have in the first rows
head(df, 70)
#see what you have in the last rows
tail(df, 70)
#1 Create plots
colnames(df)
matplot(x=df$time, y=df[,6:8], type = "l", lwd = 3, col = c(6,7,8), ylim=c(0,60),
xlab="Time,d", ylab="Amount,kg")
legend("topright", c("Mave", "M", "dmilk"), fill=c(6,7,8), cex = .8)
matplot(x=df$time, y=df[,5], type = "l", lwd = 3, col = c(5), ylim=c(0,8000),
xlab="Time,d", ylab="Amount")
legend("topright", c("Cs"), fill=c(5), cex = .8)
matplot(x=df$time, y=df[,2:3], type = "l", lwd = 3, col = c(2,3), ylim=c(0,1),
xlab="Time,d", ylab="Amount,kg")
legend("topright", c("H", "ULm"), fill=c(2,3), cex = .8)