# Day 4 V1
# Implementation of the DEBtox model to simulate survival, growth and reproduction under chemical exposure
# in this version we assume variable food and competition, but constant temperature, and continuous reproduction
#simulation settings--------------------------------------------------------------------------------
tmax <- 21 # simulation time [days]
res <- matrix(data = NA, nrow = tmax, ncol = 6) # results matrix
#total food availability mg/d
food.daily <-0.5
pop.feed<- rep(c(NA),tmax)
# Volume of environment in mL
env.size <- 1000
# exposure concentration
Cw <-rep(c(0), tmax )
#Cw [1:3]<-0.4
#Cw [7:9]<-0.4
pMoA <- 1 #Enter physiological mode of action considered for the chemCal effect, here are the ones that are implemented:
# 1: feeding/assimilation reduced
# 2: Maintanance costs increased
# 3: Cost for structure and maturation increased
# 4: Hazard during oogenesis/embryogenesis
pop.size <-1 #initial population size
# DEB parameters-----------------------------------------------------------------------------------
v <- 2.2 # energy conductance
g <- 0.92 # energy investment ratio
km <- 0.51 # somatic maintencance rate coefficalent
Lb <- 0.86 # length at birth
Lp <- 2.8 # length at puberty
R_m <- 37.8 # maximum reproduction rate
Lm <-v/(g*km)
s_shrink<- 0.8 # shrinking coefficient
Jxm <- 0.02503# scaled maximum ingestion rate
xk <- 0.0002# half saturation constant
#store parameter values at time 0 for calculating parameter changes due to stress
f<-1
km0<-km
R_m0<-R_m
f0<-1
g0<-g
#Stress function related parameters for chemical exposure---------------------------------------------------
par.kd <- 0.482 # elimination rate constant, d-1
par.z <- 0.112 # threshold (NEC) for survival (mg/L)
par.kk <- 1.21 # killing rate (L/mg/d) (SD only)
par.hb <- 0 #1.17E-02 # background hazard rate (1/d)
#par.kd <-0.482/(Lm-Lb) #elimination constant for length scaling
par.c0 <-0.0515 # no effect concentration (threshold for effect)
par.cT <-0.908 # tolerance concentration
food.total <- rep(c(food.daily),tmax)
#intialize/clear state variables vectors
ind.L <- c(NA) # initial volumetric structural length (here: start at birth)
ind.R <- c(NA) # initial cumulative reproduction (here: start at birt)
ind.Jx <- c(NA) # actual ingestion rate
ind.e <- c(NA) # scaled reserve density; e=f for constant food conditions
ind.Csi <- c(NA) # scaled internal concentration
ind.H <- c(NA) # Cumulative hazard at t=0
ind.S <- c(NA) # survival probability
ind.b_L <- c(NA) # volumetric structural length of embryo
ind.alive <- c(NA) # individual still alive?
ind.maxL <- c(NA) # maximum length over time
ind.ran <- c(NA) # random number for survival
ind.F_s <- c(NA) # stress fucntion for sublethal effects
#intial state variables of individuals in vectors
for (i in 1:pop.size ) {
ind.L[i] <- c(Lb) # initial volumetric structural length (here: start at birth)
ind.R[i] <- c(0) # initial cumulative reproduction (here: start at birt)
ind.Jx[i] <- c(0) # actual ingestion rate
ind.e[i] <- c(1) # scaled reserve density; e=f for constant food conditions
ind.Csi[i] <- c(0) # scaled internal concentration
ind.H[i] <- c(0) # Cumulative hazard at t=0
ind.S[i] <- c(1) # survival probability
ind.b_L[i] <- c(0) # volumetric structural length of embryo
ind.alive[i] <- c(TRUE) # individual still alive?
ind.maxL[i] <- c(0) # maximum length over time
ind.ran [i] <- c(sample(0:1000, 1) / 1000) # random number for survival
ind.F_s[i] <- c(0) # stress function for sublethal effects
}
#time loop---------------------------------------------------------------------------------------------
for (t in 1:tmax ) {
#number of indivduals in vector dead or alive
pop.size <- length(ind.L)
#---individual toxicokinetics and stress function---------------------------------------
for(i in 1:pop.size) {
#dose metric
ind.Csi[i]<-ind.Csi[i]+ par.kd*(Cw[t]-ind.Csi[i]) #scaled internal concentration
#ind.Csi[i]<-ind.Csi[i] + Lm/ind.L[i] *par.kd*(Cw[t]-ind.Csi[i]) #scaled internal concentration, length scaling
#The resulting stress function for a given concentration:
ind.F_s[i] <- 1/par.cT *max(0,ind.Csi[i]-par.c0)
#Parameter changes due to stress for pMoAs of fedding
if (pMoA == 1) {#feeding or assimilation from food
f <- f0 * max(0, 1-ind.F_s[i] ) }
}
#---population feeding--------------------------------------------
for(i in 1:pop.size) {
ind.Jx[i] <- (f *Jxm * ind.L[i]^ 2 * (food.total[t] / env.size)) / (xk + (food.total[t] / env.size))
}
# sum of feeding rates of all living individuals
pop.feed[t] <- sum(ind.Jx[which(ind.alive == TRUE)])
if(pop.feed[t] > food.total[t]){
for(i in 1:pop.size) {
ind.Jx[i] <- (food.total[t] / pop.feed[t]) * ind.Jx[i]
}
pop.feed[t] <- sum(ind.Jx[which(ind.alive == TRUE)])
}
#add leftovers to food of next day
if (t+1< tmax){
food.total[t+1] <- food.total[t+1] + food.total[t] - pop.feed[t] }
#---individual life histories--------------------------------------------
for(i in 1:pop.size) {
if (pMoA == 2) {#somatic and maturity maintenance
km <- km0 * (1+ind.F_s[i] )
R_m <- R_m0 / (1+ind.F_s[i] )
}
if (pMoA == 3) {#cost for structure and maturation
g <- g0 * (1+ind.F_s[i] )
km <- km0 /(1+ind.F_s[i] )
}
if (pMoA == 4) {#embryo hazard
R_m <- R_m0 * exp(-ind.F_s[i] )
}
#update scaled reserve density in iterations to minimize error
ind.f<-ind.Jx[i]/(Jxm*ind.L[i]^ 2)
for (x in 1:100 ) {ind.e[i]<-ind.e[i]+(ind.f-ind.e[i])*v/ind.L[i]/100 }
# Growth
dL<-km*g/(3*(ind.e[i]+g)) *(ind.e[i]*v/(km*g)-ind.L[i])
#maximum length over time
if (ind.maxL[i] < ind.L[i]) {ind.maxL[i]<-ind.L[i] }
#Reproduction
if (ind.L[i]>Lp) { dR<-(R_m/((v/(km*g))^3-Lp^3)) *((v/km *ind.L[i]^2+ind.L[i]^3 )*ind.e[i]/(ind.e[i]+g)-Lp^3)}
else dR<-0;
#Cumulative hazard
ind.H[i]<-ind.H[i]+par.kk*(max(0,ind.Csi[i]-par.z))+ par.hb + (s_shrink * km * max(0, ind.maxL[i] - ind.L[i])/ind.maxL[i])
#Survival probability
ind.S[i]<- exp(-ind.H[i])
#update growth and reproduction
ind.L[i]<-ind.L[i]+dL
if (dR>0) {ind.R[i]<-ind.R[i]+dR}
#individual still alive?
#if (ind.S[i]< ind.ran [i]) {ind.alive[i]<-FALSE }
#---calculate bood development-----------------
}#end individual life histories
#Store results for L and R
res[t,1]<-t
res[t,2]<-ind.S[i]
res[t,3]<-ind.L[i]
res[t,4]<-ind.R[i]
res[t,5]<-Cw[t]
res[t,6]<-ind.Csi[i]
}
#create plots for growth and reproduction---------------------------------------------------------
par(mfrow = c(3, 2))
barplot(res[,5], main="Exposure", ylab = "Concentration [mg/L]", xlab = "Time [d]")
#points(res[,6], lwd=3)
#plot scaled internal concentration
plot(res [,1], res [,6], ylab = "Concentration [mg/L]", xlab = "Time [d]", main = "Scaled internal concentration")
#plot survival curve curve
plot(res [,1], res [,2], type = "l",lwd=3, ylab = "Survival probability [-]", xlab = "Time [d]", main = "Survival", ylim=c(0, 1)) #
#points(dataS[,2],dataS[,3],pch=19,cex = 1.3)
#abline(v=8, col="red")
#plot growth curve
plot(res [,1], res [,3], type = "l",lwd=3, ylab = "Length [cm]", xlab = "Time [h]", main = "Growth",
ylim=c(0, 6) )
abline(v=8, col="red")
#lines(res [,1], Cw)
#plot reproduction
plot(res [,1], res [,4], type = "l",lwd=3, ylab = "Cumulative reproduction [#]", xlab = "Time [h]", main = "Reproduction",
ylim=c(0, 240) )
abline(v=8, col="red")
plot(res [,1], pop.feed, type = "l",lwd=3, ylab = "Feeding rate [mgC]", xlab = "Time [h]", main = "Feeding",
ylim=c(0, 1) )
