Hi Staners,
A while ago I posted about “Speeding a bioenergetic ODE model.” Thanks again for all the helpful feedback. Since then, we published a paper (paper) and now have code that routinely estimates DEB parameters for thousands of individual fish in about an hour.
That speed-up relied on assuming reproduction has a negligible effect on wet‐weight dynamics—reasonable for juveniles, but not adults. We now want to explicitly include reproduction dynamics because they matter both practically and theoretically. For each fish we monitor several times (repeated measures) length, wet weight, gonad wet weight, and the full oocyte diameter distribution.
Model idea
Energy from a reproductive buffer (as modelled by DEB; dynamic energy budget model) is mobilized into the gonad; gonad size emerges as the sum of the volumes/energies of all growing oocytes.
Assumptions
-
Off-season oocyte sizes are normally distributed with known μ and σ.
-
Potential fecundity (Fec)—the number of oocytes that will grow in the coming season—is fixed at the start of the season.
-
Oocyte growth in energy: dEdt=p E^γ. Energy maps to diameter via known constants.
-
p (growth supply) is either p_max (if the buffer can cover all oocyte demand) or is downregulated by current buffer energy (i.e., better food/temperature → larger buffer → faster growth).
-
Growth only within the reproductive season.
-
Upper size cap: oocytes stop growing at a known maximum diameter and are then spawned.
-
Discretization: we divide the diameter range into bins of dD=20 μm and evolve each bin with rules 2–6.
-
Reset cost negligible: energy to restore the off-season size distribution after spawning is assumed negligible.
Parameters of interest for reproduction are p_max, γ, Fec, and the start date of the reproductive season, with a focus on between-fish variability.
Implementation + issue
We coded this in Stan (by now, independent analysis for each fish; no paralelization; cmdstanr version 0.5.3; CmdStan version: 2.30.1), and simulated input data for length, wet weight, gonad weight, and oocyte size distributions using deSolve in R. The Stan model seems to recover patterns well (fish wet weight as example:
), but it’s extremely slow: ~1 hour for 20 iterations with several free parameters; ~10 minutes even when all but one parameter are fixed. Stiff vs. non-stiff ODE solvers behave similarly. Splitting the year into “before/after season” requires two integrals and makes it even slower. The 40-bin “brute force” discretization is likely a culprit, but we don’t have the math background to replace it with a smoother formulation.
Ask: Any advice on (a) math to avoid binning, and/or (b) Stan coding tricks to make this tractable would be hugely appreciated.
Simulated input data:
input_stan.RData (119.7 KB)
Code:
#---------------------------------------
# Fish oocyte dynamics
# Tomas and Palmer
#---------------------------------------
# Last update: Sep 2025
# -----
# 1: Loading data and libraries
# -----
remove(list=ls())
library(cmdstanr)
load("input_stan.RData")
#ls()
# [1] "age" "fish_total_length" "fish_wet_weight"
# [4] "gonad_wet_weight" "n_hat" "parms2"
dim=dim(n_hat)
# adding observation error
sd_length=0.1 # cm
sd_weight=1 # gr
sd_gonad_weight=0.1 # gr
temp=cbind(fish_total_length+rnorm(dim[1],0,sd_length),
fish_wet_weight+rnorm(dim[1],0,sd_weight),
gonad_wet_weight+rnorm(dim[1],0,sd_gonad_weight)
)
#splitting first observation from the others
# t0 for integration is set to the first observation
t0=age[1]
age=age[2:dim[1]]
obs0=temp[1,]
obs=temp[2:dim[1],]
n_obs=n_hat[2:dim[1],] # range(n_obs)
n_obs=cbind(n_obs,rep(0,dim(n_obs)[1])) # the last column represents the oocytes at the spawning diemeter
n=length(age)
parms2$n_bins=length(parms2$n_seq) # number of bins
#-----
# 2: Analysis
#-----
#----
# 2.1: STAN model
#----
sink("DEB.stan")
cat(" // first line
functions {
/////////////////////
// DEB function
/////////////////////
vector DEB(
real t, vector y,
// core parameters
real f, real E_G, real kap, real p_Am, real p_M, real k_J, real v, real kap_R, real Hp,
// temperature parameters
real Kelvin, real Temp_mean, real Temp_amp, real pi2f, real Temp_phi,
real T_A, real T_AL, real T_L, real T_AH, real T_H, real T_ref, real birth_correction, real sT_ref,
// reproduction parameters
int n_bins, vector E_seq, vector n_seq,
real Fec, real p_max, real gamma, real age0, real E_max
) {
vector[5] dydt; // derivatives: [E, V, p_sum, B, G]
// --- Precomputing repeated calculations ---
real V = y[2];
real V23 = pow(V, 2.0/3.0); // surface
vector[n_bins] E_seq_pow = E_seq^(1-gamma); // preciompute power
// --- Temperature correction ---
real Temp = Kelvin + Temp_mean + Temp_amp * sin(pi2f * (t + birth_correction) + Temp_phi);
real sT = 1 + exp(T_AL/Temp - T_AL/T_L) + exp(T_AH/T_H - T_AH/Temp);
//real sT_ref = 1 + exp(T_AL/T_ref - T_AL/T_L) + exp(T_AH/T_H - T_AH/T_ref); // moved to trasnformed data block
real cT = (sT_ref / sT) * exp(T_A/T_ref - T_A/Temp);
// Temperature-corrected parameters
real p_Am_T = cT * p_Am;
real v_T = cT * v;
real p_M_T = cT * p_M;
real k_J_T = cT * k_J;
// --- Fluxes ---
real pA = f * p_Am_T * V23; // assimilation
real pS = p_M_T * V; // maintenance
real pC = (y[1]/V) * ((v_T * E_G * V23 + pS) / (E_G + kap * y[1]/V)); // mobilization
real pG = kap * pC - pS; // growth
real pJ = k_J_T * Hp; // maturity maintenance (adult)
real pR = (1 - kap) * pC - pJ; // reproduction
// --- Reproduction ---
// Rationale and asumtions:
// 1: Oocite (future eggs) growth (with energy as common currency) is modelled as dE/dt=p*E^gamma
// 2: p can be either p_max (when the reproductive buffer is able to cover the growth of all the oocites),
// or is recalculed depending of the current energy at the reproductive buffer).
// 3: Oocytes do not grow when they reach a (known) maximum diameter; Then, they are spawn.
// 4: Oocyte diameter distribution during the non reproductive season is normal with known mu and sigma
// 5: Provided (4), the full range of possible oocyte diamters is divided into small bins (dD=20 microns)
// and the growth of each bin is monitored with (1).
// 6: The number of oocytes that will consume energy at a given reproductive season (Potential fecundity, Fec)
// is fixed at the start of the reproductive season.
// 7: The energy needed to recover the initial state of oocyte diameter distribution at the end of a
// reproductive season is considered negligible
real p_new = 0.0; real pEo = 0.0; real spawn_dt = 0.0;
vector[n_bins] spawn_old; // oocytes spawned till t
vector[n_bins] spawn_new; // oocytes sapwned till t+dt
real c_season = 1.0/(1.0 + exp(-0.1*(t-age0))); // reproductive vs non reproductive season (arbitrary slope)
real p_max_s = c_season*p_max; // p_max_s will be zero at non reproductive season
// compute Ei and Ei_old (energy per bin) (y[3] is the acumulated p from t=0)
vector[n_bins] E_i = pow(E_seq_pow + (y[3] + p_max_s)*(1-gamma), 1.0/(1-gamma));
vector[n_bins] E_i_old = pow(E_seq_pow + y[3]*(1-gamma), 1.0/(1-gamma));
// pE_i (energy needed per bin at the current dt)
vector[n_bins] pE_i = 0.5 * p_max_s * (E_i^gamma + E_i_old^gamma);
// spawning (needed for modelling the energy lost)
for (i in 1:n_bins){
spawn_old[i] = (E_i_old[i] < E_max ? 1 : 0); // spawning till t
spawn_new[i] = (E_i[i] < E_max ? 1 : 0); // spawing till t+dt when p=p_max
}
// first pass
// pEo is the total energy needed when p=p_max (the reproductive buffer can porvide the energy nedded)
pEo = 0.5 * Fec * dot_product(pE_i, n_seq .* (spawn_old + spawn_new));
real temporal = 0.5 * Fec * sum((E_i_old^gamma .* n_seq .* spawn_old) +
(E_i^gamma .* n_seq .* spawn_new));
// Comparing energy needed and energy available (y[4] is the energy at the reproductive buffer)
pEo = fmin(pEo, fmax(0, y[4]));
// Recalculating p for the cases when the energy at the buffer cannot cover the needs
if (temporal > 0.0){
p_new = fmin(p_max, fmax(0, pEo/temporal));
}
// second pass (recalculating first pass with p=p:new)
vector[n_bins] E_i_new = pow(E_seq_pow + (y[3] + p_new)*(1-gamma), 1.0/(1-gamma));
for (i in 1:n_bins){
spawn_dt += n_seq[i]*Fec * ((E_i_new[i] >= E_max) - (E_i_old[i] >= E_max));
}
// --- Derivatives ---
dydt[1] = pA - pC;
dydt[2] = pG / E_G;
dydt[3] = p_new;
dydt[4] = kap_R * pR - pEo;
dydt[5] = pEo - spawn_dt * E_max;
return dydt;
} // end of DEB function
//////////////////////
/////////////////////
// interpolation function (based on dmi3kno from stan forum)
/////////////////////
real linear_interpolation(real D_seq, vector D_i, vector n_seq_new, real th_max){
int K = rows(D_i);
vector[K] deltas = D_seq - D_i;
real ans = 0.0;
if(D_seq>D_i[1] && D_seq<th_max){
int i;
i = sort_indices_asc(abs(deltas))[1];//this is equivalent to which.min()
if(deltas[i]<=0) i -= 1;
ans = n_seq_new[i] + (n_seq_new[i + 1] - n_seq_new[i]) * (D_seq - D_i[i])/(D_i[i + 1] - D_i[i]);
}
return ans;
}
//////////////////////
/////////////////////
// oocyte diameter function (moving from energy to diameter)
/////////////////////
array [] vector diameter(
int n, array [] vector y,
int n_bins, vector E_seq, vector n_seq, vector D_seq,
real gamma, real th_max, real dD,real w_E, real mu_E, real d_egg
){
// diameter distribution (function output)
array [n] vector [n_bins] n_hat; // expected number of oocytsr per diameter bin in D_seq
// rows are sampling events; columns are diamters bins
vector[n_bins] E_seq_pow = E_seq^(1-gamma); // precompute power
for (n_event in 1:n){ // loop for each sampling event
// compute Ei
vector[n_bins] E_i = pow((E_seq_pow + y[n_event,3]*(1-gamma)), 1.0/(1.0-gamma));
// Diameters
vector [n_bins] D_i = pow((E_i*w_E/mu_E/d_egg/pi()*6e12),(1.0/3.0));
//vector [n_bins-1] n_seq_new = dD*(n_seq[1:(n_bins-1)]./(D_i[2:n_bins]-D_i[1:(n_bins-1)]));
vector [n_bins] n_seq_new;
n_seq_new [1:(n_bins-1)] = dD*(n_seq[1:(n_bins-1)]./(D_i[2:n_bins]-D_i[1:(n_bins-1)]));
n_seq_new[n_bins]=0.0;
// interpolation
for (i in 1:n_bins){
n_hat[n_event,i]=linear_interpolation(D_seq[i], D_i[], n_seq_new[], th_max);
}
} // end sampling events loop
return n_hat;
} // end of diameter function
//////////////////////
}//end bloc functions
data {
int n; // number of observations per fish
array [n,3] real obs; // observations; repeated measures (1:n),c(fish length,fish weight,gonad weight)
real t0; // initial time
array [n] real age; // times at which the system is observed
array [3] real obs0; // observations at t0; c(length,weight,gonad)
// fixed parameters (=input) for the integrate function (ODE)
real f;
//real p_Am;
real v;
real E_G;
real kap;
real p_M;
real k_J;
real kap_R;
real Hp;
// other parameters
real del_M;
real w_E;
real mu_E;
real d_V; // water fraction of body
// temperature parameters
real Kelvin;
real Temp_mean;
real Temp_amp;
real pi2f;
real Temp_phi;
real T_A;
real T_AL;
real T_L;
real T_AH;
real T_H;
real T_ref;
real birth_correction;
// reproduction parameters
real d_egg; // water fraction of eggs
int n_bins;
vector [n_bins] E_seq;
vector [n_bins] D_seq;
vector [n_bins] n_seq;
real E_max;
real th_max;
array [n,n_bins] int n_obs; // observed number of oocytes per bin
real dD;
real Fec;
real p_max;
real gamma;
real age0;
real y0_E;
real y0_V;
real y0_p_sum;
real y0_B;
real y0_G;
// priors
array [2] real prior_V0;
array [2] real prior_E0;
array [2] real prior_pAm;
array [2] real prior_v;
array [2] real prior_kap;
array [2] real prior_Fec;
array [2] real prior_p_max;
array [2] real prior_gamma;
array [2] real prior_age0;
array [2] real prior_sd_length;
array [2] real prior_sd_weight;
array [2] real prior_sd_gonad_weight;
}
transformed data {
real sT_ref = 1 + exp(T_AL/T_ref - T_AL/T_L) + exp(T_AH/T_H - T_AH/T_ref);
}
parameters {
// DEB parameters
real <lower=0> p_Am;
//real <lower=0> v;
//real <lower=0, upper=1> kap;
// reproduction parameters
//real <lower=0> Fec;
//real <lower=0> p_max;
//real <lower=0.75, upper=0.999> gamma;
//real <lower=400,upper=600> age0;
// initial state
//real <lower=0> y0_E;
//real <lower=0> y0_V;
// observation error
real <lower=0> sd_length;
real <lower=0> sd_weight;
real <lower=0> sd_gonad_weight;
}
transformed parameters {
}
model {
array [n] vector[5] y; // state variables (output numerical integration DEB function)
array [n] vector[n_bins] n_hat; // expected number of oocytes per bin (output of diamter function)
vector [n_bins] prob; // all bins
vector [5] y0; //V,E,p_cum,B,G (caution: G amd p_cum is 0.0 at non-reproduction season)
// DEB paramters
p_Am ~ normal(prior_pAm[1] , prior_pAm[2]);
//v ~ normal(prior_v[1] , prior_v[2]);
//kap ~ normal(prior_kap[1] , prior_kap[2]); //beta(prior_kap[1] , prior_kap[2]);
// reproduction parameters
//Fec ~ normal(prior_Fec[1] , prior_Fec[2]);
//p_max ~ normal(prior_p_max[1] , prior_p_max[2]);
//gamma ~ normal(prior_gamma[1] , prior_gamma[2]);
//age0 ~ normal(prior_age0[1] , prior_age0[2]);
// state variables at t0
//y0_E ~ normal(prior_E0[1], prior_E0[2]); // reserve energy
//y0_V ~ normal(prior_V0[1], prior_V0[2]); // structural volume
//y0[3] = 0.0; // p_sum at non-reproductive season CAUTION!!!
//y0[4] ~ normal(prior_B0[1], prior_B0[2]); // buffer energy CAUTION!!!
//y0[5] = 0.0; // gonad energy //t0 is at non-reproductive season CAUTION!!!
y0[1]=y0_E; // estimated
y0[2]=y0_V; // estimated
y0[3]=y0_p_sum; // 0.0 (t0 is at non reproductive season)
y0[4]=y0_B; // tentativelly asumed to be zero at t0
y0[5]=y0_G; // tentativelly asumed to be zero at t0
// measurement errors
sd_length ~ normal(prior_sd_length[1],prior_sd_length[2]);
sd_weight ~ normal(prior_sd_weight[1],prior_sd_weight[2]);
sd_gonad_weight ~ normal(prior_sd_gonad_weight[1],prior_sd_gonad_weight[2]);
obs0[1] ~ normal( ((y0[2]^(1.0/3.0))/del_M), sd_length); // length at t0
obs0[2] ~ normal( y0[2] + // structure weight at t0
(y0[1]+y0[4])*w_E/mu_E/d_V + // energy weight at t0
y0[5]*w_E/mu_E/d_egg, sd_weight); // gonad weight at t0
obs0[3] ~ normal( y0[5]*w_E/mu_E/d_egg, sd_gonad_weight); // gonad weight at t0
// numerical integration
y[1:n] = ode_rk45 (DEB, y0[], t0 , age[], //ode_bdf ode_rk45
// DEB main parameters
f,E_G,kap,p_Am,p_M,k_J,v,kap_R,Hp,
// temperature-related paramters
Kelvin,Temp_mean,Temp_amp,pi2f,Temp_phi,T_A,T_AL,T_L,T_AH,T_H,T_ref,birth_correction,sT_ref,
// reproduction
n_bins,E_seq[],n_seq[],Fec,p_max,gamma,age0,E_max
);
// estimating diameter distribution
n_hat[1:n] = diameter(n,y[],
n_bins,E_seq[],n_seq[],D_seq[],gamma,
th_max,dD,w_E,mu_E,d_egg
);
// likelihood
// 1: dEdt
// 2: dVdt
// 3: dp_sumdt
// 4: dBdt
// 5: dGdt
for (i in 1:n){
obs[i,1] ~ normal(pow(y[i,2],(1.0/3.0))/del_M,sd_length); //fish total length
obs[i,2] ~ normal(y[i,2] + // structure
(y[i,1]+y[i,4])*w_E/mu_E/d_V+ // reserve and buffer
y[i,5]*w_E/mu_E/d_egg, // gonad
sd_weight); // fish wet weight
obs[i,3] ~ normal(y[i,5]*w_E/mu_E/d_egg,sd_gonad_weight); //gonad weight
prob = to_vector(n_hat[i,1:n_bins]);
n_obs[i,1:n_bins] ~ multinomial(softmax(prob[]));
}
}
generated quantities {
array [n] vector[5] y_pre; // state variables (output numerical integration DEB function)
array [n] vector[n_bins] n_pre; // expected number of oocytes per bin (output of diamter function)
vector [5] y0_pre; //E,V,p_cum,B,G (caution: G is 0.0 at non-reproduction season)
y0_pre[1]=y0_E;
y0_pre[2]=y0_V;
y0_pre[3]=y0_p_sum;
y0_pre[4]=y0_B;
y0_pre[5]=y0_G;
y_pre[1:n] = ode_rk45 (DEB, y0_pre[], t0 , age[], //ode_bdf ode_rk45
// DEB main parameters
f,E_G,kap,p_Am,p_M,k_J,v,kap_R,Hp,
// temperature-related paramters
Kelvin,Temp_mean,Temp_amp,pi2f,Temp_phi,T_A,T_AL,T_L,T_AH,T_H,T_ref,birth_correction,sT_ref,
// reproduction
n_bins,E_seq[],n_seq[],Fec,p_max,gamma,age0,E_max
);
n_pre[1:n] = diameter(n,y_pre[],
n_bins,E_seq[],n_seq[],D_seq[],gamma,
th_max,dD,w_E,mu_E,d_egg
);
}
" # end of model code
,fill = TRUE)
sink()
#-----
# 2.2: compiling model
#-----
mod = cmdstan_model("DEB.stan")
#mod = cmdstan_model("DEB.stan",cpp_options = list(stan_threads = TRUE))
#-----
# 2.3: initializating model
#-----
n.chains = 2
initializer = function() list(
"p_Am"=parms2$p_Am,
#"v"=parms2$v,
#"kap"=parms2$kap,
#"Fec"=parms2$Fec,
#"p_max"=parms2$p_max,
#"gamma"=parms2$gamma,
#"age0"= parms2$t0_reproduction,
#"y0_E"=c((obs0[2]-(obs0[1]*parms2$del_M)^3)/(parms2$w_E/parms2$mu_E/parms2$d_V)),
#"y0_V"=(obs0[1]*parms2$del_M)^3,
"sd_length"=sd_length,
"sd_weight"=sd_weight,
"sd_gonad_weight"=sd_gonad_weight
)
inits = list()
for (chain in 1:n.chains) inits[[chain]] = initializer()
#------
# 2.4: priors
#------
prior_V0=c((obs0[1]*parms2$del_M)^3,
0.5*(obs0[1]*parms2$del_M)^3
)
prior_E0=c((obs0[2]-(obs0[1]*parms2$del_M)^3)/(parms2$w_E/parms2$mu_E/parms2$d_V),
0.5*(obs0[2]-(obs0[1]*parms2$del_M)^3)/(parms2$w_E/parms2$mu_E/parms2$d_V)
)
prior_pAm=c(parms2$p_Am,0.5*parms2$p_Am)
prior_v=c(parms2$v,0.5*parms2$v)
prior_kap=c(parms2$kap,0.1*parms2$kap)
prior_Fec=c(parms2$Fec,0.5*parms2$Fec)
prior_p_max=c(parms2$p_max,0.01*parms2$p_max)
prior_gamma=c(parms2$gamma,0.01*parms2$gamma)
prior_age0=c(parms2$t0_reproduction,0.5*parms2$t0_reproduction)
prior_sd_length=c(sd_length,0.5*sd_length)
prior_sd_weight=c(sd_weight,0.5*sd_weight)
prior_sd_gonad_weight=c(sd_gonad_weight,0.5*sd_gonad_weight)
#-----
# 2.5: running
#-----
fit = mod$sample(
data =list (
n = n, # number of observations per fish (number sampling events)
obs = obs, # observations
obs0 = obs0,
t0 = t0, # initial time
age = age, # times at which the system is observed (sampling dates)
n_obs = n_obs, # number of oocyte diameters at t0
# DEB main paramteters
f = parms2$f, # functional response
k_J = parms2$k_J,
kap_R = parms2$kap_R,
E_G = parms2$E_G,
p_M = parms2$p_M,
Hp = parms2$E_Hp,
#p_Am=parms2$p_Am,
v=parms2$v,
kap=parms2$kap,
y0_E=c((obs0[2]-(obs0[1]*parms2$del_M)^3)/(parms2$w_E/parms2$mu_E/parms2$d_V)),
y0_V=(obs0[1]*parms2$del_M)^3,
y0_p_sum=0.0,
y0_B=0.0,
y0_G=0.0,
#other parameters
del_M = parms2$del_M,
w_E = parms2$w_E,
mu_E = parms2$mu_E,
d_V = parms2$d_V,
# Temperature
Kelvin=parms2$Kelvin,
Temp_mean=parms2$Temp_mean,
Temp_amp=parms2$Temp_amp,
pi2f=parms2$pi2f,
Temp_phi=parms2$Temp_phi,
T_A=parms2$T_A,
T_AL=parms2$T_AL,
T_L=parms2$T_L,
T_AH=parms2$T_AH,
T_H=parms2$T_H,
T_ref=parms2$T_ref,
birth_correction=parms2$birth_correction,
# reproduction paramters
d_egg=parms2$d_egg,
n_bins=parms2$n_bins,
E_seq=parms2$E_seq,
D_seq=parms2$D_seq,
n_seq=parms2$n_seq,
E_max=parms2$E_max,
th_max=parms2$th_max,
dD=parms2$dD,
Fec=parms2$Fec,
p_max=parms2$p_max,
gamma=parms2$gamma,
age0=parms2$t0_reproduction,
#priors
prior_V0=prior_V0,
prior_E0=prior_E0,
prior_pAm=prior_pAm,
prior_v=prior_v,
prior_kap=prior_kap,
prior_Fec=prior_Fec,
prior_p_max=prior_p_max,
prior_gamma=prior_gamma,
prior_age0=prior_age0,
prior_sd_length=prior_sd_length,
prior_sd_weight=prior_sd_weight,
prior_sd_gonad_weight=prior_sd_gonad_weight
),
chains = n.chains,
parallel_chains = n.chains,
#threads_per_chain = 8,
iter_warmup = 10,
iter_sampling = 10,
init = inits,
max_treedepth = 10,
adapt_delta = 0.8,
refresh = 1#,
#fixed_param = TRUE
)