September 3, 2025.

Equation (8) should read $\mathrm{d\lambda }$₁/dt = (λdΓ/dt - λ₁)kτ^-1 to be consistent with (7). If $\mathrm{\Gamma (t)}$ were differentiable, this follows by writing (7) as $k$ nested convolutions of exponential densities with rate $\mathrm{k/\tau }$, and applying the Leibniz rule for integration under the integral sign. Since $\mathrm{\Gamma (t)}$ is a jump process, $\mathrm{d\Gamma /dt}$ is interpreted via a sum of Dirac delta functions at jump times, but the argument remains essentially the same.

The code is implemented correctly to be consistent with Equation (7). Equations (8) and (9) are represented in the code with nT=2 stages as

dT1dt = (lambda * dgammadt-T1)*nT/tau;
dmudt = (T1-mu)*nT/tau;


The numerical solution is coded to make an Euler approximation to the Gamma noise,

dgammadt=rgamma(dt/(sigSE*sigSE),(sigSE*sigSE)/dt);


We thank Hannes Kühnert for pointing out this issue.