# Single population : Single species model (DETERMINISTIC)
# RCN Marine Disease Modeling and Transmission Workshop
# CCPO, ODU, Norfolk, VA, USA (11-15 May 2015)
# R version 
# Bhagat Lal Dutta, Ifremer, La Tremblade, FR
# 12 May 2015
# 

require(deSolve)  # R package

## Initial conditions
Tmax <- 100 # length of the simulation in days
dT <- 1 # time increment in days to record the dynamics (NOT intergration step!!!!)


N <- 100 # Initial population size
I0 <- 1  # initial infection of 1 individual

#Indices of the variables
iT <- 1; iS <- 2; iI <- 3; iD <- 4; iP <- 5;



## Parameters

depth <- 1 # depth of water column

IPinfect <-0.003
Iinfect <-0.001
Dinfect <-0.0008
Imort <-8.0e-2
Bmort <-0
DeadDecay <-1.5
Irelease <-.015
Drelease <-1.0
IPremove <-0.001 
# parameters as a list
parms=list(IPinfect <- IPinfect, Iinfect <- Iinfect, 
           Dinfect <- Dinfect, Imort <- Imort, 
           Bmort <- Bmort, DeadDecay <- DeadDecay, 
           Irelease <- Irelease, h <- depth,
           Drelease <- Drelease, IPremove <- IPremove);

times  <- seq(0, Tmax, by=dT); # time span
abalone.init=c((N-I0), I0, 0,0); # initial population structure in terms of indivuals in each compartment


#define RHS of the model as a function
RHSabalone1 <- function(t, inipop, parms) {
  
  S <- inipop[1]  # Susceptibles
  I <- inipop[2]   # Infected and infectious
  D <- inipop[3]    # Dead
  P <- inipop[4]     # Particles in water
  
  with(as.list(parms), {
    dS <- -IPinfect*h*P*S -Iinfect*I*S - Dinfect*D*S - Bmort*S; # RHS of dS/dt in the equation
    dI <- +IPinfect*h*P*S +Iinfect*I*S + Dinfect*D*S - Imort*I; # RHS of dI/dt  "
    dD <-  Imort*I - DeadDecay*D; # RHS of dD/dt "
    dP <- Irelease*I + Drelease*D - IPremove*P; # RHS of "
    res<-c(dS, dI, dD, dP); 
    list(res)
  })
}

# RUN THE SIMULATION AND STORE THE RESULT AS DATAFRAME TABLE 
OUT <- as.data.frame(lsoda(abalone.init, times, RHSabalone1, parms));

#Store
S <- OUT[,iS]; I <- OUT[,iI]; D <- OUT[,iD]; P <- OUT[,iP];
t <- OUT[,iT];



### PLOTTING
par(mfrow=c(2,2))
par(mar=c(4,5,4,3))
plot(t, S, type="l", main="Susceptibles", 
     xlab="time (days)", ylab= expression(paste("number/", m^2, sep="")), 
     cex.lab=1.25, cex.axis=1.25, col="dark green", lwd=2)
plot(t, I, type="l", main="Infected", 
     xlab="time (days)", ylab= expression(paste("number/", m^2, sep="")), 
     cex.lab=1.25, cex.axis=1.25, col="red", lwd=2)
plot(t, D, type="l", main="Infected Dead", 
     xlab="time (days)", ylab= expression(paste("number/", m^2, sep="")), 
     cex.lab=1.25, cex.axis=1.25, col="dark red", lwd=2)
plot(t, P, type="l", main="Infectious Particles", 
     xlab="time (days)", ylab= expression(paste("number/", m^3, sep="")), 
     cex.lab=1.25, cex.axis=1.25, col="blue", lwd=2)


# Infection Rates Plot

par(mfrow=c(1,1))
par(mar=c(4,5,4,3))
infP=IPinfect * P * S;
infI=Iinfect *  I * S;
infD=Dinfect  * D * S;
plot(t,infP, type="l", lty=1, main="Infection rate (num/day)", xlab="day", 
     ylab="number/day", cex.lab=1.25, cex.axis=1.25)
lines(t, infI, lty=2)
lines(t,infD, lty=3)
legend("topright", c("P", "I", "D"), lty=c(1,2,3))