WETMASS change in DEB model of ectotherm

[drinky fan]

Hi everyone,

I am simulating body mass change by running DEB model in ectotherm(). Based on observed prey counts, I give which hour have food. And I simulated under two microclimate conditions, which have the same parameters but one warms about two celsius degree. I have two question for the simulation results. 

(1) At the start of simulation, there is no food and WETMASS is decreasing. The decrease rate is relatively high at beginning, but slow down at a time point (shown as the dashed line in the line plot). I noticed the air temperature also greatly drop at that inflection point. Can I attribute the slower WETMASS change to temperature-induced metabolism reduction? 

(2)With my data I, I found that the WETMASS and WETGUT at the end of simulate period is higher but WETGONAD, structure and reserve are smaller in warming environmental conditions. This phenomenon of WETMASS increase while structure reduction seems difficult to explain. What do you think about it?

I have put data and scripts that can reproduce above result in the attachment.

Thank you for help in advance.

Cheers,

Junqi

[drinky fan]

I can not upload the reproductive data in the topic because it exceeds the size limit. If anyone is interested in it I can email you.

[NicheMapR]

Hi Junqi,

It's hard to advise you without more details, e.g. what species you are simulating and at what life stage you are starting it at. It looks like it might be that you are starting at the egg stage which, if you just look at the wet mass  output will decline because reserve is being converted into structure (with a conversion efficiency below 1 as the 2nd law of thermodynamics necessitates). Then the mass jumps up because the animal starts feeding. I very strongly advise that you don't work with the DEB model without a full understanding of it - which can take some time and effort to acquire.

All the best,

Mike

[drinky fan]

Hi Mike,Thank you for the advice. 

I realize I did not provide enough context in previous, which led to a misunderstanding.To clarify: I am not simulating the egg stage. I am simulating an adult viper that faces a long period of fasting.The food is available only  around day 50.The decline in WETMASS I observed was the expected result of costing of reserve during the fast. I have two specific questions regarding the dynamics of this mass loss.

1. The change in rate of WETMASS loss (dW/dt)I observed two distinct phases during first fasting period (before it got food and immediately rise WETMASS). dW/dt is steady at first, then slows down significantly after a specific point (the red dashed line in my plot). I thought this might be temperature-driven,  as all energy flux follow arrhenius relationship. However, in the period before the red line, dW/dt remains constant even though the temperature fluctuates. If temperature is not the main driver here, is this sudden change in slope due to a state switch? 

2. Another question also about the calculation of WETMASS. I understand that the gut is part of Structure. The help document said "WETMASS - Wet mass total (reserve, structure, reproduction buffer, stomach contents) (g)". Can I understand that WETMASS is the sum of dV*V, WETGUT and w_E/mu_E * (E+ER) ?

Your insight would be very helpful for me.

Here is my code. The microclimate input is too large to upload in the attachment... I also can not using micro_global() to simulate a similar conditions. I put the local air tempeature in the plot previously. And I can send the microclimate file via email if you want to see it.

```

library(NicheMapR)

# parameters of ectotherm model ===
weight <- 140.6
pct_wet <- 0.02    
pct_cond <- 10
alpha_max <- 0.8
alpha_min <- 0.8 #
shape <- 1      
shape_b <- 10
T_F_min =  13.94
T_F_max = 28.29
T_B_min = 12
T_RB_min = 12
T_pref <- 25
CT_max <- 35.7
CT_min <- 1    
mindepth <- 2    
maxdepth <- 10    
shade_seek <- 1  
burrow <- 1    
climb <- 1    
minshade <- 0  
nocturn <- 1    
crepus <- 1      
diurn <- 1      
postur <- 0
M1 <- 0.000878; M2 <- 1.501; M3 <- 0.0165

# parameters of deb module ===
load('allstat.Rda')
species <- 'Gloydius.blomhoffii'
allDEB.species<-unlist(labels(allStat$allStat)) # get all the species names
allDEB.species<-allDEB.species[1:(length(allDEB.species)-2)] # last two elements are not species
species.slot <- which(allDEB.species == species)
par.names <- unlist(labels(allStat$allStat[[species.slot]]))
for(i in 1:length(par.names)){
  assign(par.names[i], unlist(allStat$allStat[[species.slot]][i]))
}
amp.param <- data.frame(param=NA, value=NA)
for(i in 1:length(par.names)) {
  amp.param[i,] <- c(par.names[i], unlist(allStat$allStat[[species.slot]][i]))
}
if(exists("E.Hj")==FALSE){E.Hj <- E.Hb}
if(exists("E.He")==FALSE){E.He <- E.Hb}
if(exists("L.j")==FALSE){L.j <- L.b}
kap.X <- 0.85
p.Xm <- p.Am / kap.X
z.mult <- 1
del_M <- 0.069345
T.L <- 1 + 273.15
T.H <- 40 + 273.15
T.AL <- 5E04
T.AH <- 9E04
ArrFunc5 <- function(x, T_A, T_AL, T_AH, T_L, T_H, T_REF, kdot_ref){
  #exp(T_A * (1 / T_REF - 1 / x)) / (1 + exp(T_AL * (1 / x - 1 / T_L)) + exp(T_AH * (1 / T_H - 1 / x))) * kdot_ref
 exp(T_A / T_REF - T_A / x) * ((1 + exp(T_AL / T_REF - T_AL / T_L) + exp(T_AH / T_H - T_AH / T_REF)) / (1 + exp(T_AL / x - T_AL / T_L) + exp(T_AH / T_H - T_AH / x))) * kdot_ref
}
Tbs <- seq(0, 50) + 273.15
TempCorr <- ArrFunc5(Tbs, T.A, T.AL, T.AH, T.L, T.H, T.ref, 1)
n.NC <- 1
n.NH <- 4/5
n.NO <- 3/5
n.NN <- 4/5
del.M <- 0.21957
V_init <- L.p^3
E_init <- E.m
E_H_init <- E.Hp+2
stage <- 2
viviparous <- 1
clutchsize <- 5
photostart <- 3
photofinish <- 2
E_sm <- 2400 / 0.3
foodmode <- 1
gutfill <- 20
J.prey <-  V_init * E_sm / 4
N.prey  <-  1

# read microclimate model ===
# we simulate the microclimate from September to following April
micro.recent <- readRDS('microclimate/2000.rds')
micro.warm <- readRDS('microclimate/2000_warm1.rds')

# hourly food ===
# based on observed count of prey of that year, we define which day the species can have prey (food_daily=1)
food.daily <- data.frame(
  yday = c(299:308, 310, 311),
  food_daily = 1
)

# convert to hourly energy that can be taken in
X <- rep(0, length=micro.recent$ndays)
X[food.daily$yday - micro.recent$soil[1,"DOY"]] <- 1
X <- X * J.prey*N.prey
X<-rep(X, each=24)

micro <- micro.recent
ecto <- ectotherm(DEB=1, starvemode = 1, intmethod=1, arrhen_mode = 1, aestivate = 1, depress = .2, aestdepth = 7, viviparous=viviparous, clutchsize=clutchsize, z.mult=z.mult, shape=shape, alpha_max=alpha_max, alpha_min=alpha_min, T_F_min=T_F_min, T_F_max=T_F_max, T_B_min=T_B_min, T_RB_min=T_RB_min, T_pref=T_pref, CT_max=CT_max, CT_min=CT_min, diurn=diurn, nocturn=nocturn, crepus=crepus, shade_seek=shade_seek, burrow=burrow, climb=climb, postur=postur, mindepth=mindepth, maxdepth=maxdepth, pct_wet=pct_wet, pct_cond=pct_cond, z=z*z.mult, E_sm = E_sm, foodmode=foodmode, gutfill=gutfill, X=X, del_M=del_M, p_Xm=p.Xm, kap_X=kap.X, v=v/24, kap=kap, p_M=p.M/24, E_G=E.G, kap_R=kap.R, k_J=k.J/24, E_Hb=E.Hb*z.mult^3, E_Hj=E.Hj*z.mult^3, E_Hp=E.Hp*z.mult^3, E_He=E.He*z.mult^3, h_a=h.a/(24^2), s_G=s.G, T_REF=T.ref, T_A=T.A, T_AL=T.AL, T_AH=T.AH, T_L=T.L, T_H=T.H, E_0=E.0*z.mult^4, f=f, d_V=d.V, d_E=d.E, d_Egg=d.E, mu_X=mu.X, mu_E=mu.E, mu_V=mu.V, mu_P=mu.P, kap_X_P=kap.P, n_X=c(n.CX,n.HX,n.OX,n.NX), n_E=c(n.CE,n.HE,n.OE,n.NE), n_V=c(n.CV,n.HV,n.OV,n.NV), n_P=c(n.CP,n.HP,n.OP,n.NP), n_M_nitro=c(n.CN,n.HN,n.ON,n.NN), L_b=L.b, V_init=V_init, E_init=E_init, E_H_init=E_H_init, stage=stage, K_egg = 0, pct_cond_egg = 100, spec_hyd_egg = 0, photostart = photostart, photofinish = photofinish, M_1=M1, M_2=M2, M_3=M3, write_input = 0)
micro <- micro.warm
ecto.warm <- ectotherm(DEB=1, starvemode = 1, intmethod=1, arrhen_mode = 1, aestivate = 1, depress = .2, aestdepth = 7, viviparous=viviparous, clutchsize=clutchsize, z.mult=z.mult, shape=shape, alpha_max=alpha_max, alpha_min=alpha_min, T_F_min=T_F_min, T_F_max=T_F_max, T_B_min=T_B_min, T_RB_min=T_RB_min, T_pref=T_pref, CT_max=CT_max, CT_min=CT_min, diurn=diurn, nocturn=nocturn, crepus=crepus, shade_seek=shade_seek, burrow=burrow, climb=climb, postur=postur, mindepth=mindepth, maxdepth=maxdepth, pct_wet=pct_wet, pct_cond=pct_cond, z=z*z.mult, E_sm = E_sm, foodmode=foodmode, gutfill=gutfill, X=X, del_M=del_M, p_Xm=p.Xm, kap_X=kap.X, v=v/24, kap=kap, p_M=p.M/24, E_G=E.G, kap_R=kap.R, k_J=k.J/24, E_Hb=E.Hb*z.mult^3, E_Hj=E.Hj*z.mult^3, E_Hp=E.Hp*z.mult^3, E_He=E.He*z.mult^3, h_a=h.a/(24^2), s_G=s.G, T_REF=T.ref, T_A=T.A, T_AL=T.AL, T_AH=T.AH, T_L=T.L, T_H=T.H, E_0=E.0*z.mult^4, f=f, d_V=d.V, d_E=d.E, d_Egg=d.E, mu_X=mu.X, mu_E=mu.E, mu_V=mu.V, mu_P=mu.P, kap_X_P=kap.P, n_X=c(n.CX,n.HX,n.OX,n.NX), n_E=c(n.CE,n.HE,n.OE,n.NE), n_V=c(n.CV,n.HV,n.OV,n.NV), n_P=c(n.CP,n.HP,n.OP,n.NP), n_M_nitro=c(n.CN,n.HN,n.ON,n.NN), L_b=L.b, V_init=V_init, E_init=E_init, E_H_init=E_H_init, stage=stage, K_egg = 0, pct_cond_egg = 100, spec_hyd_egg = 0, photostart = photostart, photofinish = photofinish, M_1=M1, M_2=M2, M_3=M3, write_input = 0)


# Question 1: why does the wetmass loss rate slow down around  Day 30?
# is this because the sudden drop in temperature led to a decrease in metabolism？
plot(x = ecto$debout[, "DAY"],
     y = ecto$debout[, "WETMASS"],
     type = "l",
     xlab = "Day",  
     ylab = "Wet Mass (g)",
     cex.lab = 2,        
     cex.axis = 1.5        
)
abline(v = 32, lty = 2, col = "red", lwd = 2)

# body temperature
plot(x = ecto$environ[, "DAY"],
     y = ecto$environ[, "TC"],
     type = "l",
     xlab = "Day",  
     ylab = "Body Temperature",
     cex.lab = 2,        
     cex.axis = 1.5        
)
abline(v = 32, lty = 2, col = "red", lwd = 2)

# environment temperature
metout <- micro.recent$metout
metout <- cbind(metout, DAY=rep(1:242, each=24))
plot(x = metout[,"DAY"],
     y = metout[, "TALOC"],
     type = 'l')
abline(v = 32, lty = 2, col = "red", lwd = 2)

mean.ta <- aggregate(
  x = metout[, "TALOC"],      
  by = list(DAY = metout[, "DAY"]),
  FUN = mean                  
)
plot(x =mean.ta$DAY,
    y = mean.ta$x,
    type = 'l',
    xlab = "Day",  
    ylab = "Daily averaged TALLOC",
    cex.lab = 2,        
    cex.axis = 1.5)
abline(v = 32, lty = 2, col = "red", lwd = 2)

```

Best regards,

Junqi

[NicheMapR]

Hi Junqi,

That extra detail helps and you are doing some clever things with the model. You are right about how the total wet mass is arrived at - structure converted from volume to wet mass, reserve converted from an energy density (J/cm^3) to mass, energy in the stomach and the reproduction buffer converted to wet mass, all summed. I suggest you try plotting the structure, reserve, reproduction buffer and gut contents to help understand what's happening (like the plot in the Shiny app) and additionally it might be worth you setting up a simpler simulation under constant environmental conditions just using the 'DEB_var' function. You should be able to work out what's going on from that. 

The distinction between structure and reserve is subtle in that it is a sub-cellular concept (and an abstraction), so it's not quite right to say that the gut is part of the structure. That said, during starvation the structure can eventually start being consumed for maintenance and this is what a lot of reptiles do in practice by reducing their gut and thus saving maintenance costs. Perhaps you'll see that the transition in rate of decline in mass is because a certain threshold has been reached where structure is being consumed.

Note that there are some aspects of the DEB model as implemented in the ectotherm function that are still a work in progress since the last release of the package, which I've been doing for some snake work of my own to allow food to be provided as periodic meals. It's not ready yet though. And I've been working on different starvation modules where reproduction buffer either not used at all (reproduction is the top priority), or only to avoid catabolising structure (reproduction is intermediate in priority), or always used whenever reserve dips down (long-term survival a priority). Again, not fully tested so not necessarily working properly.

All the best,

Mike

[drinky]

Hi Make, thanks for your advice.

After plotting the wetmass, reserve, structure and reproduction buffer together I can see that the growth flux and reproduction flux decrease (as |dV/dt| and |dE_R/dt| decrease) when reserve reached a threshold (in the red dashed line). Does this starve mode use structure and reproduction buffer simultaneously when reserve drop to a point? 

And I try to manually calculate the WETMASS using reserve, structure, reproduction buffer and stomach content. 

```

d.V < 0.3; w.E <- 23.9; mu.E <- 550000 # coming from AMP

R <- ecto$debout[,'E'] # ecto is the return of ectotherm(DEB=1,.....)
V<- ecto$debout[, 'V']
WETGUT <- ecto$debout[, 'WETGUT']
EH <- ecto$debout[,'E_H'] # maturity state
ER <- ecto$debout[, 'E_R']

w.E/mu.E * ((R+ER+EB)*structure) +V*d.V + WETGUT

```

But the sum is much lower than WETMASS. Sum of w.E/mu.E * ((ER+EB)*structure) is also lower than WETGONAD. Have I made mistake at any step?

Thank you for you help in advance.

Cheers,

Junqi

NicheMapR <nich...@googlegroups.com> 于2025年12月22日周一 12:28写道：

--
You received this message because you are subscribed to the Google Groups "NicheMapR" group.
To unsubscribe from this group and stop receiving emails from it, send an email to nichemapr+...@googlegroups.com.
To view this discussion, visit https://groups.google.com/d/msgid/nichemapr/121c6e8c-8dff-4395-b1be-13b18d0bd3b5n%40googlegroups.com.

[NicheMapR]

Hi Junqi,

You are getting close. d_V is the dry mass fraction of structure but that is for going from structure as energy to structure as wet mass via moles. When structure is considered in the dimension of volume it is assumed to have a density of water, so ~1 g / cm^3.

So you have 

d_E <- 0.3 # density of reserve, dimensionless

w_E <- 23.9 # molar mass of reserve, g/mol

mu_E <- 550000 # chemical potential of reserve, J/mol

E <- ecto$debout[,'E'] # reserve density, J/cm^3 (reserve per volume of structure)
V<- ecto$debout[, 'V'] # structural volume, cm^3
EB <- ecto$debout[,'E_B'] # reproduction batch, J
ER <- ecto$debout[, 'E_R'] # reproduction buffer, J

WETGUT <- ecto$debout[, 'WETGUT'] # wet mass of food in gut, g

To get the wet mass of reserve, where I'll use the DEB notation of reserve density with the square brackets:

E = [E] * V # J, convert reserve from energy density to energy

M_E = E / mu_E # moles, energy of reserve to moles of reserve

w_E_dry = M_E * w_E # g dry, dry mass of reserve

w_E_wet = w_E_dry / d_E # g wet, wet mass of reserve

I'll leave you to continue with the remaining calculations. Note that animals may have a lower water content of their reproduction buffer which would affect your overall calculation of wet mass. Reptiles often lay eggs with less water than the hatchlings need, and allow the water to be taken up from the soil by the eggs, thus saving space in their bodies so they can fit more eggs in. But in your case it's a viviparous snake I think so that won't be happening.

Of course with these calculations you need to pay close attention to the dimensions and units and it's easy to make a mistake. In the forthcoming Julia version of NicheMapR, everything must be assigned to dimensions and units, the code will error if you make a dimensional error and it will automatically convert between units.

All the best,
Mike

[drinky]

Hi Mike,

I think I'm almost there, but I am still seeing a small discrepancy between my manual calculations and the model’s direct output. My manually calculated initial wet mass is approximately 1g higher than the value returned by the model. I calculated it as "(w.E/mu.E * ((E_init + E_R_init + E_B_init)*V_init))/d.E + V_init*1", which is approximately 1g larger than the model's returned value. Do you think this difference come from estimation for structure density, or is there another component I might be missing in the manual summation?

And noticed that the returned WETGONAD variable is non-zero, even though I set the initial states E_R_init and E_B_init to 0. Since the help documentation defines WETGONAD as the mass of the reproduction and batch buffers, I am unsure why this value is not zero at the start.

Hope you have a wonderful Christmas with your families.

Best regards,

Junqi

NicheMapR <nich...@googlegroups.com> 于2025年12月23日周二 06:06写道：

You received this message because you are subscribed to a topic in the Google Groups "NicheMapR" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/nichemapr/SF3GC1K-Le8/unsubscribe.
To unsubscribe from this group and all its topics, send an email to nichemapr+...@googlegroups.com.
To view this discussion, visit https://groups.google.com/d/msgid/nichemapr/af215a61-a521-4be5-9f20-9a26465796f3n%40googlegroups.com.

[NicheMapR]

Hi Junqi,

That's great. Re the non-zero reproduction buffers - there could be a bug so I'll take a look. Re the 1 g discrepancy, have a look at the way it's calculated in the DEB.R code or the ectotherm Fortran code and you may find the reason.

All the best,
Mike

[drinky]

Hi Mike. Happy new year! Just let you know I have figured out figured out small discrepancy between my manually calculated wet mass (from reserve, structure, and reproduction buffer) and the output. DEB module in NicheMapR convert input molar energy of reserve to stoichiometric number and then to molar weight, whereas I mistakenly used w.E provided by AmP before. Thank you for you patience and time.

Best regards,

Junqi

NicheMapR <nich...@googlegroups.com> 于2025年12月28日周日 11:15写道：

To view this discussion, visit https://groups.google.com/d/msgid/nichemapr/c9e79082-4a85-4438-9294-af9c96da2f25n%40googlegroups.com.