## ====================================================================
## Mammary Model Neal and Thornley, 1983
## No substrate
## No milking pulse function m(t) is Rc
## Milk yield only
##
##
##
## ====================================================================
# 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 (q in paper)
Rc = 40) #milk removal by suckling calf kg/d (umlkcr=1)
#Specify initial values for State Variables
qinit = c(H = 1, #hormone effector of udder enzymes kg/m3 (Lhor)
Cs = 520, #number of udder enzymes (uenz)
tvmlk = 0, #total milk production kg
M = 0, #quantity of milk in cow kg (Umilk)
Mave = 0) #time avg or Umilk kg (Umave)
# 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
#eq8:rate of secretion of milk kg/d is proportional to number of udder enzymes Cs
Kminh = (Mm - M)/(Mm - M + KR) #milk retention effects
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, dCs, tvmlk, dM, 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$dmilk = (df$M/(parms["KM"]+df$M))*parms["Rc"]
#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[,5:7], type = "l", lwd = 3, col = c(5,6,7), ylim=c(0,40),
xlab="Time,d", ylab="Amount,kg")
legend("topright", c("M", "Mave", "dmilk"), fill=c(5,6,7), cex = .8)
matplot(x=df$time, y=df[,3:4], type = "l", lwd = 3, col = c(3,4), ylim=c(0,8000),
xlab="Time,d", ylab="Amount")
legend("topright", c("Cs","tvmlk"), fill=c(3,4), cex = .8)