Unlimited growth scenarios
This section simulates several Lemna growth scenarios under unlimited growth conditions using the Klein et al. model. Model dynamics and numerical results are compared to results of the Schmitt et al. model.
Seven days exposure, seven days recovery
Figure 2 in Schmitt et al. (2013) presents several scenarios which simulate the growth of Lemna exposed to several concentrations of metsulfuron-methyl. Exposure lasts for seven days, followed by a seven day recovery period. Symbols show observed data and lines as calculated with fitted model parameters.
# simulate the 7d exposure, 7d recovery experiment
simulate_7d <- function(model=c("schmitt","klein"), ...) {
model <- match.arg(model)
time <- seq(0, 14, 0.2) # output time points
levels <- c(0, 0.32, 0.56, 1, 1.8, 3.2, 5.6) # simulated exposure levels
# set parameters as needed to replicate Figure 2 from Schmitt et al.
if(model == "schmitt") {
param$k_phot_fix <- TRUE # unlimited growth
param$k_phot_max <- 0.4
param$k_resp <- 0
} else if(model == "klein") {
param <- metsulfuron$param
param$k_photo_fixed <- TRUE
param$k_photo_max <- 0.4
param$k_loss <- 0
envir <- metsulfuron$envir
}
result <- data.frame()
for(conc in levels) {
exposure <- data.frame(time=c(0, 7, 7.01, 14), c=c(conc, conc, 0, 0))
if(model == "schmitt") {
param$Conc <- exposure
out <- calcgrowth(time, c(BM=0.0012, E=1, M_int=0), param=param, ...)
} else if(model == "klein") {
envir$conc <- exposure
out <- lemna(times=time, init=c(BM=0.0012, M_int=0), param=param, envir=envir, ...)
}
out$level = as.character(conc)
result <- bind_rows(result, out)
}
result
}
# simulate using the Schmitt et al. model
df_s <- simulate_7d(model="schmitt")
# recreate figure 2 from paper, including observed data
ggplot(df_s)+
geom_line(aes(time,FrondNo,color=level))+
geom_point(aes(time,FrondNo,color=level,shape=level), data=schmitt77)+
scale_y_log10(breaks=10^seq(-1,6))+
scale_x_continuous(breaks=seq(0,14,2))+
scale_shape_manual(values=c(0,4,2,8,1,3,5), name="Conc. (ug/L)")+
coord_cartesian(ylim=c(10,10000))+
geom_vline(xintercept=7, linetype="dotted", color="#888888", size=1)+
geom_segment(aes(x=0,y=10000,xend=6.9,yend=10000),arrow=arrow(length=unit(0.3,"cm"),type="closed"),linejoin="mitre",size=1,color="#888888")+
geom_segment(aes(x=7.1,y=10000,xend=14,yend=10000),arrow=arrow(length=unit(0.3,"cm"),type="closed"),linejoin="mitre",size=1,color="#888888")+
geom_text(aes(x=3.5,y=8000,label="Exposure",hjust="middle"),color="#888888")+
geom_text(aes(x=10.5,y=8000,label="Recovery",hjust="middle"),color="#888888")+
theme_bw()+
labs(x="Time (days)", y="Number of Fronds", color="Conc. (ug/L)",
title="Schmitt et al. model: 7d exposure, 7d recovery")
#> Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
#> ℹ Please use `linewidth` instead.
#> This warning is displayed once every 8 hours.
#> Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
#> generated.
The depicted model dynamics of the Schmitt et al. are identical to published results. Simulating the same scenarios using the Klein et al. model yields identical dynamics:
# simulate using the Klein et al. model
df_k <- simulate_7d(model="klein")
# recreate figure 2 from paper, including observed data
ggplot(df_k)+
geom_line(aes(time,FrondNo,color=level))+
geom_point(aes(time,FrondNo,color=level,shape=level), data=schmitt77)+
scale_y_log10(breaks=10^seq(-1,6))+
scale_x_continuous(breaks=seq(0,14,2))+
scale_shape_manual(values=c(0,4,2,8,1,3,5), name="Conc. (ug/L)")+
coord_cartesian(ylim=c(10,10000))+
geom_vline(xintercept=7, linetype="dotted", color="#888888", size=1)+
geom_segment(aes(x=0,y=10000,xend=6.9,yend=10000),arrow=arrow(length=unit(0.3,"cm"),type="closed"),linejoin="mitre",size=1,color="#888888")+
geom_segment(aes(x=7.1,y=10000,xend=14,yend=10000),arrow=arrow(length=unit(0.3,"cm"),type="closed"),linejoin="mitre",size=1,color="#888888")+
geom_text(aes(x=3.5,y=8000,label="Exposure",hjust="middle"),color="#888888")+
geom_text(aes(x=10.5,y=8000,label="Recovery",hjust="middle"),color="#888888")+
theme_bw()+
labs(x="Time (days)", y="Number of Fronds", color="Conc. (ug/L)",
title="Klein et al. model: 7d exposure, 7d recovery")
A comparison of numerical values between the two models shows that simulated biomass deviates at most by 0.15%. A plot of relative errors describes small but consistent deviations in biomass results:
df_s %>%
mutate(error = 100 * (df_k$BM/BM - 1)) %>%
ggplot(aes(time, error, color=level)) +
geom_point() +
theme_bw() +
labs(x="Time (days)", y="Relative error (%)", color="Conc. (ug/L)",
title="Deviation in biomass in Klein model relative to Schmitt")
The deviations are a result of imprecisions introduced by the
numerical ODE solver. Numerical precision can be increased by, for
example, decreasing the solver’s step length in time. This can be
achieved by decreasing the solver parameter hmax from the
model’s default value of hmax=0.1 to
e.g. hmax=0.01. The value of 0.01 is
equivalent to a step length of about 15 minutes. This decreases the
observed deviations significantly to a maximum of about 0.007%:
df_s <- simulate_7d(model="schmitt", hmax=0.01)
df_k <- simulate_7d(model="klein", hmax=0.01)
df_s %>%
mutate(error = 100 * (df_k$BM/BM - 1)) %>%
ggplot(aes(time, error, color=level)) +
geom_point() +
theme_bw() +
labs(x="Time (days)", y="Relative error (%)", color="Conc. (ug/L)",
title="Deviation in biomass in Klein model relative to Schmitt, hmax=0.01")
It can be concluded that model dynamics and numerical results are identical as far as numerical precision allows.
Two days of exposure, twelve days of recovery
The right-hand side of Figure 3 in Hommen et al. (2015) presents additional scenarios which simulate the growth of Lemna exposed to several concentrations of metsulfuron-methyl. Exposure lasts for two days, followed by a twelve day recovery period. Symbols show observed data and lines as calculated with fitted model parameters:
# simulate the 2d exposure, 12d recovery experiment
simulate_2d <- function(model=c("schmitt","klein"), ...) {
model <- match.arg(model)
time <- seq(0, 14, 0.2) # output time points
levels <- c(0, 0.35, 0.875, 1.75, 3.5) # simulated exposure levels
# set parameters as needed to replicate Figure 3 from Hommen et al.
if(model == "schmitt") {
param$k_phot_fix <- TRUE # unlimited growth
param$k_phot_max <- 0.295
param$k_resp <- 0.00
param$mass_per_frond <- 0.0004
} else if(model == "klein") {
param <- metsulfuron$param
param$k_photo_fixed <- TRUE
param$k_photo_max <- 0.295
param$k_loss <- 0
param$r_DW_FN <- 0.0004
envir <- metsulfuron$envir
}
result <- data.frame()
for(conc in levels) {
exposure <- data.frame(time=c(0, 2, 2.01, 14), c=c(conc, conc, 0, 0))
if(model == "schmitt") {
param$Conc <- exposure
out <- calcgrowth(time, c(BM=12 * param$mass_per_frond, E=1, M_int=0), param=param, ...)
} else if(model == "klein") {
envir$conc <- exposure
out <- lemna(times=time, init=c(BM=12 * param$r_DW_FN, M_int=0), param=param, envir=envir, ...)
}
out$level = as.character(conc)
result <- bind_rows(result, out)
}
result
}
# simulate using the Klein et al. model
df_k <- simulate_2d(model="klein")
# recreate figure 3 from paper, including observed data
ggplot(df_k)+
geom_line(aes(time,FrondNo,color=level))+
geom_point(aes(time,FrondNo,color=level,shape=level), data=hommen212)+
scale_y_log10(breaks=10^seq(-1,6))+
scale_x_continuous(breaks=seq(0,14,2))+
scale_shape_manual(values=c(0,4,2,8,1,3,5), name="Conc. (ug/L)")+
coord_cartesian(ylim=c(10,1000))+
geom_vline(xintercept=2, linetype="dotted", color="#888888", size=1)+
geom_segment(aes(x=0,y=1000,xend=1.9,yend=1000),arrow=arrow(length=unit(0.3,"cm"),type="closed"),linejoin="mitre",size=1,color="#888888")+
geom_segment(aes(x=2.1,y=1000,xend=14,yend=1000),arrow=arrow(length=unit(0.3,"cm"),type="closed"),linejoin="mitre",size=1,color="#888888")+
geom_text(aes(x=1,y=800,label="Exposure",hjust="middle"),color="#888888")+
geom_text(aes(x=8,y=800,label="Recovery",hjust="middle"),color="#888888")+
theme_bw()+
labs(x="Time (days)", y="Number of Fronds", color="Conc. (ug/L)",
title="Klein et al. model: 2d exposure, 12d recovery")
Observed model dynamics of the Klein et al. model are identical to published results and fit equally well with the observed frond numbers reported by Hommen et al. A comparison of numerical values of simulated biomass yields:
# simulate using the Schmitt et al. model
df_s <- simulate_2d(model="schmitt", hmax=0.01)
# simulate using the Klein et al. model
df_k <- simulate_2d(model="klein", hmax=0.01)
# calculate relative
df_s %>%
mutate(error = 100 * (df_k$BM/BM - 1)) %>%
pull(error) %>%
abs() %>%
max()
#> [1] 0.004476714Numerical biomass values deviate at most by about 0.004% between the two models. It can be concluded that model dynamics and numerical results of both models are identical for the presented scenarios.
Complex exposure patterns
Figure 4 in Schmitt et al. (2013) depicts the growth of Lemna exposed to three exposure patterns of several different concentrations. The patterns include constant exposure, repeating intervals of four days out of seven of exposure, as well as repeating intervals of two days out of seven of exposure.
First, the constant exposure pattern (Figure 4a) is simulated:
# simulate a specific exposure pattern
simulate_pattern <- function(model=c("schmitt","klein"), levels, expo_fun, ...) {
model <- match.arg(model)
time <- seq(0, 50, 0.2) # output time points
# set parameters as needed to replicate Figure 2 from Schmitt et al.
if(model == "schmitt") {
param$k_phot_fix <- TRUE # unlimited growth
param$k_phot_max <- 0.295
param$EC50 <- 0.39
param$k_resp <- 0
} else if(model == "klein") {
param <- metsulfuron$param
param$k_photo_fixed <- TRUE
param$k_photo_max <- 0.295
param$EC50_int <- 0.39
param$k_loss <- 0
envir <- metsulfuron$envir
}
result <- data.frame()
for(conc in levels) {
exposure <- expo_fun(conc)
if(model == "schmitt") {
param$Conc <- exposure
out <- calcgrowth(time, c(BM=0.0006, E=1, M_int=0), param=param, hmax=0.01)
out$Area <- out$BM * param$AperBM
} else if(model == "klein") {
envir$conc <- exposure
out <- lemna(times=time, init=c(BM=0.0006, M_int=0), param=param, envir=envir, hmax=0.01)
out$Area <- out$BM * param$r_A_DW
}
out$level = as.character(conc)
result <- bind_rows(result, out)
}
result
}
# constant exposure
expo_const <- function(conc=1) {
data.frame(time=c(0, 50), c=c(conc, conc))
}
# exposure concentrations used for constant patterns
levels_const <- c(0, 0.1, 0.25, 0.5, 1)
# simulate constant exposure patterns
df_p1 <- simulate_pattern(model="klein", levels=levels_const, expo_fun=expo_const)
ggplot(df_p1)+
geom_line(aes(time,Area,color=level))+
scale_y_log10(breaks=10^seq(-1,7),minor_breaks=NULL)+
scale_x_continuous(breaks=seq(0,50,10))+
scale_shape_manual(values=c(0,4,2,8,1,3,5))+
coord_cartesian(ylim=c(0.1,10^7))+
theme_bw()+
labs(x="Time (days)", y="Frond area (cm²)", color="Conc. (ug/L)",
title="Klein et al. model: constant exposure")
Next, the four out of seven days exposure pattern (Figure 4b) is simulated:
# exposure pattern: repeats 4d exposure, 3d recovery for 42 days,
# followed by recovery for 8 days
expo_4of7 <- function(conc=1) {
tim <- c(0, 4, 4.01, 6.99)
exp <- c(conc, conc, 0, 0)
tim_all <- c()
exp_all <- c()
for(i in seq(0,35, 7)) {
tim_all <- c(tim_all, tim+i)
exp_all <- c(exp_all, exp)
}
data.frame("time"=c(tim_all, 50), "conc"=c(exp_all, 0))
}
# exposure concentrations used for 4 out of 7d patterns
levels_4of7 <- c(0, 0.175, 0.438, 0.875, 1.75)
df_p2 <- simulate_pattern(model="klein", levels=levels_4of7, expo_fun=expo_4of7)
ggplot(df_p2)+
geom_line(aes(time,Area,color=level))+
scale_y_log10(breaks=10^seq(-1,7),minor_breaks=NULL)+
scale_x_continuous(breaks=seq(0,50,10))+
scale_shape_manual(values=c(0,4,2,8,1,3,5))+
coord_cartesian(ylim=c(0.1,10^7))+
theme_bw()+
labs(x="Time (days)", y="Frond area (cm²)", color="Conc. (ug/L)",
title="Klein et al. model: 4 out of 7d of exposure")
Finally, the two out of seven days exposure pattern (Figure 4c) is simulated:
# exposure pattern: repeats 2d exposure, 5d recovery for 42 days,
# followed by recovery for 8 days
expo_2of7 <- function(conc=1) {
tim <- c(0, 2, 2.01, 6.99)
exp <- c(conc, conc, 0, 0)
tim_all <- c()
exp_all <- c()
for(i in seq(0,35, 7)) {
tim_all <- c(tim_all, tim+i)
exp_all <- c(exp_all, exp)
}
data.frame("time"=c(tim_all, 50), "conc"=c(exp_all, 0))
}
# exposure concentrations used for 2 out of 7d patterns
levels_2of7 <- c(0, 0.35, 0.875, 1.75, 3.5)
df_p3 <- simulate_pattern(model="klein", levels=levels_2of7, expo_fun=expo_2of7)
ggplot(df_p3)+
geom_line(aes(time,Area,color=level))+
scale_y_log10(breaks=10^seq(-1,7),minor_breaks=NULL)+
scale_x_continuous(breaks=seq(0,50,10))+
scale_shape_manual(values=c(0,4,2,8,1,3,5))+
coord_cartesian(ylim=c(0.1,10^7))+
theme_bw()+
labs(x="Time (days)", y="Frond area (cm²)", color="Conc. (ug/L)",
title="Klein et al. model: 2 out of 7d of exposure")
Comparing numerical results of the Klein et al. model with the Schmitt et al. implementation yields the following relative errors in percent:
# relative error (%) for the constant exposure pattern
simulate_pattern(model="schmitt", levels=levels_const, expo_fun=expo_const) %>%
mutate(error = 100 * (df_p1$BM/BM - 1)) %>%
pull(error) %>%
abs() %>%
max()
#> [1] 0.0005302434
# relative error (%) for the 4 out of 7d exposure pattern
simulate_pattern(model="schmitt", levels=levels_4of7, expo_fun=expo_4of7) %>%
mutate(error = 100 * (df_p2$BM/BM - 1)) %>%
pull(error) %>%
abs() %>%
max()
#> [1] 0.03496461
# relative error (%) for the 2 out of 7d exposure pattern
simulate_pattern(model="schmitt", levels=levels_2of7, expo_fun=expo_2of7) %>%
mutate(error = 100 * (df_p3$BM/BM - 1)) %>%
pull(error) %>%
abs() %>%
max()
#> [1] 0.02197463It can be observed that the Klein et al. model reproduces the same model dynamics and biomass amounts as reported by Schmitt et al. under unlimited growth conditions.
