## ====================================================================
## 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)