Bayesian Parameter Balancing Enables Robust and Consistent Estimation of Kinetic Parameter Uncertainty

We study two benchmark pathways, glycolysis and the TCA cycle, because both are well characterized and contain sufficient thermodynamic and kinetic information for controlled parameter balancing experiments. For each pathway, we constructed a thermodynamically consistent reference parameter set by combining curated literature values from pathway-specific kinetic datasets with equilibrium information from eQuilibrator (Flamholz et al., 2012; Noor et al., 2014; Beber et al., 2022). The resulting reference values, shown in Figures 1 and 2 (blue bars), serve as the baseline parameter sets for downstream simulation and validation.

The kinetic context is a metabolic reaction network with $N{}_{\mathrm{m}\mathrm{e}\mathrm{t}}$ metabolites and $N{}_{\mathrm{r}\mathrm{x}\mathrm{n}}$ reactions. Let $𝑿\in ℝ{}^{N{}_{\mathrm{m}\mathrm{e}\mathrm{t}}}$ denote the metabolite concentration vector, $𝑺\in ℝ{}^{N{}_{\mathrm{m}\mathrm{e}\mathrm{t}}\times N{}_{\mathrm{r}\mathrm{x}\mathrm{n}}}$ the stoichiometric matrix, and $𝒗(𝑿,𝜽)\in ℝ{}^{N{}_{\mathrm{r}\mathrm{x}\mathrm{n}}}$ the reaction-rate vector, where $𝜽\in ℝ{}^P$ collects the kinetic parameters governing the rates. The corresponding dynamics can be written as

$$\frac{d𝑿}{dt}=𝑺𝒗(𝑿,𝜽).$$

In typical biochemical models, the entries of $𝜽$ include catalytic rate constants, Michaelis constants, enzyme abundances, metabolite concentrations, and thermodynamic quantities such as standard chemical potentials or equilibrium constants. Although our inference does not fit the ordinary differential equation dynamics directly, this mechanistic setting determines which parameters are linked and which physical constraints must be respected.

Following the classical parameter balancing literature (Liebermeister et al., 2010; Lubitz et al., 2010; Lubitz and Liebermeister, 2019), we work on a transformed scale on which the relevant parameter dependencies become linear. We denote the full collection of transformed model quantities by $𝒙\in ℝ{}^M,$ and a lower-dimensional set of independent basis quantities by $𝒒\in ℝ{}^D,D\le M.$ These are linked through the dependency map

$𝒙=𝑸𝒒,$

(1)

where $𝑸\in ℝ{}^{M\times D}$ is determined by the pathway structure and the chosen thermodynamic parameterization. Observed quantities are extracted through the corresponding observation map $𝑸{}_{\mathrm{o}\mathrm{b}\mathrm{s}}\in ℝ{}^{M{}_{\mathrm{o}\mathrm{b}\mathrm{s}}\times D},$ yielding the curated data vector $𝒙{}^{\ast }\in ℝ{}^{M{}_{\mathrm{o}\mathrm{b}\mathrm{s}}}.$ In practice, $𝑸{}_{\mathrm{o}\mathrm{b}\mathrm{s}}$ is obtained from $𝑸$ by selecting, deleting, and, when needed, duplicating rows to align the dependency structure with the available observations. The associated observation-noise covariance is $𝑪{}_{\mathrm{o}\mathrm{b}\mathrm{s}}\in ℝ{}^{M{}_{\mathrm{o}\mathrm{b}\mathrm{s}}\times M{}_{\mathrm{o}\mathrm{b}\mathrm{s}}}.$ This construction ensures that inferred parameter sets remain thermodynamically consistent by design.

The role of $𝑸$ is to encode biochemical relationships among kinetic and thermodynamic quantities. Forward and reverse catalytic constants, Michaelis constants, and equilibrium constants cannot in general be chosen independently; they must satisfy thermodynamic identities such as Haldane relationships, which link kinetics to equilibrium thermodynamics. As a result, naive assembly of literature values can produce parameter sets that are numerically plausible in isolation but physically incompatible when combined.

Parameter balancing resolves this by reparameterizing the model in terms of independent basis quantities $𝒒$, from which all derived quantities $𝒙$ are obtained through $𝑸$. In practice, $𝒒$ includes transformed rate, affinity, binding, and thermodynamic terms, while $𝒙$ contains the full set of kinetic and thermodynamic quantities of interest. After transformation, these dependencies are linear, which makes it possible to formulate parameter balancing as a constrained statistical inference problem while preserving biochemical plausibility. Additional kinetic background, including an illustrative reversible reaction, convenience kinetics, and the induced Haldane constraint, is provided in Supplementary Sections C.2, C.3 and C.4.

We formulate parameter balancing as Bayesian inference on the basis parameters $𝒒$. In classical parameter balancing, one first assembles a curated data vector $𝒙{}^{\ast }$ containing all available numerical information about the model quantities, including literature values for $K{}_M$, $k{}_{\mathrm{c}\mathrm{a}\mathrm{t}}$, metabolite concentrations, and equilibrium-related quantities. Some entries of $𝒙{}^{\ast }$ may be missing, while others may summarize multiple measurements with associated uncertainty. The aim is to infer a complete, self-consistent parameter vector $𝒙$ that matches the observed entries as closely as possible while respecting the structural constraints encoded by $𝑸$.

The forward model is defined as in Equation 1 and the observed quantities are linked to the basis parameters through the observation map $𝑸{}_{\mathrm{o}\mathrm{b}\mathrm{s}}$. The Gaussian baseline model is

$𝒙{}^{\ast }\mid 𝒒$

$\sim 𝒩(𝑸{}_{\mathrm{o}\mathrm{b}\mathrm{s}}𝒒,𝑪{}_{\mathrm{o}\mathrm{b}\mathrm{s}}),𝒒\sim 𝒩(\overline{𝒒}{}_{\mathrm{p}\mathrm{r}\mathrm{i}\mathrm{o}\mathrm{r}},𝑪{}_{\mathrm{p}\mathrm{r}\mathrm{i}\mathrm{o}\mathrm{r}}),$

where $𝑪{}_{\mathrm{o}\mathrm{b}\mathrm{s}}$ encodes observation uncertainty and $(\overline{𝒒}{}_{\mathrm{p}\mathrm{r}\mathrm{i}\mathrm{o}\mathrm{r}},𝑪{}_{\mathrm{p}\mathrm{r}\mathrm{i}\mathrm{o}\mathrm{r}})$ specify prior biochemical plausibility on the transformed scale. Equivalently, the observation model can be written in regression form as

$$𝒙{}^{\ast }=𝑸{}_{\mathrm{o}\mathrm{b}\mathrm{s}}𝒒+\mathit{ϵ},\mathit{ϵ}\sim 𝒩(\mathrm{𝟎},𝑪{}_{\mathrm{o}\mathrm{b}\mathrm{s}}),$$

so that parameter balancing is a Bayesian linear regression problem in the latent basis vector $𝒒$, with design matrix $𝑸{}_{\mathrm{o}\mathrm{b}\mathrm{s}}$, response vector $𝒙{}^{\ast }$, and Gaussian prior on the regression coefficients. This makes explicit the probabilistic structure that is only implicit in classical parameter balancing.

Under this linear-Gaussian model, the posterior is again Gaussian and is available in closed form:

	$𝑪{}_{\mathrm{p}\mathrm{o}\mathrm{s}\mathrm{t}}$	$=(𝑪{}_{\mathrm{p}\mathrm{r}\mathrm{i}\mathrm{o}\mathrm{r}}{}^{-1}+𝑸{}_{\mathrm{o}\mathrm{b}\mathrm{s}}{}^{\top }𝑪{}_{\mathrm{o}\mathrm{b}\mathrm{s}}{}^{-1}𝑸{}_{\mathrm{o}\mathrm{b}\mathrm{s}}){}^{-1},$
	$\overline{𝒒}{}_{\mathrm{p}\mathrm{o}\mathrm{s}\mathrm{t}}$	$=𝑪{}_{\mathrm{p}\mathrm{o}\mathrm{s}\mathrm{t}}(𝑸{}_{\mathrm{o}\mathrm{b}\mathrm{s}}{}^{\top }𝑪{}_{\mathrm{o}\mathrm{b}\mathrm{s}}{}^{-1}𝒙{}^{\ast }+𝑪{}_{\mathrm{p}\mathrm{r}\mathrm{i}\mathrm{o}\mathrm{r}}{}^{-1}\overline{𝒒}{}_{\mathrm{p}\mathrm{r}\mathrm{i}\mathrm{o}\mathrm{r}}).$

The posterior mean $\overline{𝒒}{}_{\mathrm{p}\mathrm{o}\mathrm{s}\mathrm{t}}$ gives the balanced basis parameter estimates, while $𝑪{}_{\mathrm{p}\mathrm{o}\mathrm{s}\mathrm{t}}$ quantifies their uncertainty. Posterior summaries for the full parameter vector follow by propagation through $𝑸$, giving

$\overline{𝒙}{}_{\mathrm{p}\mathrm{o}\mathrm{s}\mathrm{t}}$

$=𝑸\overline{𝒒}{}_{\mathrm{p}\mathrm{o}\mathrm{s}\mathrm{t}},𝑪{}_{𝒙,\mathrm{p}\mathrm{o}\mathrm{s}\mathrm{t}}=𝑸𝑪{}_{\mathrm{p}\mathrm{o}\mathrm{s}\mathrm{t}}𝑸{}^{\top }.$

Thus, the balanced parameter set is obtained together with uncertainty for both basis and derived quantities.

Parameters without direct observations remain prior-informed but constrained through $𝑸$, while conflicting measurements are reconciled according to their uncertainties. The posterior predictive distribution also provides a basis for held-out calibration checks, consistent with the broader Bayesian workflow perspective on model checking and evaluation (Gabry et al., 2019; Gelman et al., 2026).

To address misspecification, outliers, and source heterogeneity, we extended the Gaussian baseline with Student-$t$ errors for heavy tails, skew-normal errors for asymmetric noise, and source-level random effects for between-source variation. The robust likelihoods use conditional Gaussian representations, while the hierarchical term preserves the thermodynamic structure encoded by $𝑸$.

For the raw-data analyses, let $𝒙{}^{\ast }=(x{}_1{}^{\ast },\mathrm{\dots },x{}_N{}^{\ast }){}^{\top }\in ℝ{}^N$ denote the log-transformed records, and let $𝑸{}^{\ast }\in ℝ{}^{N\times D}$ be the raw-observation design matrix whose $n$-th row maps record $n$ to the corresponding derived quantity. We write

$$x{}_n{}^{\ast }=(𝑸{}^{\ast }𝒒){}_n+\eta {}_{g(n)}+ϵ{}_n,$$

where $g(n)$ denotes the publication or curated dataset from which record $n$ was extracted. The source-specific random intercept $\eta {}_{g(n)}$ captures broad publication- or dataset-level shifts without introducing a separate random effect for every source–parameter combination. We assume

$$\eta {}_g\sim 𝒩(0,\tau {}^2),g=1,\mathrm{\dots },G,$$

independently across sources. The Student-$t$ model uses a degrees-of-freedom parameter $\nu >2$, while the skew-normal model uses shape parameter $\alpha$. In raw BRENDA analyses, $\nu$ was fixed at a weakly heavy-tailed value and $\alpha$ was profiled within each training fold; in synthetic experiments, the robust fits used the corresponding data-generating or sweep values described in Supplementary Sections D.3 and F.1. The thermodynamic block is included in every fold through the same Gaussian contribution used in the baseline model.

We evaluated predictive calibration on raw BRENDA records using leakage-free fold-wise cross-validation. Let $c(n)$ denote the fold key for record $n$, which may correspond to the record index, the source label $g(n)$, the Enzyme Commission classification of the associated enzyme, or the derived-quantity index, depending on the experiment. The mapping $n\mapsto c(n)$ partitions the observation indices into disjoint folds, and each fold is evaluated by holding out one partition at a time and fitting the model on its complement. For fold $c$, define the held-out and training sets by

$$ℐ{}_c=\{n\in \{1,\mathrm{\dots },N\}:c(n)=c\},ℐ{}_{-c}=\{1,\mathrm{\dots },N\}\setminus ℐ{}_c.$$

All fitting, hyperparameter learning, and robustness-parameter handling used only the training records indexed by $ℐ{}_{-c}$, together with the always-included thermodynamic block.

For each held-out record $n\in ℐ{}_c$, let $g{}^{\ast }=g(n)$ denote its source. If $g{}^{\ast }$ appears in the training fold, prediction conditions on the sampled source effect $\eta {}_{g{}^{\ast }}{}^{(s)}$. If $g{}^{\ast }$ is absent from training, as in leave-one-source-out cross-validation, we instead draw

$$\eta {}_{\mathrm{n}\mathrm{e}\mathrm{w}}{}^{(s)}\sim 𝒩(0,\widehat{\tau }{}^2)$$

independently across posterior draws, where $\widehat{\tau }{}^2$ is the fitted between-source variance. This is the appropriate hierarchical predictive rule for unseen sources.

Given a posterior draw $𝒒{}^{(s)}$ and the corresponding source effect, the predictive mean for held-out record $n$ is

$$\mu {}_n{}^{(s)}=(𝑸{}^{\ast }𝒒{}^{(s)}){}_n+\eta {}_n{}^{(s)},$$

where $\eta {}_n{}^{(s)}$ denotes either $\eta {}_{g{}^{\ast }}{}^{(s)}$ or $\eta {}_{\mathrm{n}\mathrm{e}\mathrm{w}}{}^{(s)}$, depending on whether the source is represented in training. We then generated predictive draws from the same observation family used in the fitted model, using the fitted robustness parameters $\widehat{\nu }$ or $\widehat{\alpha }$ where applicable, together with the record-specific variance $\sigma {}_n{}^2$. For the Gaussian model with an unseen source, the source effect can be integrated out analytically, yielding

$$x{}_n{}^{\ast }\mid 𝒒{}^{(s)}\sim 𝒩((𝑸{}^{\ast }𝒒{}^{(s)}){}_n,\sigma {}_n{}^2+\widehat{\tau }{}^2)$$

in the non-contaminated case. For Student-$t$ and skew-normal models, we estimated the predictive law by Monte Carlo over the source effect and the likelihood-specific latent variables.

From the $S$ posterior predictive draws $\{x{}_n{}^{\ast ,(s)}\}{}_{s=1}{}^S$, we compute central predictive intervals, probability integral transform (PIT) values, continuous ranked probability score (CRPS), weighted interval score (WIS), and log predictive density (LPD); see Supplementary Sections E.4, E.6 and E.5. For nominal level $p$, if $[L{}_{n,p},U{}_{n,p}]$ is the central $100p\%$ predictive interval, then

$$C{}_{n,p}=𝕀\{L{}_{n,p}\le x{}_n{}^{\ast }\le U{}_{n,p}\}.$$

Empirical coverage averages $C{}_{n,p}$ over held-out records.

We evaluated inference under controlled misspecification by generating synthetic datasets from the reference pathway parameter sets. Let $𝒒{}_{\mathrm{t}\mathrm{r}\mathrm{u}\mathrm{e}}$ denote the reference basis vector. Synthetic observations were generated from

$$𝒙{}^{\ast (r)}=𝑸{}_{\mathrm{o}\mathrm{b}\mathrm{s}}𝒒{}_{\mathrm{t}\mathrm{r}\mathrm{u}\mathrm{e}}+\mathit{ϵ}{}^{(r)},r=1,\mathrm{\dots },R,$$

under three data-generating families: Gaussian, Student-$t$, and skew-normal noise. This allowed us to compare well-specified settings with heavy-tailed and asymmetric misspecification. For each simulated dataset, we fitted the Gaussian baseline and, where appropriate, the corresponding robust alternative. The main synthetic recovery experiments use the variance-matched robust families defined in Supplementary Section D.3. The controlled misspecification sweeps are deliberately more diagnostic: the degrees-of-freedom sweep varies the canonical Student-$t$ tail parameter while matching the robust comparator to the corresponding sweep variance, and the skewness sweep uses the location-parameterized skew-normal family with the robust comparator fixed at the corresponding sweep value.

We evaluated both parameter recovery and predictive calibration. For synthetic experiments, we summarized point estimation by bias and root mean squared error, and uncertainty calibration by empirical coverage of nominal credible intervals across simulation replicates. For held-out prediction, we used leave-one-out and leave-one-source-out cross-validation with posterior-predictive intervals to assess whether uncertainty was well calibrated on unseen observations. In the main text, we focus on the primary calibration and predictive findings; full interval definitions are provided in Supplementary Sections E.2 and E.4, predictive scoring rules in Supplementary Section E.5, additional coverage summaries in Supplementary Section F.1, and computational scaling results in Supplementary Section F.2.

For the Gaussian baseline, posterior means and covariances were computed analytically. Robust and hierarchical fits used conditional Gaussian representations with Gibbs-type MCMC updates for latent variables and source effects, embedded in the fold-wise MCEM procedure Wei and Tanner (1990); Levine and Casella (2001); see Supplementary Sections D.3, D.4, D.6 and D.7.

We first examined interval calibration under controlled synthetic misspecification. When data were generated under the Gaussian model, the classical Gaussian formulation achieved close-to-nominal coverage, as expected. However, when the data-generating mechanism was heavy-tailed or skewed, the Gaussian model systematically under-covered. This pattern was observed in both glycolysis and the TCA cycle: credible intervals from the classical model became too narrow as tails thickened or skew increased.

This pattern is shown clearly in Figure 3. For glycolysis, the Gaussian intervals remain close to the nominal 95% target when the data are nearly Gaussian, but coverage drops steadily as the tails thicken (i.e., as the degrees of freedom decrease) in panel (a) or as skew increases ($\alpha$ diverges from zero) in panel (b). The same qualitative behaviour appears for the TCA cycle in Figure 4. The decline is not subtle: under strong misspecification, the Gaussian baseline can fall well below the nominal level, indicating systematic overconfidence rather than mere Monte Carlo fluctuation.

By contrast, robust likelihoods substantially improved calibration. Student-$t$ likelihoods restored coverage under heavy-tailed noise, and skew-normal likelihoods improved performance under asymmetric noise. The main effect was not a dramatic change in point estimates, but a correction of posterior uncertainty so that the intervals behaved much more nearly as advertised.

Across both pathways, bias and root mean squared error (RMSE) were broadly similar across Gaussian and robust fits in well-specified and mildly misspecified settings. Thus, classical and robust parameter balancing often produced comparable point estimates even when their uncertainty quantification differed sharply. This separation between point accuracy and interval calibration is important for interpretation. A method can remain competitive in RMSE while still providing misleadingly narrow intervals.

The complete numerical summaries for the 200-replicate misspecification sweeps are reported in Supplementary Section F.1, together with the coverage curves in Figures 3 and 4 and the point-accuracy diagnostics in Supplementary Figure S2. These summaries confirm that robust fits can substantially improve calibration under heavy-tailed or asymmetric misspecification while preserving broadly similar point accuracy in well-specified and mildly misspecified regimes.

We next turn to held-out predictive assessment based on leave-one-record-out and leave-one-source-out analyses. Detailed numerical summaries for these experiments, including source-level 95% predictive coverage and predictive RMSE ranges, are reported in Supplementary Section F.1; the main text focuses on the diagnostic patterns in the source-level panels.

For glycolysis, the first source-level result is shown in Figure 5. The fixed-effects models under-cover across most nominal predictive levels, especially in the upper tail, indicating that the fitted predictive intervals are again too narrow. Once source-level random effects are introduced, all three observation families move markedly closer to the ideal diagonal. This is an important empirical finding, because it shows that a substantial share of the miscalibration arises not merely from heavy tails or skewness, but from systematic variation across laboratories, assays, or data sources. In this setting, accounting explicitly for source heterogeneity yields a larger gain in predictive reliability than changing the likelihood family alone.

Figure 6 makes the same point from a grouped perspective. Coverage at the nominal 95% level is shown by derived-quantity group together with exact binomial acceptance bands. Random-effects formulations improve calibration most clearly for the small and medium groups, where source-level variability is more visible relative to within-group sample size. This grouped view supports the interpretation from Figure 5: the calibration gain is systematic and not confined to a single subset of the pathway.

Full diagnostic panels for glycolysis and the TCA cycle are reported in Supplementary Figures S9 and S10. These panels combine multilevel predictive coverage, interval-width heatmaps, CRPS, WIS, and nominal 95% coverage summaries. They show that average scoring rules change only modestly across fitted models, whereas calibration improves substantially once source-level random effects are included. The interval-width heatmaps further indicate that the added uncertainty is concentrated in a restricted subset of kinetic quantities rather than spread uniformly across the full parameter vector.

Taken together, Figures 5 and 6 and Supplementary Figures S9 and S10 support the same substantive conclusion as the synthetic experiments. Gaussian parameter balancing can remain superficially competitive in average error or RMSE while still failing in predictive calibration. By contrast, robust likelihoods and especially hierarchical source effects deliver the more reliable uncertainty quantification needed for downstream kinetic modelling.

To keep the main text focused, we defer the full derivations of the robust and hierarchical models to Supplementary Section D; detailed definitions of interval coverage, predictive coverage, PIT, and predictive scoring rules to Supplementary Sections E.2, E.4, E.6 and E.5; and extended empirical results, including combined held-out diagnostic panels, additional group-wise and multilevel summaries, computational scaling benchmarks, and further pathway-specific sensitivity analyses, to Supplementary Sections F and F.2.

Overall, the empirical evidence is highly consistent across both pathways and across both simulation-based and held-out evaluations. Classical Gaussian parameter balancing works well when its assumptions are approximately correct, but can become overconfident under heavy tails, skewness, or source heterogeneity. Robust likelihoods and hierarchical random effects preserve the thermodynamic consistency of parameter balancing while delivering substantially better-calibrated uncertainty.

Kinetic models of cellular metabolism require many quantities, including Michaelis constants, turnover rates, elementary kinetic constants, equilibrium constants, and related thermodynamic terms. Literature and database values from resources such as BRENDA, SABIO-RK, and eQuilibrator are valuable but often sparse, noisy, heterogeneous, and mutually inconsistent (Schomburg et al., 2004; Wittig et al., 2011; Flamholz et al., 2012; Noor et al., 2014; Beber et al., 2022). Parameter balancing addresses this problem by working with independent basis quantities $𝒒$ and a dependency matrix $𝑸$, so that the complete derived vector $𝒙=𝑸𝒒$ satisfies the thermodynamic and algebraic constraints implied by the pathway. This preserves the key strength of classical parameter balancing (Lubitz et al., 2010; Lubitz and Liebermeister, 2019): missing and conflicting measurements can be reconciled without abandoning physical consistency.

The Bayesian formulation makes this reconciliation explicit. Observed log-kinetic quantities are linked to $𝒒$ through $𝑸{}_{\mathrm{o}\mathrm{b}\mathrm{s}}𝒒$, while priors $(\overline{𝒒}{}_{\mathrm{p}\mathrm{r}\mathrm{i}\mathrm{o}\mathrm{r}},𝑪{}_{\mathrm{p}\mathrm{r}\mathrm{i}\mathrm{o}\mathrm{r}})$ encode biochemical plausibility and dependency structure. In the Gaussian case, the resulting posterior is the conjugate linear update underlying classical parameter balancing, but written as a fully specified generative model. This view clarifies what uncertainty means: unobserved quantities remain constrained by $𝑸$, conflicting measurements are weighted by their uncertainty, and posterior or posterior-predictive intervals can be checked against held-out data. The detailed kinetic motivation and notation are given in Sections B and C.

The robust and hierarchical extensions address two common departures from the Gaussian fixed-effects idealization. Student-$t$ and skew-normal likelihoods target record-level errors that are heavy-tailed or asymmetric, respectively. Source-level random effects target systematic differences across publications, laboratories, assay protocols, organisms, temperatures, pH values, and extraction procedures. These extensions modify the observation model but keep the same thermodynamic dependency map $𝑸$, so the balanced parameters remain constrained by the biochemical structure.

The fold-wise predictive construction is important for interpreting the raw BRENDA analyses. All fitting and hyperparameter updates are performed using only the training records and the always-included thermodynamic block. Held-out records do not enter the training-fold updates. If a held-out record comes from a source seen during training, prediction can condition on the fitted source effect; if the whole source is held out, prediction integrates over a new source effect from the fitted between-source distribution. The mathematical details of these likelihoods, source effects, and predictive rules are provided in Section D.

The main figures should be read primarily as evidence about uncertainty calibration rather than as evidence that the central balanced parameter values change dramatically. In the synthetic pathway studies, all model specifications recover the broad thermodynamic and kinetic structure, and most posterior means remain similar across likelihood families. The main differences appear in interval width and coverage: Gaussian intervals behave well under near-Gaussian noise but under-cover when the data are heavy-tailed or skewed, whereas the Student-$t$ and skew-normal models improve calibration under their corresponding departures. In the raw BRENDA held-out analyses, source-level random effects provide the strongest calibration gain, especially for leave-one-source-out prediction, indicating that between-source heterogeneity is a major component of predictive uncertainty at pathway scale.

The practical message is therefore that parameter balancing should be judged not only by plausible point estimates, but also by calibrated uncertainty. The present hierarchical model uses one shared intercept per source, which is a parsimonious first-order representation of broad source-level shifts. Richer source-by-quantity structures or mechanistic condition corrections may be useful when different assay conditions affect parameter classes in different ways. The supplementary sections below give the formal diagnostics and empirical detail needed to assess these choices: evaluation metrics in Section E, computational and predictive results in Section F, and scalable exact-ELBO structured variational fitting in Section D.9.

The symbols in Table S1 are used primarily in the supplementary derivations for the raw-data observation model, variance-matched robust likelihoods, latent-variable augmentations, and hierarchical source effects; see especially Sections D.2, D.3, D.5 and D.6.

The symbols in Table S2 are used in the supplementary definitions of simulation metrics, held-out predictive assessment, and calibration diagnostics; see Sections E.2, E.3, E.4, E.5, E.6 and E.7.

Consider a metabolic reaction network, for example a pathway such as glycolysis or the tricarboxylic acid (TCA) cycle. Under standard well-mixed assumptions, the network dynamics can be written as a system of ordinary differential equations (ODEs)

$$\frac{d𝑿}{dt}=𝑺𝒗(𝑿,𝜽),$$

where $𝑿$ is the vector of metabolite concentrations, $𝑺$ is the stoichiometric matrix, $𝒗(𝑿,𝜽)$ is the vector of reaction rates, and $𝜽$ denotes the collection of kinetic parameters governing those rates. In biochemical kinetics, these parameters typically include catalytic constants, Michaelis constants, enzyme abundances, metabolite concentrations, and thermodynamic quantities such as equilibrium constants or standard chemical potentials.

Although the present work does not fit the ODE dynamics directly, this mechanistic setting determines which biochemical quantities must be estimated and which constraints they must satisfy. In particular, a realistic parameter set must be compatible not only with the observed data, but also with the pathway stoichiometry and the laws of thermodynamics. Parameter balancing addresses this problem by inferring a complete, thermodynamically consistent parameter vector from incomplete, noisy, and heterogeneous measurements.

To illustrate the central ideas, consider a reversible reaction of the form

$$\text{Citrate}\rightleftharpoons \text{Aconitate}+\text{H}{}_2\text{O}.$$

A simple mass-action description is

$$v=k{}^{+}[\text{Citrate}]-k{}^{-}[\text{Aconitate}][\text{H}{}_2\text{O}],$$

where $k{}^{+}$ and $k{}^{-}$ are the forward and reverse rate constants. This form is useful for intuition, but enzyme-catalyzed reactions are more commonly represented by saturation-type laws such as Michaelis–Menten or related reversible kinetics.

Following Liebermeister et al. (2010), we adopt convenience kinetics as a flexible reversible rate law that remains biochemically plausible while allowing a systematic thermodynamic parameterization. For the reaction above, one such form is

$$v=e{}_0\frac{k{}_{\mathrm{c}\mathrm{a}\mathrm{t},f}{\displaystyle \frac{[\text{Citrate}]}{K{}_{M,\text{Citrate}}}}-k{}_{\mathrm{c}\mathrm{a}\mathrm{t},r}{\displaystyle \frac{[\text{Aconitate}]}{K{}_{M,\text{Aconitate}}}}{\displaystyle \frac{[\text{H}{}_2\text{O}]}{K{}_{M,\text{H}{}_2\text{O}}}}}{1+{\displaystyle \frac{[\text{Citrate}]}{K{}_{M,\text{Citrate}}}}+{\displaystyle \frac{[\text{Aconitate}]}{K{}_{M,\text{Aconitate}}}}+{\displaystyle \frac{[\text{H}{}_2\text{O}]}{K{}_{M,\text{H}{}_2\text{O}}}}+{\displaystyle \frac{[\text{Aconitate}][\text{H}{}_2\text{O}]}{K{}_{M,\text{Aconitate}}K{}_{M,\text{H}{}_2\text{O}}}}},$$

where $e{}_0$ is the enzyme concentration, $k{}_{\mathrm{c}\mathrm{a}\mathrm{t},f}$ and $k{}_{\mathrm{c}\mathrm{a}\mathrm{t},r}$ are the forward and reverse catalytic constants, and the $K{}_{M,\cdot }$ terms are Michaelis constants. This example already shows the main challenge: a realistic kinetic description typically involves several interacting parameters, many of which are only partially observed in practice.

At equilibrium, the net flux must vanish. For reversible enzyme kinetics, this requirement induces algebraic constraints linking catalytic constants, Michaelis constants, and the reaction equilibrium constant. For the example above, the associated Haldane relationship is

$$\frac{k{}_{\mathrm{c}\mathrm{a}\mathrm{t},f}}{k{}_{\mathrm{c}\mathrm{a}\mathrm{t},r}}\cdot \frac{K{}_{M,\text{H}{}_2\text{O}}K{}_{M,\text{Aconitate}}}{K{}_{M,\text{Citrate}}}=K{}_{\mathrm{e}\mathrm{q}},$$

(2)

where $K{}_{\mathrm{e}\mathrm{q}}$ is the equilibrium constant, which may be derived from thermodynamic quantities such as standard Gibbs free energies (Flamholz et al., 2012).

Equation Equation 2 shows that the kinetic quantities in Section C.2 are not independent. If one were to draw $k{}_{\mathrm{c}\mathrm{a}\mathrm{t}}$ and $K{}_M$ values independently from literature ranges, the resulting parameter set would generally fail to satisfy Equation 2, and hence would be thermodynamically inconsistent. In small examples this may already lead to implausible kinetics; in larger pathways and genome-scale networks, where many such constraints must hold simultaneously, naive combination of literature values becomes even less reliable.

This is the main motivation for parameter balancing: the task is not merely to smooth noisy measurements, but to infer a complete parameter set that is statistically regularized and physically admissible.

Parameter balancing resolves these dependencies by introducing a set of independent basis quantities $𝒒$, from which all derived quantities $𝒙$ are computed. For the reversible reaction above, one may introduce a reparameterization in terms of independent transformed quantities such as a reaction-specific rate parameter and species-specific affinity or binding terms. For example, one may write

	$k{}_{\mathrm{c}\mathrm{a}\mathrm{t},f}$	$=\overline{r}\rho {}_{\text{Citrate}},$		(3)
	$k{}_{\mathrm{c}\mathrm{a}\mathrm{t},r}$	$=\overline{r}\rho {}_{\text{Aconitate}}\rho {}_{\text{H}{}_2\text{O}},$		(4)
	$K{}_{M,\text{Citrate}}$	$={\displaystyle \frac{\rho {}_{\text{Citrate}}}{K{}_{\text{Citrate}}}},$		(5)
	$K{}_{M,\text{Aconitate}}$	$={\displaystyle \frac{\rho {}_{\text{Aconitate}}}{K{}_{\text{Aconitate}}}},$		(6)
	$K{}_{M,\text{H}{}_2\text{O}}$	$={\displaystyle \frac{\rho {}_{\text{H}{}_2\text{O}}}{K{}_{\text{H}{}_2\text{O}}}},$		(7)

where $\overline{r}$ is a reaction-specific rate parameter and the $\rho$ and $K$ terms encode transformed affinity or binding information. Substituting Equation 3 into Equation 2 yields

$$\frac{k{}_{\mathrm{c}\mathrm{a}\mathrm{t},f}}{k{}_{\mathrm{c}\mathrm{a}\mathrm{t},r}}\cdot \frac{K{}_{M,\text{H}{}_2\text{O}}K{}_{M,\text{Aconitate}}}{K{}_{M,\text{Citrate}}}=\frac{K{}_{\text{Citrate}}}{K{}_{\text{Aconitate}}K{}_{\text{H}{}_2\text{O}}},$$

so that the equilibrium constraint is automatically satisfied when

$$K{}_{\mathrm{e}\mathrm{q}}=\frac{K{}_{\text{Citrate}}}{K{}_{\text{Aconitate}}K{}_{\text{H}{}_2\text{O}}}.$$

In the bond-graph and modular-kinetics literature (Gawthrop and Crampin, 2014; Pan et al., 2021), $\overline{r}$ is often interpreted as a transformed rate parameter, whereas the $\rho$ terms capture species-dependent energetic contributions and the $K$ terms encode effective binding information. The exact interpretation depends on the chosen kinetic formalism, but the common idea is that the original kinetic parameters are replaced by a smaller set of independent quantities from which they can be reconstructed in a thermodynamically consistent way.

The exact form of the reparameterization depends on the kinetic formalism and network representation; see Liebermeister et al. (2010); Lubitz et al. (2010) for details. For the present work, the key point is that after appropriate transformation, the resulting dependencies become linear. We therefore write

$$𝒙=𝑸𝒒,$$

where $𝒒$ collects the basis quantities, $𝒙$ collects the full set of derived kinetic and thermodynamic quantities, and $𝑸$ is a constant dependency matrix determined by the pathway structure and the chosen thermodynamic parameterization.

This linear form is what makes parameter balancing computationally convenient. For example, one row of $𝑸$ may encode that $\log k{}_{\mathrm{c}\mathrm{a}\mathrm{t},f}$ is the sum of $\log \overline{r}$ and $\log \rho {}_{\text{Citrate}}$, whereas another row may encode a different combination for $\log k{}_{\mathrm{c}\mathrm{a}\mathrm{t},r}$ or $\log K{}_{M,\text{Citrate}}$. More generally, whenever multiplicative kinetic relationships become additive on the transformed scale, they can be represented by suitable rows of $𝑸$.

The main-text Bayesian formulation builds directly on Section C.4. Observed or curated quantities are represented by the data vector $𝒙{}^{\ast }$, and the observation map $𝑸{}_{\mathrm{o}\mathrm{b}\mathrm{s}}$ selects the relevant rows of $𝑸$. The Gaussian baseline model is

	$𝒙{}^{\ast }\mid 𝒒$	$\sim 𝒩(𝑸{}_{\mathrm{o}\mathrm{b}\mathrm{s}}𝒒,𝑪{}_{\mathrm{o}\mathrm{b}\mathrm{s}}),$
	$𝒒$	$\sim 𝒩(\overline{𝒒}{}_{\mathrm{p}\mathrm{r}\mathrm{i}\mathrm{o}\mathrm{r}},𝑪{}_{\mathrm{p}\mathrm{r}\mathrm{i}\mathrm{o}\mathrm{r}}),$

so that parameter balancing becomes Bayesian inference on the basis quantities $𝒒$. The posterior distribution combines the biochemical constraints encoded by $𝑸$, the uncertainty structure encoded by $𝑪{}_{\mathrm{o}\mathrm{b}\mathrm{s}}$, and prior biochemical plausibility encoded by $(\overline{𝒒}{}_{\mathrm{p}\mathrm{r}\mathrm{i}\mathrm{o}\mathrm{r}},𝑪{}_{\mathrm{p}\mathrm{r}\mathrm{i}\mathrm{o}\mathrm{r}})$.

Posterior summaries are then propagated back to the derived scale through $𝑸$. In this way, the thermodynamic and kinetic consistency of the full parameter set is inherited automatically from the basis representation, rather than being imposed after the fact.

The discussion above illustrates the central rationale of parameter balancing. The goal is not simply to interpolate missing values, but to infer a complete parameter vector that is numerically plausible, statistically regularized, and biochemically admissible. The dependency matrix $𝑸$ provides the bridge between these requirements:

1. 1.

   it encodes pathway-specific biochemical structure,

2. 2.

   it enforces thermodynamic consistency through the choice of basis quantities, and

3. 3.

   it enables statistical inference through a linear observation model.

This structure is what makes both the classical Gaussian model and the robust Bayesian extensions in the main text computationally tractable. The robust likelihoods modify the observation model, and the hierarchical extensions add source-specific variation, but the biochemical dependency structure itself remains encoded through the same linear map $𝑸$. As a result, the framework can incorporate richer statistical modelling without losing the thermodynamic coherence that motivates parameter balancing in the first place.

The Gaussian parameter balancing model is attractive because it yields closed-form posterior updates and calibrated intervals when the observation model is correctly specified. In practice, however, curated biochemical datasets often contain outliers, heavy-tailed deviations, asymmetric measurement errors, and source-specific biases. Under such misspecification, the Gaussian model can become overconfident: posterior means may remain reasonable, but the corresponding credible intervals may under-cover because the assumed likelihood understates the probability of large or asymmetric deviations.

This issue is especially relevant in parameter balancing, where data are often assembled from multiple laboratories, assay protocols, and reporting standards. In the linear-Gaussian setting, posterior intervals coincide with classical confidence intervals and therefore have the expected coverage under correct specification. This property need not survive under misspecification, which is why empirical coverage and posterior-predictive calibration are central diagnostics in our study. The robust likelihoods considered here are designed to address these effects directly. The Student-$t$ model protects against rare but extreme observations through heavier tails, whereas the skew-normal model accommodates directional asymmetry in the residual distribution.

For the raw kinetic-data analysis, let

$$𝒙{}^{\ast }=(x{}_1{}^{\ast },\mathrm{\dots },x{}_N{}^{\ast }){}^{\top }\in ℝ{}^N$$

denote the vector of log-transformed records, and let $𝑸{}^{\ast }\in ℝ{}^{N\times D}$ be the raw-observation design matrix whose $n$-th row maps record $n$ to the corresponding derived quantity. Each record has reported log-scale variance $\sigma {}_n{}^2>0$, and we write

$$𝑪{}_{\mathrm{o}\mathrm{b}\mathrm{s}}=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}(\sigma {}_1{}^2,\mathrm{\dots },\sigma {}_N{}^2).$$

In addition to the raw BRENDA records, we retain the thermodynamic block

$$𝒚{}_0=𝑸{}_0𝒒+\mathit{ϵ}{}_0,\mathit{ϵ}{}_0\sim 𝒩(\mathrm{𝟎},𝚺{}_0),$$

which is present in every fold. Its information contribution is

$$𝑨{}_0=𝑸{}_0{}^{\top }𝚺{}_0{}^{-1}𝑸{}_0,𝒃{}_0=𝑸{}_0{}^{\top }𝚺{}_0{}^{-1}𝒚{}_0.$$

To model between-source heterogeneity, let $g(n)\in \{1,\mathrm{\dots },G\}$ denote the source label of record $n$, and introduce source-specific random intercepts

$$\eta {}_g\sim 𝒩(0,\tau {}^2),g=1,\mathrm{\dots },G,$$

independently. The latent mean for record $n$ is then

$$\mu {}_n=(𝑸{}^{\ast }𝒒){}_n+\eta {}_{g(n)}.$$

Hence $\tau {}^2$ quantifies systematic between-source variation beyond the reported within-record measurement variances. This source-effect formulation is the same object used in the raw-data cross-validation protocol, where held-out sources in leave-one-source-out cross-validation are predicted by integrating over a new intercept $\eta {}_{\mathrm{n}\mathrm{e}\mathrm{w}}\sim 𝒩(0,\tau {}^2)$ rather than by conditioning on any training-fold estimate.

We consider three observation families for the within-record residual $\mathit{ϵ}{}_n=x{}_n{}^{\ast }-\mu {}_n,$ each parameterized so that the variance matches $\sigma {}_n{}^2$. This ensures that differences in performance are not artifacts of trivial rescaling.

We evaluated predictive calibration on raw BRENDA records using leakage-free fold-wise cross-validation. Let $c(n)$ denote a fold key. Depending on the experiment, $c(n)$ may be the record index, the source label $g(n)$, the Enzyme Commission classification of the associated enzyme, or the derived-quantity index. For a fold level $c$, define

$$ℐ{}_c=\{n\in \{1,\mathrm{\dots },N\}:c(n)=c\},ℐ{}_{-c}=\{1,\mathrm{\dots },N\}\setminus ℐ{}_c.$$

All model fitting, hyperparameter learning, and robustness-parameter handling are performed using only the training records indexed by $ℐ{}_{-c}$, together with the always-included thermodynamic block. No statistic computed from $\{x{}_n{}^{\ast }:n\in ℐ{}_c\}$ enters any training-fold update.

For Student-$t$ and skew-normal models, the robustness-shape parameter is handled within the training fold and never estimated from held-out records. In the raw BRENDA analyses, the Student-$t$ degrees of freedom were fixed at $\nu =50$, giving a weakly heavy-tailed comparator to the Gaussian model. In the synthetic Student-$t$ recovery experiments, the fitted Student-$t$ model used the prespecified value $\nu =5$, matching the data-generating design. In the controlled degrees-of-freedom sweep, the Student-$t$ comparator used the current sweep value. For skew-normal analyses, $\alpha$ was initialized at the symmetric value and then updated by a one-dimensional profile step within the training fold, except in the controlled skewness sweep, where the comparator used the current sweep value to isolate the effect of asymmetric misspecification. Conditional on the current or fixed robustness-shape value, the fold-specific empirical-Bayes hyperparameters are estimated by MCEM.

Viewed at the level of exact objective updates, this block-wise scheme is closest to an expectation conditional-maximization (ECM) or expectation/conditional maximization either (ECME) construction: the location and variance hyperparameters $(𝒒{}_{\mathrm{p}\mathrm{r}\mathrm{i}\mathrm{o}\mathrm{r}},\lambda ,\tau {}^2)$ are updated by conditional maximisation of an augmented criterion, whereas an estimated robustness-shape parameter is handled by a one-dimensional observed-likelihood or profile-likelihood step. Exact ECM and ECME updates inherit the monotone-ascent property of expectation–maximization (EM) when each conditional step is solved against the appropriate complete- or observed-data objective (Meng and Rubin, 1993; Liu and Rubin, 1994; Meng and van Dyk, 1997); analogous shape updates are standard in maximum-likelihood fitting of Student-$t$ models (Liu and Rubin, 1995). Because our expectation-step (E-step) quantities are approximated by Monte Carlo, this monotonicity should be interpreted as an exact-algorithm template rather than as a finite-$S$ guarantee for the MCEM procedure.

The robust likelihoods admit latent-variable augmentations under which the complete-data likelihood is Gaussian in $(𝒒,𝜼)$. Let $\mathrm{\ell }$ denote the collection of latent variables introduced by the chosen family. Conditional on $\mathrm{\ell }$, the training-fold likelihood can be written in the unified quadratic form

$$p(𝒙{}_{-c}{}^{\ast }\mid 𝒒,𝜼,\mathrm{\ell })\propto \exp \{-\frac{1}{2}\sum {}_{n\in ℐ{}_{-c}}w{}_n(\mathrm{\ell })(\tilde{x}{}_n{}^{\ast }(\mathrm{\ell })-(𝑸{}^{\ast }𝒒){}_n-\eta {}_{g(n)}\left){}^2\right\},$$

for some effective responses $\tilde{x}{}_n{}^{\ast }(\mathrm{\ell })$ and positive weights $w{}_n(\mathrm{\ell })$.

For the Gaussian model,

$$w{}_n(\mathrm{\ell })=\sigma {}_n{}^{-2},\tilde{x}{}_n{}^{\ast }(\mathrm{\ell })=x{}_n{}^{\ast }.$$

For the Student-$t$ model,

$$w{}_n(\mathrm{\ell })=\frac{w{}_n}{s{}_n{}^2},\tilde{x}{}_n{}^{\ast }(\mathrm{\ell })=x{}_n{}^{\ast }.$$

For the skew-normal model,

$$w{}_n(\mathrm{\ell })=\frac{1}{\omega {}_n{}^{2}c{}_{\alpha }{}^2},\tilde{x}{}_n{}^{\ast }(\mathrm{\ell })=x{}_n{}^{\ast }-\xi {}_n-\omega {}_n\delta {}_{\alpha }U{}_n.$$

Let $𝒁{}_{-c}\in \{0,1\}{}^{N{}_{\mathrm{t}\mathrm{r}}\times G{}_{\mathrm{t}\mathrm{r}}}$ denote the source-incidence matrix on the training fold, so that $(𝒁{}_{-c}𝜼){}_n=\eta {}_{g(n)}$ for $n\in ℐ{}_{-c}$. Let $𝑾=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}(w{}_n:n\in ℐ{}_{-c})$. Combining Section D.5, the thermodynamic block, and the priors

$$𝒒\mid 𝒒{}_{\mathrm{p}\mathrm{r}\mathrm{i}\mathrm{o}\mathrm{r}},\lambda \sim 𝒩(𝒒{}_{\mathrm{p}\mathrm{r}\mathrm{i}\mathrm{o}\mathrm{r}},\lambda 𝑪{}_0),\eta {}_g\mid \tau {}^2\sim 𝒩(0,\tau {}^2),$$

gives the conditional posterior

$$(𝒒,𝜼)\mid 𝒙{}_{-c}{}^{\ast },\mathrm{\ell },𝒒{}_{\mathrm{p}\mathrm{r}\mathrm{i}\mathrm{o}\mathrm{r}},\lambda ,\tau {}^2\sim 𝒩(𝒎,𝑷{}^{-1}),$$

with block precision

$$𝑷=\left(\begin{array}{cc}𝑷{}_{𝒒𝒒}& 𝑷{}_{𝒒\eta }\\ 𝑷{}_{\eta 𝒒}& 𝑫{}_{\eta \eta }\end{array}\right),$$

where

	$𝑷{}_{𝒒𝒒}$	$=\lambda {}^{-1}𝑪{}_0{}^{-1}+𝑨{}_0+(𝑸{}_{-c}{}^{\ast }){}^{\top }𝑾𝑸{}_{-c}{}^{\ast },$
	$𝑷{}_{𝒒\eta }$	$=(𝑸{}_{-c}{}^{\ast }){}^{\top }𝑾𝒁{}_{-c},$
	$𝑫{}_{\eta \eta }$	$=\tau {}^{-2}𝑰{}_{G{}_{\mathrm{t}\mathrm{r}}}+𝒁{}_{-c}{}^{\top }𝑾𝒁{}_{-c}.$

The corresponding right-hand side is

$$𝒓=\left(\begin{array}{c}\lambda {}^{-1}𝑪{}_0{}^{-1}𝒒{}_{\mathrm{p}\mathrm{r}\mathrm{i}\mathrm{o}\mathrm{r}}+𝒃{}_0+(𝑸{}_{-c}{}^{\ast }){}^{\top }𝑾\widetilde{𝒙}{}_{-c}{}^{\ast }\\ 𝒁{}_{-c}{}^{\top }𝑾\widetilde{𝒙}{}_{-c}{}^{\ast }\end{array}\right).$$

Because each training record belongs to exactly one source, $𝒁{}_{-c}{}^{\top }𝑾𝒁{}_{-c}$ is diagonal, and therefore $𝑫{}_{\eta \eta }$ is diagonal as well. This diagonal structure is the key to the collapsed Gibbs update.

The direct Gibbs update samples $(𝒒,𝜼)$ jointly from the Gaussian conditional above. We instead integrate out $𝜼$ analytically. Since $𝑫{}_{\eta \eta }$ is diagonal and positive definite, the Schur complement gives the exact marginal precision of $𝒒$:

$$𝑷{}_{𝒒}{}^{\mathrm{m}\mathrm{a}\mathrm{r}\mathrm{g}}=𝑷{}_{𝒒𝒒}-𝑷{}_{𝒒\eta }𝑫{}_{\eta \eta }{}^{-1}𝑷{}_{\eta 𝒒}.$$

The corresponding marginal right-hand side is

$$𝒓{}_{𝒒}{}^{\mathrm{m}\mathrm{a}\mathrm{r}\mathrm{g}}=𝒓{}_{𝒒}-𝑷{}_{𝒒\eta }𝑫{}_{\eta \eta }{}^{-1}𝒓{}_{\eta }.$$

Hence

$$𝒒\mid 𝒙{}_{-c}{}^{\ast },\mathrm{\ell },𝒒{}_{\mathrm{p}\mathrm{r}\mathrm{i}\mathrm{o}\mathrm{r}},\lambda ,\tau {}^2\sim 𝒩((𝑷{}_{𝒒}{}^{\mathrm{m}\mathrm{a}\mathrm{r}\mathrm{g}}){}^{-1}𝒓{}_{𝒒}{}^{\mathrm{m}\mathrm{a}\mathrm{r}\mathrm{g}},(𝑷{}_{𝒒}{}^{\mathrm{m}\mathrm{a}\mathrm{r}\mathrm{g}}){}^{-1}).$$

Conditional on $𝒒$, the source effects become independent:

$$\eta {}_g\mid 𝒒,𝒙{}_{-c}{}^{\ast },\mathrm{\ell },\tau {}^2\sim 𝒩(m{}_g(𝒒),v{}_g),v{}_g=(\tau {}^{-2}+\sum {}_{n\in ℐ{}_{-c}:g(n)=g}w{}_n){}^{-1},$$

with

$$m{}_g(𝒒)=v{}_g\sum {}_{n\in ℐ{}_{-c}:g(n)=g}w{}_n(\tilde{x}{}_n{}^{\ast }-(𝑸{}^{\ast }𝒒){}_n).$$

Sampling first from the exact marginal law of $𝒒$ and then from the exact conditional law of $𝜼\mid 𝒒$ therefore produces the same joint conditional distribution as the uncollapsed Gaussian block update.

This is a collapsed Gaussian block update in the sense of collapsed Gibbs sampling (Liu, 1994): $𝜼$ is integrated out analytically when sampling $𝒒$, and is then reconstructed from its exact conditional distribution. This differs from partially collapsed Gibbs samplers, where marginalised and conditional updates may be interleaved and their order must be chosen carefully to preserve the target distribution (van Dyk and Park, 2008). For the present Gaussian block, the marginal-then-conditional composition samples from the same joint conditional distribution as the uncollapsed update. Collapsing tractable blocks is often associated with improved operator or autocorrelation behaviour under covariance-ordering comparisons for augmentation schemes and collapsed Gibbs samplers (Liu et al., 1994; Liu, 1994), although we do not require a universal ordering for the full sampler.

Within each training fold, empirical-Bayes hyperparameters are estimated by MCEM. The Gaussian prior used by parameter balancing is

$$𝒒\mid 𝒒{}_{\mathrm{p}\mathrm{r}\mathrm{i}\mathrm{o}\mathrm{r}},\lambda \sim 𝒩(𝒒{}_{\mathrm{p}\mathrm{r}\mathrm{i}\mathrm{o}\mathrm{r}},\lambda 𝑪{}_0),$$

where $𝒒{}_{\mathrm{p}\mathrm{r}\mathrm{i}\mathrm{o}\mathrm{r}}$ is the fold-specific empirical-Bayes prior center and $\lambda >0$ is a global prior-scale multiplier. To regularize this update, we place the hyperprior

$$𝒒{}_{\mathrm{p}\mathrm{r}\mathrm{i}\mathrm{o}\mathrm{r}}\sim 𝒩(𝒒{}_0,\kappa {}_0𝑪{}_0),$$

where $𝒒{}_0\in ℝ{}^D$ is the fixed thermodynamic prior center used in the classical parameter-balancing initialization, and $\kappa {}_0>0$ controls the strength of shrinkage toward this baseline. Smaller $\kappa {}_0$ keeps $𝒒{}_{\mathrm{p}\mathrm{r}\mathrm{i}\mathrm{o}\mathrm{r}}$ close to $𝒒{}_0$, whereas $\kappa {}_0\to \mathrm{\infty }$ gives the unregularized update $𝒒{}_{\mathrm{p}\mathrm{r}\mathrm{i}\mathrm{o}\mathrm{r}}\leftarrow \overline{𝒒}$.

Given Monte Carlo draws $𝒒{}^{(1)},\mathrm{\dots },𝒒{}^{(S)}$, define

$$\overline{𝒒}=\frac{1}{S}\sum {}_{s=1}{}^S𝒒{}^{(s)}.$$

The part of the MCEM objective depending on $𝒒{}_{\mathrm{p}\mathrm{r}\mathrm{i}\mathrm{o}\mathrm{r}}$, up to additive constants, is

$$-\frac{1}{2\lambda S}\sum {}_{s=1}{}^S(𝒒{}^{(s)}-𝒒{}_{\mathrm{p}\mathrm{r}\mathrm{i}\mathrm{o}\mathrm{r}}){}^{\top }𝑪{}_0{}^{-1}(𝒒{}^{(s)}-𝒒{}_{\mathrm{p}\mathrm{r}\mathrm{i}\mathrm{o}\mathrm{r}})-\frac{1}{2\kappa {}_0}(𝒒{}_{\mathrm{p}\mathrm{r}\mathrm{i}\mathrm{o}\mathrm{r}}-𝒒{}_0){}^{\top }𝑪{}_0{}^{-1}(𝒒{}_{\mathrm{p}\mathrm{r}\mathrm{i}\mathrm{o}\mathrm{r}}-𝒒{}_0).$$

Differentiating with respect to $𝒒{}_{\mathrm{p}\mathrm{r}\mathrm{i}\mathrm{o}\mathrm{r}}$ gives

$$\lambda {}^{-1}𝑪{}_0{}^{-1}(\overline{𝒒}-𝒒{}_{\mathrm{p}\mathrm{r}\mathrm{i}\mathrm{o}\mathrm{r}})+\kappa {}_0{}^{-1}𝑪{}_0{}^{-1}(𝒒{}_0-𝒒{}_{\mathrm{p}\mathrm{r}\mathrm{i}\mathrm{o}\mathrm{r}})=0.$$

Since $𝑪{}_0$ is positive definite, this is equivalent to

$$(\lambda {}^{-1}+\kappa {}_0{}^{-1})𝒒{}_{\mathrm{p}\mathrm{r}\mathrm{i}\mathrm{o}\mathrm{r}}=\lambda {}^{-1}\overline{𝒒}+\kappa {}_0{}^{-1}𝒒{}_0.$$

Therefore, the MCEM update for the empirical-Bayes prior center is

$$𝒒{}_{\mathrm{p}\mathrm{r}\mathrm{i}\mathrm{o}\mathrm{r}}\leftarrow \frac{\lambda {}^{-1}\overline{𝒒}+\kappa {}_0{}^{-1}𝒒{}_0}{\lambda {}^{-1}+\kappa {}_0{}^{-1}},$$

and the inverse-gamma maximum a posteriori (MAP) update for $\lambda$ is

$$\lambda \leftarrow \frac{b{}_{\lambda }+\frac{1}{2}𝒟{}_{\lambda }}{a{}_{\lambda }+1+\frac{D}{2}}.$$

For the source-variance parameter $\tau {}^2$, a direct MCEM update would use

$$S{}_{\eta \eta }=\frac{1}{S}\sum {}_{s=1}{}^S\parallel 𝜼{}^{(s)}\parallel {}^2.$$

Instead, we use the Rao–Blackwellized estimator

$$S{}_{\eta \eta }{}^{\mathrm{R}\mathrm{B}}=\frac{1}{S}\sum {}_{s=1}{}^S𝔼[\parallel 𝜼\parallel {}^2\mid 𝒒{}^{(s)},𝒙{}_{-c}{}^{\ast },\mathrm{\ell }{}^{(s)}]=\frac{1}{S}\sum {}_{s=1}{}^S\sum {}_{g=1}{}^{G{}_{\mathrm{t}\mathrm{r}}}(m{}_g(𝒒{}^{(s)}){}^2+v{}_g).$$

By the law of total variance,

$$\mathrm{V}\mathrm{a}\mathrm{r}\left(𝔼\right[\parallel 𝜼\parallel {}^2\mid 𝒒,𝒙{}_{-c}{}^{\ast },\mathrm{\ell }\left]\right)\le \mathrm{V}\mathrm{a}\mathrm{r}(\parallel 𝜼\parallel {}^2),$$

so the Rao–Blackwellized estimator removes the extra Monte Carlo variability generated by sampling $𝜼$ conditionally on $𝒒$. More precisely, for a fixed target distribution, the conditional expectation has no larger ordinary variance than the raw plug-in statistic $\parallel 𝜼{}^{(s)}\parallel {}^2$, by the Rao–Blackwell inequality (Gelfand and Smith, 1990; Casella and Robert, 1996). For Markov-chain output, the relevant object is the asymptotic variance of the ergodic average; Rao–Blackwellisation can reduce this quantity under additional transition-kernel and reversibility or covariance-ordering conditions (Liu et al., 1994; Robert and Roberts, 2021), but such dominance is not automatic for arbitrary non-reversible or adaptive chains (Geyer, 1992). We therefore use $S{}_{\eta \eta }{}^{\mathrm{R}\mathrm{B}}$ as a variance-reduction device for the $\tau {}^2$ update, while treating the finite-$S$ MCEM run as a monitored Monte Carlo approximation. The resulting MAP update is

$$\tau {}^2\leftarrow \frac{b{}_{\tau }+\frac{1}{2}S{}_{\eta \eta }{}^{\mathrm{R}\mathrm{B}}}{a{}_{\tau }+1+\frac{G{}_{\mathrm{t}\mathrm{r}}}{2}}.$$

If contamination modelling is enabled, letting

$$Z=\frac{1}{S}\sum {}_{s=1}{}^S\sum {}_{n\in ℐ{}_{-c}}z{}_n{}^{(s)},$$

the posterior-mean update for $\pi$ is

$$\pi \leftarrow \frac{a{}_{\pi }+Z}{a{}_{\pi }+b{}_{\pi }+N{}_{\mathrm{t}\mathrm{r}}}.$$