lumpedAgeModel.R

lumpedAgeModel.R

library(deSolve)



mod <- function(t,y,par){
  with(as.list(c(y,par)),{
    
    # lagged state variables
    if((t - TO) < 0){
      S_TO <- S_init
    } else {
      S_TO <- lagvalue((t - TO),2)
    }
    if((t - TO - TL) < 0){
      L_TL <- L_init
      S_TL <- S_init
    } else {
      L_TL <- lagvalue((t - TO - TL),1)
      S_TL <- lagvalue((t - TO - TL),2)
    }
    if((t - TO - TL - TP) < 0){
      S_TP <- S_init
    } else {
      S_TP <- lagvalue((t - TO - TL - TP),2)
    }
    if((t - TP) < 0){
      omega_TP <- omega_init
    } else {
      omega_TP <- lagvalue((t - TP),3)
    }
    if((t - TL) < 0){
      L_TL_o <- L_init
    } else {
      L_TL_o <- lagvalue((t - TL),1)
    }
     
    # dx/dt
    dL <- (lambda*S_TO*theta_O) - (lambda*S_TL*theta_O*theta_L*omega) - (mu_L*L) - ((gamma*L)*L)
    dS <- (lambda*S_TP*theta_O*theta_L*omega_TP*theta_P) - (mu_S*S)
    dOmega <- omega*((gamma*L_TL_o) - (gamma*L))
    
    # return dx/dt vector
    list(c(dL,dS,dOmega))
  })
}

# equilibrium

# fixed parameters
TO <- 1
TL <- 14
TP <- 1
mu_O <- 0.05
mu_L <- 0.1
mu_P <- 0.05
mu_S <- 0.1
lambda <- 30
gamma <- 0.001

theta_O <- exp(-TO*mu_O)
theta_L <- exp(-TL*mu_L)
theta_P <- exp(-TP*mu_P)

Lambda <- (lambda*theta_O*theta_L*theta_P)/mu_S

omega_init <- (1/TL)*log(Lambda)
L_init <- log(Lambda)/(gamma*TL)
S_init <- (L_init*(mu_L + (gamma*L_init))) / ((lambda*theta_O) - (mu_S/theta_P))

y_init <- c("L"=L_init,"S"=S_init,"omega"=omega_init)
par <- c("T0"=TO,"TL"=TL,"TP"=TP,"mu_O"=mu_O,"mu_L"=mu_L,"mu_P"=mu_P,"lambda"=lambda,"gamma"=gamma,
         "theta_O"=theta_O,"theta_L"=theta_L,"theta_P"=theta_P)

out <- dede(y = y_init,times = 0:800,func = mod,parms = par)

par(mfrow=c(2,1))
matplot(out[,2:3],type="l",lty=1,col=c("blue","red"))
plot(out[,"S"],type="l",col="red")
par(mfrow=c(1,1))