search

kdorheim / MENDplus

MEND R package and materials for the mechanistic reaction modeling analysis.
https://kdorheim.github.io/MENDplus/
0 stars 1 forks source link

Validate MEND configuration #3

Open kdorheim opened 4 years ago

kdorheim commented 4 years ago

Run a series of tests to

kdorheim commented 4 years ago

Hi @jianqiuz just wanted to give you an update I think I might have MEND up and running but when I try to reproduce figure 5 from Wang et al. 2013 the values returned for the DOC and Q pools were much larger than I expected. Please see the attachment some initial results. I will go through and double check the code and parameters. But if you have any ideas about some tests to try and run please let me know!

Reproduce_MEND.pdf

kdorheim commented 4 years ago

Ah yes so in our default case there are some carbon pools that end up being negative 😞 which seems wrong. However it is unclear if this is a problem with the parameters, the initial values, a numerical artifact, or something else. @jianqiuz have you seen something like this behavior in MEND before?

image

jianqiuz commented 4 years ago

My first try LOL

first try
kdorheim commented 4 years ago

Hmmm thanks, that helps! Glad to know with other configurations of MEND and these parameter values the M pool ends up being negative. What does the IC pool correspond to?

jianqiuz commented 4 years ago

I think my P and M are not coded correctly. IC is the inorganic carbon ( CO2) I wasn't able to run your code to get time series. Also curious what ODE solver you're using. Another note on the model, in Wang's paper, this setup we're using is for model spin up. the simulations plotted in Figure are based on at least different initial conditions

kdorheim commented 4 years ago

Awesome I will include IC in this as well. I am using the default ode from deSovle which allows me to use a Runge-Kutta method.

jianqiuz commented 4 years ago

rm(list = ls ())

library(deSolve)

solving MEND model with B, D, P, Q, M, EP, EM and IC (CO2)

Decom <- function(times, states, parameters, flux_function) { with(as.list(c(states, parameters)), { F1 <- (1/Ec) (Vd + mR) B D /( kD + D) ## microbial uptake F2 <- Vp EP P / (kP + P) # POM decomposition F3 <- Vm EM M / (kM + M) #MOAM decomposition F4 <- (1/Ec -1) Vd B D /( kD + D) #microbial growth respiration F5 <- (1/Ec -1) mR B D /( kD + D) ## microbial maintenance respiration F6 <- Kads D (1- Q/ Qmax) ##adsorption F7 <- Kdes Q/Qmax ##desorption F8 <- (1- pEP- pEM) mR B #microbial biomass decay F9ep <- pEP mR B #enzyme production F9em <- pEM mR B F10ep <- rEP EP #enzyme decay F10em <- rEM EM

dP <- Ip + (1 - gD) * F8- F2
dM <- (1 - fD) * F2 - F3
dQ <- F6 - F7
dB <- F1- (F4 + F5) - F8 - (F9em + F9em)
dD <- Id + fD * F2 + gD * F8 + F3 + (F10em + F10ep)- F1 - (F6 - F7)
dEP <- F9ep - F10ep
dEM <- F9em - F10em
dIC <- F4 + F5   #CO2 fllux
dTot <- Ip + Id -(F4 + F5) 
return(list(c(dP, dM, dQ, dB, dD, dEP, dEM, dIC, dTot)))

}) }

times <- seq(0, 8760, 1) ##per hour 24hour by 365 days =8760

B<-2 D<-1 P<-10 Q<-0.1 M<-5 EP<-0.00001 EM<- 0.00001 IC<- 0 Tot<- 18.10002

states <- c(B = B, D = D, P = P, Q = Q, M = M, EP = EP, EM = EM, IC = IC, Tot= Tot) parameters <- c(Ec = 0.47, Vd = 5e-4, mR = 2.8e-4, kD = 0.26, Vp = 2.5, kP = 50, Vm = 0.1, kM = 250, Kads = 0.006, Kdes = 0.001, Qmax = 1.7,pEP=0.01, pEM=0.01, rEP=1e-3, rEM=1e-3, Ip = 8e-5, Id = 8e-5, fD=0.5, gD=0.5)

DecomOut = ode(y = states, parms = parameters, times = times, func = Decom, method = "radau") print("done") plot(DecomOut)

jianqiuz commented 4 years ago

Here's my quick MEND code I wrote last night, as I was trying to copy your runs. Strangely the P and M simulations are off. Maybe you can help me trouble-shooting

kdorheim commented 4 years ago

Hmmm how does this look? This wasn't really what I was expecting but there aren't any negative pools any more which I guess is good. But the shape of the lines seem unusual.

image