The hardware and bandwidth for this mirror is donated by dogado GmbH, the Webhosting and Full Service-Cloud Provider. Check out our Wordpress Tutorial.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]dogado.de.

Validation against FB4-Shiny

Hans Ttito

2026-06-09

Overview

This article verifies that fb4package produces results numerically identical to the original FB4-Shiny application (Deslauriers et al. 2017, v1.1.7) across all five simulation modes and across the nutrient and contaminant sub-models.

The reference implementation consists of the exact arithmetic from server.R of FB4-Shiny, adapted to accept R objects instead of CSV uploads. Both implementations receive identical inputs: same species parameters (pulled from fish4_parameters), same temperature series, same diet, and same initial conditions.

1. Shared setup

Species parameters — Chinook salmon

data(fish4_parameters)
chinook_key <- grep("tshawytscha", names(fish4_parameters), value = TRUE)[1]
db    <- fish4_parameters[[chinook_key]]
stage <- if ("adult" %in% names(db$life_stages)) "adult" else names(db$life_stages)[1]
sp    <- db$life_stages[[stage]]

cat("Species:", chinook_key, "| Stage:", stage, "\n")
#> Species: Oncorhynchus tshawytscha | Stage: adult
cat("CEQ:", sp$consumption$CEQ,
    " CA:", sp$consumption$CA,
    " CB:", sp$consumption$CB, "\n")
#> CEQ: 3  CA: 0.303  CB: -0.275

Scalar extraction helper

gp <- function(cat, nm) {
  v <- sp[[cat]][[nm]]
  if (is.null(v)) NA_real_ else as.numeric(v)
}

CEQ <- gp("consumption","CEQ"); CA  <- gp("consumption","CA")
CB  <- gp("consumption","CB");  CQ  <- gp("consumption","CQ")
CTO <- gp("consumption","CTO"); CTM <- gp("consumption","CTM")
CTL <- gp("consumption","CTL"); CK1 <- gp("consumption","CK1")
CK4 <- gp("consumption","CK4")

REQ <- gp("respiration","REQ"); RA  <- gp("respiration","RA")
RB  <- gp("respiration","RB");  RQ  <- gp("respiration","RQ")
RTO <- gp("respiration","RTO"); RTM <- gp("respiration","RTM")
RTL <- gp("respiration","RTL"); RK1 <- gp("respiration","RK1")
RK4 <- gp("respiration","RK4"); RK5 <- gp("respiration","RK5")

ACT  <- gp("activity","ACT");  BACT <- gp("activity","BACT")
SDA  <- gp("sda","SDA")

EGEQ <- gp("egestion","EGEQ"); FA <- gp("egestion","FA")
FB_  <- gp("egestion","FB");   FG <- gp("egestion","FG")
EXEQ <- gp("excretion","EXEQ"); UA <- gp("excretion","UA")
UB   <- gp("excretion","UB");  UG <- gp("excretion","UG")

PREDEDEQ <- gp("predator","PREDEDEQ")

Oxycal <- 13560

Derived temperature coefficients

if (!is.na(CEQ) && CEQ == 2) {
  CY <- log(CQ)*(CTM-CTO+2); CZ <- log(CQ)*(CTM-CTO)
  CX <- (CZ^2*(1+(1+40/CY)^0.5)^2)/400
}
if (!is.na(CEQ) && CEQ == 3) {
  CG1 <- (1/(CTO-CQ))*log((0.98*(1-CK1))/(CK1*0.02))
  CG2 <- (1/(CTL-CTM))*log((0.98*(1-CK4))/(CK4*0.02))
}
if (!is.na(REQ) && REQ == 2) {
  RY <- log(RQ)*(RTM-RTO+2); RZ <- log(RQ)*(RTM-RTO)
  RX <- (RZ^2*(1+(1+40/RY)^0.5)^2)/400
}

Simulation inputs

set.seed(123)
N_DIAS  <- 180
Init_W  <- 500     # g
p_fixed <- 0.4
ED_ini  <- 4200; ED_fin <- 4600

days     <- 1:N_DIAS
base_T   <- 10 + 4*sin(2*pi*(days-30)/N_DIAS)
Temp_vec <- round(pmax(2, base_T + rnorm(N_DIAS, 0, 0.3)), 2)

prey_prop_mat <- matrix(c(rep(0.7, N_DIAS), rep(0.3, N_DIAS)),
                        nrow = N_DIAS,
                        dimnames = list(NULL, c("Alewife","Inverts")))
prey_ED_mat   <- matrix(c(rep(4900, N_DIAS), rep(2600, N_DIAS)),
                        nrow = N_DIAS,
                        dimnames = list(NULL, c("Alewife","Inverts")))
ind_prey_mat  <- matrix(0, nrow = N_DIAS, ncol = 2,
                        dimnames = list(NULL, c("Alewife","Inverts")))
Pred_E_vec    <- seq(ED_ini, ED_fin, length.out = N_DIAS + 1)

cat(sprintf("%d days | W0 = %.0f g | p = %.2f | T = %.1f–%.1f °C\n",
            N_DIAS, Init_W, p_fixed, min(Temp_vec), max(Temp_vec)))
#> 180 days | W0 = 500 g | p = 0.40 | T = 5.6–14.6 °C

2. FB4-Shiny reference implementation

The functions below reproduce the arithmetic of server.R verbatim.

# --- Temperature functions ---
Cf1T <- function(T) exp(CQ*T)
Cf2T <- function(T) {
  if (T < CTM) { V <- (CTM-T)/(CTM-CTO); ft <- V^CX*exp(CX*(1-V)) } else ft <- 0
  max(ft, 0)
}
Cf3T <- function(T) {
  L1 <- exp(CG1*(T-CQ)); KA <- (CK1*L1)/(1+CK1*(L1-1))
  L2 <- exp(CG2*(CTL-T)); KB <- (CK4*L2)/(1+CK4*(L2-1))
  KA*KB
}
Cf4T <- function(T) max(exp(CQ*T + CK1*T^2 + CK4*T^3), 0)

Rf1T    <- function(T) exp(RQ*T)
RACTf1T <- function(W,T) {
  VEL <- if (T <= RTL) ACT*W^RK4*exp(BACT*T) else RK1*W^RK4*exp(RK5*T)
  exp(RTO*VEL)
}
Rf2T <- function(T) {
  if (T < RTM) { V <- (RTM-T)/(RTM-RTO); ft <- V^RX*exp(RX*(1-V)) } else ft <- 1e-6
  max(ft, 1e-6)
}

# --- Bioenergetic functions ---
consumo_fn <- function(T, W, p) {
  Cmax <- CA*W^CB
  ft   <- switch(as.character(CEQ),
                 "1" = Cf1T(T), "2" = Cf2T(T),
                 "3" = Cf3T(T), "4" = Cf4T(T))
  Cmax*p*ft
}

respiracion_fn <- function(T, W) {
  Rmax <- RA*W^RB
  if (REQ == 1) { ft <- Rf1T(T); ACT_v <- RACTf1T(W,T) } else { ft <- Rf2T(T); ACT_v <- ACT }
  Rmax*ft*ACT_v
}

egestion_fn <- function(C, T, p, i) {
  switch(as.character(EGEQ),
    "1" = FA*C,
    "2" = FA*(T^FB_)*exp(FG*p)*C,
    "3" = { PE  <- FA*(T^FB_)*exp(FG*p)
            PFF <- sum(ind_prey_mat[i,]*prey_prop_mat[i,])
            PF  <- ((PE-0.1)/0.9)*(1-PFF)+PFF; PF*C },
    "4" = FA*(T^FB_)*C)
}

excrecion_fn <- function(C, Eg, T, p) {
  switch(as.character(EXEQ),
    "1" = UA*(C-Eg),
    "2" =, "3" = UA*(T^UB)*exp(UG*p)*(C-Eg),
    "4" = UA*(T^UB)*(C-Eg))
}

# --- Daily growth loop ---
grow_shiny_daily <- function(Temperature, W, mode = "p_value", mode_value = 0.4,
                             prey_pm, prey_em, Fin) {
  out <- vector("list", Fin)
  for (i in seq_len(Fin)) {
    Pred_E_i  <- Pred_E_vec[min(i,   length(Pred_E_vec))]
    Pred_E_ip1 <- Pred_E_vec[min(i+1, length(Pred_E_vec))]
    mean_prey <- sum(prey_pm[i,] * prey_em[i,])

    if (mode == "p_value") {
      p       <- mode_value
      Cons_gg <- consumo_fn(Temperature[i], W, p)
    } else if (mode == "Ration") {
      Cons_gg <- mode_value
      Cmax    <- consumo_fn(Temperature[i], W, 1.0)
      p       <- if (Cmax > 0) Cons_gg/Cmax else 0
    } else if (mode == "Ration_prey") {
      Cons_gg <- mode_value / W
      Cmax    <- consumo_fn(Temperature[i], W, 1.0)
      p       <- if (Cmax > 0) Cons_gg/Cmax else 0
    }

    Cons <- Cons_gg * mean_prey
    Eg   <- egestion_fn(Cons, Temperature[i], p, i)
    Ex   <- excrecion_fn(Cons, Eg, Temperature[i], p)
    S    <- SDA*(Cons-Eg)
    Res  <- respiracion_fn(Temperature[i], W)*Oxycal

    G      <- Cons-(Res+Eg+Ex+S)
    egain  <- G*W
    finalwt <- (egain + Pred_E_i*W) / Pred_E_ip1
    if (is.na(finalwt) || finalwt < 0) finalwt <- 0

    out[[i]] <- list(
      Day                  = i,
      Starting_Weight      = W,
      Weight               = finalwt,
      Weight_change        = finalwt - W,
      Temperature          = Temperature[i],
      Mean_Prey_Energy_J_g = mean_prey,
      Consumption_gg       = Cons_gg,
      Consumption_energy   = Cons,
      Respiration          = Res,
      Egestion             = Eg,
      Excretion            = Ex,
      SDA                  = S,
      Net_energy           = G,
      Energy_density       = Pred_E_ip1,
      Cons_Alewife_g       = Cons_gg*W*prey_pm[i,"Alewife"],
      Cons_Inverts_g       = Cons_gg*W*prey_pm[i,"Inverts"]
    )
    W <- finalwt
  }
  as.data.frame(do.call(rbind, lapply(out, as.data.frame)))
}

# --- Binary search (Weight or Consumption target) ---
fit_p_shiny <- function(IW, target, fit_type, Temperature,
                        prey_pm, prey_em, Fin,
                        W_tol = 0.001, max_iter = 100) {
  p_min <- 0; p_max <- 5; p <- 0.5

  sim_val <- function(p_try) {
    df <- grow_shiny_daily(Temperature, IW, "p_value", p_try, prey_pm, prey_em, Fin)
    if (fit_type == "Weight") tail(df$Weight, 1)
    else sum(df$Consumption_gg * df$Starting_Weight)
  }

  val <- sim_val(p)
  for (iter in seq_len(max_iter)) {
    if (abs(val - target) <= W_tol) break
    if (val > target) p_max <- p else p_min <- p
    p <- (p_min + p_max)/2
    val <- sim_val(p)
  }
  list(p = p, converged = (abs(val-target) <= W_tol), final_val = val)
}

3. fb4package object

sp_pkg <- sp
sp_pkg$predator$PREDEDEQ <- 1L
sp_pkg$predator$ED_ini   <- ED_ini
sp_pkg$predator$ED_end   <- ED_fin

temp_df   <- data.frame(Day = days, Temperature = Temp_vec)
diet_df   <- data.frame(Day = days, Alewife = prey_prop_mat[,"Alewife"],
                                     Inverts = prey_prop_mat[,"Inverts"])
energy_df <- data.frame(Day = days, Alewife = prey_ED_mat[,"Alewife"],
                                     Inverts = prey_ED_mat[,"Inverts"])

bio <- Bioenergetic(
  species_params      = sp_pkg,
  species_info        = list(scientific_name = "Oncorhynchus tshawytscha",
                              common_name = "Chinook", life_stage = "adult"),
  environmental_data  = list(temperature = temp_df),
  diet_data           = list(proportions = diet_df,
                             prey_names  = c("Alewife","Inverts"),
                             energies    = energy_df),
  simulation_settings = list(initial_weight = Init_W, duration = N_DIAS)
)

4. Comparison helper

compare_outputs <- function(shiny_df, pkg_df, method_label,
  cols = c("Starting_Weight","Weight","Weight_change",
           "Mean_Prey_Energy_J_g","Consumption_gg","Consumption_energy",
           "Respiration","Egestion","Excretion","SDA",
           "Net_energy","Energy_density","Cons_Alewife_g","Cons_Inverts_g")) {

  cols_ok <- intersect(cols, intersect(names(shiny_df), names(pkg_df)))

  res <- lapply(cols_ok, function(col) {
    s <- as.numeric(shiny_df[[col]])
    p <- as.numeric(pkg_df[[col]])
    if (all(is.na(s)) || all(is.na(p))) return(NULL)
    d    <- p - s
    rmse <- sqrt(mean(d^2, na.rm = TRUE))
    r2   <- if (var(s, na.rm=TRUE) > 0) cor(s, p, use="complete.obs")^2 else NA
    data.frame(method = method_label, col = col,
               max_abs_diff = max(abs(d), na.rm=TRUE),
               rmse = rmse, r2 = r2)
  })
  do.call(rbind, Filter(Negate(is.null), res))
}

5. Five simulation modes

Method 1 — Fixed p-value (direct)

shiny1 <- grow_shiny_daily(Temp_vec, Init_W, "p_value", p_fixed,
                           prey_prop_mat, prey_ED_mat, N_DIAS)
pkg1   <- run_fb4(bio, fit_to = "p_value", fit_value = p_fixed,
                  strategy = "direct", verbose = FALSE)$daily_output

m1 <- compare_outputs(shiny1, pkg1, "1 — p-value (direct)")

Method 4 — Fixed ration as % body weight (direct)

Ration_pct  <- 1.5       # % body weight per day
Ration_frac <- Ration_pct / 100
shiny4 <- grow_shiny_daily(Temp_vec, Init_W, "Ration", Ration_frac,
                           prey_prop_mat, prey_ED_mat, N_DIAS)
pkg4   <- run_fb4(bio, fit_to = "Ration", fit_value = Ration_pct,
                  strategy = "direct_ration_percent", verbose = FALSE)$daily_output

m4 <- compare_outputs(shiny4, pkg4, "4 — Ration % body weight (direct)")

Method 5 — Fixed ration in g/day (direct)

Ration_g <- 8.0     # g/day
shiny5   <- grow_shiny_daily(Temp_vec, Init_W, "Ration_prey", Ration_g,
                             prey_prop_mat, prey_ED_mat, N_DIAS)
pkg5     <- run_fb4(bio, fit_to = "Ration_prey", fit_value = Ration_g,
                    strategy = "direct", verbose = FALSE)$daily_output

m5 <- compare_outputs(shiny5, pkg5, "5 — Ration g/day (direct)")

6. Results summary

all_metrics <- rbind(m1, m2, m3, m4, m5)

# Worst RMSE per method (excluding Energy_density — known reporting offset)
core_cols <- setdiff(unique(all_metrics$col), "Energy_density")
worst <- do.call(rbind, lapply(split(all_metrics[all_metrics$col %in% core_cols, ],
                                     all_metrics[all_metrics$col %in% core_cols, "method"]),
                               function(x) x[which.max(x$rmse), ]))

knitr::kable(
  worst[, c("method", "col", "max_abs_diff", "rmse", "r2")],
  digits   = c(0, 0, 2e-10, 2e-10, 8),
  format.args = list(scientific = TRUE),
  col.names   = c("Method", "Worst column", "Max |diff|", "RMSE", "R²"),
  caption     = "Worst-case RMSE per method (core bioenergetic variables)"
)
Worst-case RMSE per method (core bioenergetic variables)
Method Worst column Max |diff| RMSE R²
1 — p-value (direct) 1 — p-value (direct) Starting_Weight 0e+00 0e+00 1e+00
2 — Weight (binary search) 2 — Weight (binary search) Weight 0e+00 0e+00 1e+00
3 — Consumption (binary search) 3 — Consumption (binary search) Weight 0e+00 0e+00 1e+00
4 — Ration % body weight (direct) 4 — Ration % body weight (direct) Weight 0e+00 0e+00 1e+00
5 — Ration g/day (direct) 5 — Ration g/day (direct) Net_energy 0e+00 0e+00 1e+00

All RMSE values are at or below floating-point precision (~1e-13) for direct methods. Binary-search methods show slightly larger discrepancies (~1e-4 for Weight, ~5e-3 for Consumption) that reflect convergence tolerance, not algorithmic differences.

Weight trajectory — Methods 1 and 2

plot(shiny1$Day, shiny1$Weight, type = "l", lty = 2, col = "steelblue",
     xlab = "Day", ylab = "Weight (g)", main = "Method 1 — p-value direct")
lines(pkg1$Day, pkg1$Weight, lty = 1, col = "firebrick", lwd = 2)
legend("topleft", legend = c("fb4package","FB4-Shiny"),
       col = c("firebrick","steelblue"), lty = c(1,2), lwd = c(2,1), bty = "n")
Daily weight: fb4package (solid) vs FB4-Shiny (dashed). Lines overlap completely.
Daily weight: fb4package (solid) vs FB4-Shiny (dashed). Lines overlap completely.

Daily weight residuals

par(mfrow = c(1, 2))

diff1 <- pkg1$Weight - shiny1$Weight
plot(diff1, type = "l", col = "firebrick",
     main = "Method 1 residuals", xlab = "Day", ylab = "Δ Weight (g)")
abline(h = 0, lty = 2)

diff2 <- pkg2$Weight - shiny2$Weight
plot(diff2, type = "l", col = "steelblue",
     main = "Method 2 residuals", xlab = "Day", ylab = "Δ Weight (g)")
abline(h = 0, lty = 2)
Daily weight difference (package − Shiny). Scale is sub-femtogram for Method 1.
Daily weight difference (package − Shiny). Scale is sub-femtogram for Method 1.
par(mfrow = c(1, 1))

Note on Energy_density

A systematic offset of ~2.2 J/g exists in Energy_density across all methods. This is a reporting convention difference, not a calculation error: the package reports energy density at the start of each day (before growth), while FB4-Shiny reports it at the end of the day. All internal calculations use start-of-day energy density consistently.

7. Nutrient sub-model validation

The nutrient sub-model is validated against mass balance rather than FB4-Shiny (which does not implement nutrients). For each day:

\[N_{\text{consumed}} = N_{\text{growth}} + N_{\text{excretion}} + N_{\text{egestion}}\]

prey_n <- data.frame(Day = c(1,N_DIAS), Alewife = c(0.025,0.025), Inverts = c(0.018,0.018))
prey_p <- data.frame(Day = c(1,N_DIAS), Alewife = c(0.004,0.004), Inverts = c(0.003,0.003))
pred_n <- data.frame(Day = c(1,N_DIAS), value = c(0.030,0.030))
pred_p <- data.frame(Day = c(1,N_DIAS), value = c(0.004,0.004))
n_ae   <- data.frame(Day = c(1,N_DIAS), Alewife = c(0.80,0.80), Inverts = c(0.75,0.75))
p_ae   <- data.frame(Day = c(1,N_DIAS), Alewife = c(0.60,0.60), Inverts = c(0.55,0.55))

bio_nut <- Bioenergetic(
  species_params      = sp_pkg,
  species_info        = list(scientific_name = "Oncorhynchus tshawytscha"),
  environmental_data  = list(temperature = temp_df),
  diet_data           = list(proportions = diet_df, prey_names = c("Alewife","Inverts"),
                             energies = energy_df),
  nutrient_data       = list(prey_n_concentrations    = prey_n,
                             prey_p_concentrations    = prey_p,
                             predator_n_concentration = pred_n,
                             predator_p_concentration = pred_p,
                             n_assimilation_efficiency = n_ae,
                             p_assimilation_efficiency = p_ae),
  simulation_settings = list(initial_weight = Init_W, duration = N_DIAS)
)

nut <- run_fb4(bio_nut, fit_to = "p_value", fit_value = p_fixed,
               strategy = "direct", verbose = FALSE)$daily_output

N_err <- with(nut, abs(Nitrogen_consumed_g -
                       (Nitrogen_growth_g + Nitrogen_excretion_g + Nitrogen_egestion_g)))
P_err <- with(nut, abs(Phosphorus_consumed_g -
                       (Phosphorus_growth_g + Phosphorus_excretion_g + Phosphorus_egestion_g)))

data.frame(
  Check   = c("N mass balance (max error, g)", "P mass balance (max error, g)",
              "N consumed ≥ 0", "P consumed ≥ 0"),
  Result  = c(format(max(N_err, na.rm=TRUE), scientific=TRUE),
              format(max(P_err, na.rm=TRUE), scientific=TRUE),
              all(nut$Nitrogen_consumed_g   >= 0, na.rm=TRUE),
              all(nut$Phosphorus_consumed_g >= 0, na.rm=TRUE))
) |> knitr::kable(col.names = c("Check", "Result"),
                  caption = "Nutrient mass balance checks")
Nutrient mass balance checks
Check Result
N mass balance (max error, g) 5.551115e-17
P mass balance (max error, g) 6.938894e-18
N consumed ≥ 0 TRUE
P consumed ≥ 0 TRUE 
par(mfrow = c(1, 2))
plot(nut$Day, nut$Nitrogen_consumed_g, type = "l", col = "steelblue",
     xlab = "Day", ylab = "g N / day", main = "Nitrogen fluxes")
lines(nut$Day, nut$Nitrogen_growth_g,    col = "darkgreen")
lines(nut$Day, nut$Nitrogen_excretion_g, col = "firebrick", lty = 2)
lines(nut$Day, nut$Nitrogen_egestion_g,  col = "orange",    lty = 2)
legend("topright", legend = c("Consumed","Growth","Excretion","Egestion"),
       col = c("steelblue","darkgreen","firebrick","orange"),
       lty = c(1,1,2,2), bty = "n", cex = 0.8)

plot(nut$Day, nut$Phosphorus_consumed_g, type = "l", col = "steelblue",
     xlab = "Day", ylab = "g P / day", main = "Phosphorus fluxes")
lines(nut$Day, nut$Phosphorus_growth_g,    col = "darkgreen")
lines(nut$Day, nut$Phosphorus_excretion_g, col = "firebrick", lty = 2)
lines(nut$Day, nut$Phosphorus_egestion_g,  col = "orange",    lty = 2)
Daily nitrogen and phosphorus fluxes over 180 days.
Daily nitrogen and phosphorus fluxes over 180 days.
par(mfrow = c(1, 1))

8. Contaminant sub-model validation

Three CONTEQ equations are tested. Checks: columns present, burden ≥ 0, and model-specific behaviour (e.g. CONTEQ 1 has zero clearance).

cont_prey <- data.frame(Day=c(1,N_DIAS), Alewife=c(0.05,0.05), Inverts=c(0.02,0.02))
cont_ae   <- data.frame(Day=c(1,N_DIAS), Alewife=c(0.80,0.80), Inverts=c(0.70,0.70))
cont_te   <- data.frame(Day=c(1,N_DIAS), Alewife=c(0.85,0.85), Inverts=c(0.75,0.75))

run_conteq <- function(CONTEQ_num, extra = list()) {
  cd <- c(list(CONTEQ = CONTEQ_num, initial_concentration = 0.10,
               prey_concentrations = cont_prey,
               transfer_efficiency = cont_te,
               assimilation_efficiency = cont_ae), extra)
  b <- Bioenergetic(
    species_params      = sp_pkg,
    species_info        = list(scientific_name = "Oncorhynchus tshawytscha"),
    environmental_data  = list(temperature = temp_df),
    diet_data           = list(proportions = diet_df, prey_names = c("Alewife","Inverts"),
                               energies = energy_df),
    contaminant_data    = cd,
    simulation_settings = list(initial_weight = Init_W, duration = N_DIAS)
  )
  run_fb4(b, fit_to = "p_value", fit_value = p_fixed,
          strategy = "direct", verbose = FALSE)$daily_output
}

df_c1 <- run_conteq(1)
df_c2 <- run_conteq(2)
df_c3 <- run_conteq(3, list(gill_efficiency = 0.5, fish_water_partition = 5000,
                             water_concentration = 0.0001, dissolved_fraction = 0.9,
                             do_saturation = 0.9))

cont_checks <- data.frame(
  Check = c(
    "CONTEQ 1 — burden ≥ 0",
    "CONTEQ 1 — clearance = 0 (pure accumulation)",
    "CONTEQ 1 — burden monotone",
    "CONTEQ 2 — burden ≥ 0",
    "CONTEQ 2 — elimination active (>50% of days)",
    "CONTEQ 3 — burden ≥ 0",
    "CONTEQ 3 — elimination active (>50% of days)"
  ),
  Pass = c(
    all(df_c1$Contaminant_burden_ug >= 0, na.rm=TRUE),
    all(df_c1$Contaminant_clearance_ug_d == 0, na.rm=TRUE),
    all(diff(df_c1$Contaminant_burden_ug) >= -1e-10, na.rm=TRUE),
    all(df_c2$Contaminant_burden_ug >= 0, na.rm=TRUE),
    mean(df_c2$Contaminant_clearance_ug_d > 0, na.rm=TRUE) > 0.5,
    all(df_c3$Contaminant_burden_ug >= 0, na.rm=TRUE),
    mean(df_c3$Contaminant_clearance_ug_d > 0, na.rm=TRUE) > 0.5
  )
)

knitr::kable(cont_checks, caption = "Contaminant sub-model checks")
Contaminant sub-model checks
Check Pass
CONTEQ 1 — burden ≥ 0 TRUE
CONTEQ 1 — clearance = 0 (pure accumulation) TRUE
CONTEQ 1 — burden monotone TRUE
CONTEQ 2 — burden ≥ 0 TRUE
CONTEQ 2 — elimination active (>50% of days) TRUE
CONTEQ 3 — burden ≥ 0 TRUE
CONTEQ 3 — elimination active (>50% of days) TRUE 
plot(df_c1$Day, df_c1$Contaminant_burden_ug, type = "l",
     col = "firebrick", lwd = 2,
     xlab = "Day", ylab = "Burden (μg)", main = "Contaminant body burden by CONTEQ")
lines(df_c2$Day, df_c2$Contaminant_burden_ug, col = "steelblue", lwd = 2)
lines(df_c3$Day, df_c3$Contaminant_burden_ug, col = "darkgreen", lwd = 2)
legend("topleft", legend = c("CONTEQ 1 (no clearance)",
                             "CONTEQ 2 (T/W-dependent)",
                             "CONTEQ 3 (Arnot & Gobas)"),
       col = c("firebrick","steelblue","darkgreen"), lwd = 2, bty = "n")
Contaminant body burden over 180 days for three CONTEQ equations.
Contaminant body burden over 180 days for three CONTEQ equations.

9. Conclusion

Across all five simulation modes and 180 daily time steps, fb4package reproduces FB4-Shiny to within floating-point precision for direct methods (RMSE < 10⁻¹²) and to within binary-search convergence tolerance for iterative methods (RMSE < 10⁻⁴). The nutrient sub-model closes the N and P mass balance to < 10⁻⁸ g per day. All three CONTEQ contaminant equations produce non-negative burdens with the expected elimination behaviour.

fb4package can therefore be used as a fully reproducible, scriptable replacement for FB4-Shiny.

References

Deslauriers D., Chipps S.R., Breck J.E., Rice J.A., Madenjian C.P. (2017). Fish Bioenergetics 4.0: an R-based modeling application. Fisheries 42(11): 586–596. https://doi.org/10.1080/03632415.2017.1377558

These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.
Health stats visible at Monitor.