#Background
#Bormann e al 1977 investigated nitrogen (N) dynamics as one of the most important limiting nutrients for ecosystem production of the Northern hardwood forest ecosystem. Modeling N cycle further elucidates the biogeochemistry of N in a temperate forest ecosystem and the ecological response in aggrading vegetation considering the fluxes between different pools.
#Modeling structures
#Attempts were made to represent the model structure of N in the study area and the interactive effects among different pools. In particular, model equations of fluxes and parameterization of the Bormann model were done to study the effect of in- and -out fluxes on the patterns of vegetation, available and bound N pools over time.
#Modeling stages (Shown in R codes)Bormann original model#???   Constraint (vegetation carrying capacity correction) ???    Calibration (optimization and parameter sensitivity) ???    What happens if the system is disturbed?
#Base_model digram

getwd()
## [1] "C:/Users/Paul/Documents/Rwork"
##Load DiagrammeR
library(DiagrammeR)
# a digraph is a directed graph
# the model is written in the "dot" language
Bormann_model <-
  "digraph{
graph[overlap=true,fontsize=10]

Vegetation_532 -> AvailableN_26 [label='7.4'];

Vegetation_532-> Boundpool_4700[label='63.3'];

AvailableN_26-> Vegetation_532 [label='79.6'];
Boundpool_4700-> AvailableN_26 [label='69.6'];
Boundpool_4700-> External_output_0 [label='0.1'];
External_output_0->Boundpool_4700 [label='14.2'];
External_output_0 -> AvailableN_26 [label='6.5'];
AvailableN_26 -> External_output_0 [label='3.9'];
}"
grViz(Bormann_model)
#Base_model diagram in the second assignment link
##The response shows that vegetation growth is indefinite which doesn't represent the actual system. So, next is to add a constraint to check the time of aggrading of the vegetation pool.
#code with carrying capacity
library (deSolve)
borman_base <- function(t,y,p){with(as.list(c(y,p)),
                                    {dV.dt<-a12*A*V*(1-(V/K))-a21*V-a31*V
                                    dA.dt<-F20 + a21*V + a23*B - a12*A*V - a02*A
                                    dB.dt<-F30 + a31*V - a23*B - a03*B
                                    return(list(c(dV.dt,dA.dt,dB.dt)))})}
#Define p,y,t
#p
p <- c(F20=6.5,
       F30=14.2,
       a12=79.6/(26*532),
       a31=63.6/532,
       a21=7.4/532,
       a23=69.6/4700,
       a02=3.9/26,
       a03=0.1/4700,
       A=26,
        B=4700,
       K=1000)
y <- c(V=500,A=20,B=4000) #inital state of the system
#t <- c(0,1,10,100)#timepoints we want output from
t=seq(0,500,by=0.5)
out2 <- ode(y=y ,time=t, func=borman_base,parms=p)
plot(out2)

##From here we used bormann.logistic as an input function from file BormannLogistic.R

bormann.logistic <- function(t, y, p) {
  
  with(as.list( c(y, p) ), {
    ### ODEs
    ## We hypothesize that the growth function is all that Bormann et al. included. 
    ## Therefore, our logistic inhibition term acts on the entire forest.
    dV.dt <- (a12 * A * V  - a21 * V - a31 * V) * (1 - V/K )
    ## 
    
    ## A still loses N to vegetation, so we leave the loss function the same.
    dA.dt <- ( a20 + a23 * B ) - ( a12 * V * A + a02 * A)
    
    ## 
    dB.dt <- ( a30 + a31 * V ) - ( a23 * B + a03 * B )
    loss <- a02*A + a03*B
    # returning a list whose first (and only) element 
    # is a vector of the ODE values
    return( 
      list(
        c( dV.dt, dA.dt, dB.dt ), loss=loss
      )
    )
  }
  )
}
One percent deviations in some parameters cause substantial variation.

One percent deviations in some parameters cause substantial variation. 

##Next tested sensitivity of model
source("BormannLogistic.R") #Sourcing the file mentioned above
bormann.logistic #Pull the function from the above file into the workspace
## function(t, y, p) {
##   
##   with(as.list( c(y, p) ), {
##     ### ODEs
##     ## We hypothesize that the growth function is all that Bormann et al. included. 
##     ## Therefore, our logistic inhibition term acts on the entire forest.
##     dV.dt <- (a12 * A * V  - a21 * V - a31 * V) * (1 - V/K )
##     ## 
##     
##     ## A still loses N to vegetation, so we leave the loss function the same.
##     dA.dt <- ( a20 + a23 * B ) - ( a12 * V * A + a02 * A)
##     
##     ## 
##     dB.dt <- ( a30 + a31 * V ) - ( a23 * B + a03 * B )
##     loss <- a02*A + a03*B
##     # returning a list whose first (and only) element 
##     # is a vector of the ODE values
##     return( 
##       list(
##         c( dV.dt, dA.dt, dB.dt ), loss=loss
##       )
##     )
##   }
##   )
## }
#Running the model using following pyt
p <- c( a20 = 6.5, a30 = 14.2, a12 = 79.6/(532*26), a21 = (6.6 + 0.8)/532, 
        a31 = (54.2 + 2.7 + 6.2 + 0.1)/532,
        a23 = 69.6/4700, a02 = 3.9/26, a03 = 0.1/4700,
        K = 1000,
          A=26,
        B=4700)

y <- c(V = 532, A = 26, B = 4700)

t <- seq(from = 0, to = 500, by = 1)

out3 <- ode(y = y, times = t, func = bormann.logistic, parms = p)

plot( out3 ) #plot initial model
One percent deviations in some parameters cause substantial variation.

One percent deviations in some parameters cause substantial variation. 

##The following is to find a steady state to test the changes against
#Our model -- looks pretty steady-state-ish to me. We will drop the first 500 years to create our reference."}

y <- c(V = 532, A = 26, B = 4700)
t <- seq(from = 0, to = 1000, by = 1)
out4 <- ode(y = y, times = t, func = bormann.logistic, parms = p)
plot(out4)
One percent deviations in some parameters cause substantial variation.

One percent deviations in some parameters cause substantial variation. 

Reference <- as.data.frame(out4-(1:500))
## Warning in out4 - (1:500): longer object length is not a multiple of
## shorter object length
## Make small deviations to see the difference in parameters
#One percent deviations in some parameters cause substantial variation."}
## we make a clean COPY of the parameters and make sure p is not a list, 
## in case we created it as a list in the first place....
pp <- unlist( p ) 

## proportional amount by which we change parameters, i.e. 1%
tiny <- 0.01 

## Make a graphics device with 3 rows and 3 columns of figures
par( mfrow=c(3,3) )

## number of parameters
npar <- length(p)

## do the following for each parameter, i
for (i in 1:npar ) {
  
  ## make a slightly bigger parameter value
  p[i]  <- pp[i]*(1+tiny)
  
  ## run the model with a bigger p[i], dropping the first 500 obs.
  outp <- as.data.frame(
    ode(y = y, times=t, fun=bormann.logistic, parms = p) 
  )[ -(1:500), ]
  
  # run the model with a smaller p[i], dropping the first 500 obs.
  p[i]  <- pp[i]*(1-tiny)
  outm <- as.data.frame(
    ode(y = y, times=t, fun=bormann.logistic, parms = p) 
  )[ -(1:500), ]
  my.max <- 5900 
  my.min  <- 5700
  
  par(mar = c(2,2,2,2)) #changes margins of figures
  
  ## plot the result for parameter i
  plot( x = Reference$time, y = Reference$B, xlab="year", ylab="kg/m2", 
        type="l", main = names(p)[i], ylim = c(my.min, my.max))
  
  ## add a blue dashed line for plus tiny 
  lines(outp$time, outp$B, lty = 2, col = "blue")
  
  ## add a red dotted line for minus tiny
  lines(outm$time, outm$B, lty=3, col="red") 
  
  ## We RESET p[i] to its original value using our clean copy
  p[i] <- pp[i]
  
} # This is the loop to go through all parameters
One percent deviations in some parameters cause substantial variation.

One percent deviations in some parameters cause substantial variation. 

#Create reference to consider a numerical approach to sensitivity

out5 <- as.data.frame(
  ode(y = y, times=t, fun=bormann.logistic, parms = p) 
)
reference <- out5[-(1:500),]


pp <- unlist(p) 
state   <- y
yRef    <- reference$B
pp      <- unlist(p)
nout    <- length(yRef)
npar    <- length(p)

## now we use a truly small value: 0.000001%
tiny    <- 1e-8
dp      <- pp*tiny # the relative changes

## a matrix to hold our sensitivity values 
Sens    <- matrix(NA, nrow=nout, ncol=npar)
## name the columns by the parameters
colnames(Sens) <- names(p)

for (i in 1:npar)
{
  ## change the parrameter value
  dval    <- pp[i] + dp[i]
  p[i] <- dval
  
  ## run the model, and drop the first 500
  yPert   <- as.data.frame(ode(state, t, fun = bormann.logistic, p))$B[-(1:500)] 
  
  ##Calculates the sensitivities at all time points
  Sens[,i]<- (yPert-yRef)/dp[i]*pp[i] 
  
  ##...AND return p to its original values
  p[i] <- pp[i]
}

## Name the rows as the time points.
rownames(Sens) <- t[-(1:500)]

## show examples of sensitivities at arbitrary time points
format(as.data.frame(Sens[c(1, 61, 121, 181, 241),]), digits=2)
##      a20  a30 a12   a21  a31   a23  a02 a03    K A B
## 500 1342 3548 737 -2105 3331 -4403 -730 -34 1152 0 0
## 560 1402 3713 769 -2250 3634 -4773 -762 -36 1357 0 0
## 620 1443 3834 790 -2365 3932 -5127 -784 -38 1579 0 0
## 680 1466 3913 802 -2450 4224 -5464 -797 -40 1817 0 0
## 740 1473 3953 805 -2506 4507 -5785 -800 -41 2070 0 0
## make a picture
matplot(as.numeric(row.names(Sens)), Sens, type = 'l', 
        xlab = "Days", ylab = "Sensitivities")
legend("right", legend = colnames(Sens), lty =  1:5, col = 1:6  )

## Calculate quanititative summaries
mabs <- colMeans( abs(Sens) ) # or apply( Sens, 2, function(x) mean( abs( x ) ) )
msqr <- sqrt( colSums( Sens * Sens ) / nout ) 
format(data.frame(msqr,mabs),digits=2)
##     msqr mabs
## a20 1429 1428
## a30 3855 3854
## a12  781  781
## a21 2432 2429
## a31 4576 4528
## a23 5833 5785
## a02  776  776
## a03   40   40
## K   2243 2153
## A      0    0
## B      0    0
## Make a picture of the summed quantities
One percent deviations in some parameters cause substantial variation.

One percent deviations in some parameters cause substantial variation. 

##Optimized model for accurate representation of system
library(deSolve)
source("BormannLogistic.R")
bormann.logistic
## function(t, y, p) {
##   
##   with(as.list( c(y, p) ), {
##     ### ODEs
##     ## We hypothesize that the growth function is all that Bormann et al. included. 
##     ## Therefore, our logistic inhibition term acts on the entire forest.
##     dV.dt <- (a12 * A * V  - a21 * V - a31 * V) * (1 - V/K )
##     ## 
##     
##     ## A still loses N to vegetation, so we leave the loss function the same.
##     dA.dt <- ( a20 + a23 * B ) - ( a12 * V * A + a02 * A)
##     
##     ## 
##     dB.dt <- ( a30 + a31 * V ) - ( a23 * B + a03 * B )
##     loss <- a02*A + a03*B
##     # returning a list whose first (and only) element 
##     # is a vector of the ODE values
##     return( 
##       list(
##         c( dV.dt, dA.dt, dB.dt ), loss=loss
##       )
##     )
##   }
##   )
## }
p <- c( a20 = 6.5, a30 = 14.2, a12 = 79.6/(532*26), a21 = (6.6 + 0.8)/532,
        a31 = (54.2 + 2.7 + 6.2 + 0.1)/532,
        a23 = 69.6/4700, a02 = 3.9/26, a03 = 0.1/4700,
        K = 600)
y <- c(V = 532, A = 26, B = 4700)
t <- seq(from = 0, to = 500, by = 1)
out6 <- ode(y = y, times = t, func = bormann.logistic, parms = p)
plot( out6 )
One percent deviations in some parameters cause substantial variation.

One percent deviations in some parameters cause substantial variation. 

sse.bormann.mineralization <- function(params, data) {
  p['a23'] <- params
  ## run the model
  out7 <- ode(y = y, times = t, func = bormann.logistic, parms = p)
  ## store the last values of the state variables
  nr <- nrow(out7)
  model.output <- out7[nr,c("V", "A", "B")]
  ## Calculate the sum of the squared differences
  ## squaring the differences makes them positives and
  ## weights big differences even more heavily
  diffs <- model.output - data
  diffs2 <- diffs^2
  sse <- sum( diffs2 )
  ## Return the SSE
  sse
}
data = c(V = 532, A=26, B=4700)
fit0 <- optimize(f = sse.bormann.mineralization, interval = c(0.01, 0.5), data=data)
fit0
## $minimum
## [1] 0.018149
## 
## $objective
## [1] 4625.243
## $minimum
## [1] 0.018149
##
## $objective
## [1] 4625.243
p.opt <- p
p.opt['a23'] <- fit0$minimum
out.original <- as.data.frame( ode(y = y, times = 0:500, func = bormann.logistic, parms = p))
out.opt <- as.data.frame( ode(y = y, times = 0:500, func = bormann.logistic, parms = p.opt))
layout(matrix(1:4, nrow=2) )
plot(out.original$time, out.original$V, type="l")
lines(out.opt[,1], out.opt[,2], lty=2, lwd=1, col=2)
abline(h=data['V'], lty=3, lwd=1, col=4)
plot(out.original$time, out.original$A, type="l", ylim = c(23,33))
lines(out.opt[,1], out.opt[,3], lty=2, lwd=1, col=2)
abline(h=data['A'], lty=3, lwd=1, col=4)
plot(out.original$time, out.original$B, type="l", ylim=c(1000, 6000))
lines(out.opt[,1], out.opt[,4], lty=2, lwd=1, col=2)
abline(h=data['B'], lty=3, lwd=1, col=4)
sse.bormann.mineralization2 <- function(params, data) {
  ## an objective function whose arguments are
  ## params - the parameter of interest
  ## data - data from Bormannm et al.
  ## Assign a new value of the parameter to our set of parameters
  p['a23'] <- params[1]
  ## run the model
  out8 <- ode(y = y, times = t, func = bormann.logistic, parms = p)
  ## store the last values of the state variables
  nr <- nrow(out8)
  model.output <- out8[nr,c("V", "A", "B")]
  ## Calculate the sum of the squared differences
  ## squaring the differences makes them positives and
  ## weights big differences even more heavily
  diffs <- log(model.output, 10) - log(data, 10)
  diffs2 <- diffs^2
  sse <- sum( diffs2 )
  ## Return the SSE
  sse
}
fit0.log <- optimize(f = sse.bormann.mineralization2,
                     interval = c(0.001, 0.1), data=data)
fit0.log
## $minimum
## [1] 0.0181761
## 
## $objective
## [1] 0.002796935
## $minimum
## [1] 0.0181761
##
## $objective
## [1] 0.002796935
fit0
## $minimum
## [1] 0.018149
## 
## $objective
## [1] 4625.243
## $minimum
## [1] 0.018149
##
## $objective
## [1] 4625.243
p.opt.log <- p
p.opt.log['a23'] <- fit0.log$minimum
out.opt.log <- as.data.frame( ode(y = y, times = 0:500, func = bormann.logistic, parms = p.opt.log))
rbind(data=data, Bormann = out.opt[501, 2:4], raw=out.opt[501,2:4], log=out.opt.log[501,2:4])
##                V        A        B
## data    532.0000 26.00000 4700.000
## Bormann 599.9441 25.49684 4702.931
## raw     599.9441 25.49684 4702.931
## log     599.9447 25.49695 4695.946
## V A B
## data 532.0000 26.00000 4700.000
## Bormann 599.7254 25.46203 5753.202
## raw 599.9441 25.49684 4702.931
## log 599.9447 25.49695 4695.946

sse.bormann.mineralization3 <- function(params, data, vars) {
  
  p['a23'] <- params[1]
  
  ## run the model 
  out9 <- ode(y = y, times = t, func = bormann.logistic, parms = p)
  
  ## store the last values of the state variables
  nr <- nrow(out9)
  model.output <- out9[nr,c("V", "A", "B")]
  
  ## Calculate the sum of the squared differences
  ## squaring the differences makes them positives and 
  ## weights big differences even more heavily
  diffs <- log(model.output, 10) - log(data, 10)
  diffs2 <- diffs^2 
  wdiffs2 <- diffs2/vars
  sse <- sum( wdiffs2 )
  
  ## Return the SSE
  sse
}
vars <- c(532/10, 26/10, 4700*10)
fitv <- optimize(f = sse.bormann.mineralization3, interval = c(0.001, 0.1), data=data, vars=vars)
cbind(fit0, fit0.log, fitv)
##           fit0     fit0.log    fitv        
## minimum   0.018149 0.0181761   0.02453208  
## objective 4625.243 0.002796935 7.808268e-05
##           fit0     fit0.log    fitv        
## minimum   0.018149 0.0181761   0.02453208  
## objective 4625.243 0.002796935 7.808268e-05

######### Multiple Parameters ####
sse.b.m.p2 <- function(params, data) {
  ## an objective function whose arguments are 
  ## params - the parameter of interest
  ## data - data from Bormannm et al.
  
  ## Assign a new value of the parameter to our set of parameters
  p['a23'] <- params[1]
  p['a31'] <- params[2]
  ## run the model 
  outnp <- ode(y = y, times = t, func = bormann.logistic, parms = p)
  
  ## store the last values of the state variables
  nr <- nrow(outnp)
  model.output <- outnp[nr,c("V", "A", "B")]
  
  ## Calculate the sum of the squared differences
  ## squaring the differences makes them positives and 
  ## weights big differences even more heavily
  diffs <- log(model.output, 10) - log(data, 10)
  diffs2 <- diffs^2 
  sse <- sum( diffs2 )
  
  ## Return the SSE
  sse
}
params <- c(a32 = 69.6/4700, a31 = 0.19)
fit2 <- optim(par = params, fn = sse.b.m.p2, data=data)
fit2
## $par
##        a32        a31 
## 0.01854462 0.12182748 
## 
## $value
## [1] 0.002724496
## 
## $counts
## function gradient 
##       65       NA 
## 
## $convergence
## [1] 0
## 
## $message
## NULL
## $par
##        a32        a31 
## 0.01854462 0.12182748 
## 
## $value
## [1] 0.002724496
## 
## $counts
## function gradient 
##       65       NA 
## 
## $convergence
## [1] 0
## 
## $message
## NULL

########Variance weighted
sse.bormann.mineralization3 <- function(params, data, vars) {
  ## an objective function whose arguments are 
  ## params - the parameter of interest
  ## data - data from Bormannm et al.
  
  ## Assign a new value of the parameter to our set of parameters
  p['a23'] <- params[1]
  
  ## run the model 
  outm <- ode(y = y, times = t, func = bormann.logistic, parms = p)
  
  ## store the last values of the state variables
  nr <- nrow(outm)
  model.output <- outm[nr,c("V", "A", "B")]
  
  ## Calculate the sum of the squared differences
  ## squaring the differences makes them positives and 
  ## weights big differences even more heavily
  diffs <- log(model.output, 10) - log(data, 10)
  diffs2 <- diffs^2 
  wdiffs2 <- diffs2/vars
  sse <- sum( wdiffs2 )
  
  ## Return the SSE
  sse
}
vars <- c(532, 26, 4700)
fitv <- optimize(f = sse.bormann.mineralization3, 
                 interval = c(0.001, 0.1), data=data, vars=vars)
cbind(fit0, fit0.log, fitv)
##           fit0     fit0.log    fitv        
## minimum   0.018149 0.0181761   0.01831603  
## objective 4625.243 0.002796935 7.888169e-06
##           fit0     fit0.log    fitv        
## minimum   0.018149 0.0181761   0.01831603  
## objective 4625.243 0.002796935 7.888169e-06
One percent deviations in some parameters cause substantial variation.

One percent deviations in some parameters cause substantial variation. 

## Finally added in fire as disturbance
source("BormannLogistic.R")
p <- c( a20 = 6.5, a30 = 14.2, a12 = 79.6/(532*26), a21 = (6.6 + 0.8)/532, 
        a31 =  0.12182485, 
        a23 = 0.01823085, # previously optimized
        a02 = 3.9/26, a03 = 0.1/4700,
        K = 600)

y <- c(V = 532, A = 26, B = 4700)
t <- seq(from = 0, to = 500, by = 1)
outf <- ode(y = y, times = t, func = bormann.logistic, parms = p)
plot( outf )
One percent deviations in some parameters cause substantial variation.

One percent deviations in some parameters cause substantial variation. 

## event with data frame, 1 way to add fire
rows <- 5
eventdat <- data.frame(var = rep("V", rows), 
                       time = seq(100, 500, length=rows),
                       value = rep(100, rows),
                       method = rep("rep", rows)
)
eventdat
##   var time value method
## 1   V  100   100    rep
## 2   V  200   100    rep
## 3   V  300   100    rep
## 4   V  400   100    rep
## 5   V  500   100    rep
##Integrate the above
t <- 0:max(eventdat$time)
initial = y
outt <- ode(initial, times=t, fun=bormann.logistic, parms=p, 
            events=list(data=eventdat) # EVENTS as a data frame
)
plot(outt)
One percent deviations in some parameters cause substantial variation.

One percent deviations in some parameters cause substantial variation. 

##different plot of the above
matplot(out1[,1], out1[,-1], type="l", xlab="Time", ylab="N", log="y")
legend('right', colnames(out1[,-1]), lty=1:5, col=1:5, bty="n") 

## Add fire as an event
event.func <- function (t,y,p){
  V.burn <- y[1]*3/5*0.5 
  V.air <- V.burn/3
  V.bound <- V.burn/3
  V.avail <- V.burn/3
  B.burn <- y[3]*1/3*0.5
  B.avail <- B.burn/3
  A.burn <- y[2]*0.5
  A.air <- A.burn/1
  y[1] <- y[1] - V.burn
  y[2] <- y[2] + V.avail + B.avail - A.air
  y[3] <- y[3] + V.bound - B.burn
  return(y)}
fire <- round( seq(from=200, to=1000, by=200),0)
t <- 0:1000
initial <- c(V=100, A=30, B=4700)
outevent <- ode(initial, times=t, fun=bormann.logistic, parms=p,events=list(func=event.func, time=fire))
plot(outevent)
One percent deviations in some parameters cause substantial variation.

One percent deviations in some parameters cause substantial variation.

One percent deviations in some parameters cause substantial variation.

One percent deviations in some parameters cause substantial variation. 

#Conclusion
#The key mechanisms of N dynamics in the study site revealed the different fluxes among the pools and also check for balance of N for a steady state ecosystem