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