Maximum likelihood estimation of log-affine models using detailed-balanced reaction networks

TL;DR

This work tackles maximum likelihood estimation for log-affine models by realizing the MLE as the unique positive steady state of a mass-action network G_Lambda,c, constructed from a finite spanning set Lambda of ker(A) that spans ker(A). It generalizes Gopalkrishnan's approach to arbitrary spanning sets and proves that using a Markov basis eliminates relevant boundary steady states, guaranteeing global convergence to the MLE within the positive class. By leveraging Birch's theorem and toric geometry, the paper clarifies how the spanning set shapes deficiency and convergence properties, with explicit analyses on the independence model and Hardy–Weinberg example illustrating fast, input-independent convergence when Lambda approaches a Markov basis. The results offer a principled, dynamical pathway to compute MLEs in log-affine models via detailed-balanced networks, while emphasizing a practical trade-off between network size and convergence robustness for molecular implementations.

Abstract

A fundamental question in the field of molecular computation is what computational tasks a biochemical system can carry out. In this work, we focus on the problem of finding the maximum likelihood estimate (MLE) for log-affine models. We revisit a construction due to Gopalkrishnan of a mass-action system with the MLE as its unique positive steady state, which is based on choosing a basis for the kernel of the design matrix of the model. We extend this construction to allow for any finite spanning set of the kernel, and explore how the choice of spanning set influences the dynamics of the resulting network, including the existence of boundary steady states, the deficiency of the network, and the rate of convergence. In particular, we prove that using a Markov basis as the spanning set guarantees global stability of the MLE steady state.

Figures (5)

• Figure 1: Illustration of Birch's Theorem for the Hardy--Weinberg model (\ref{['ex:hardy_weinberg']}). The model $\mathcal{M}_{\mathbf{A}}$ (orange) is the intersection of $\mathcal{X}_{\mathbf{A},{\boldsymbol{c}}}$ (blue) and the probability simplex $\Delta_2$ (gray). Birch's Theorem guarantees that for any observed distribution $\bar{{\boldsymbol{u}}}\in\Delta_{2}$, there is a unique point $\widehat{{\boldsymbol{p}}}$ in the intersection between $\overline{\mathcal{X}}_{\mathbf{A}}$ and $\bar{{\boldsymbol{u}}} + \ker(\mathbf{A})$ (green), and that this is the MLE.
• Figure 2: Illustration of \ref{['ex-DYN-stoichiometric_compatibility_class-a']}, showing the set of positive steady states $\mathcal{E}_{G,{\boldsymbol{k}}}$ (blue), an initial state ${\boldsymbol{x}}_0$ and the corresponding stoichiometric compatibility class ${\boldsymbol{x}}_0+S$ (green), and the unique positive steady state $\widehat{{\boldsymbol{x}}}$ (orange), and the phase portrait restricted to $\Delta_2$ (light green). Note that $\dim(S) = 1$.
• Figure 3: Illustration of \ref{['ex:boundary_steady_states']}, showing trajectories and steady states for several initial conditions under the canonical projection to $(x_1,x_2,x_3)$. Red boxes: The observed distribution $\bar{{\boldsymbol{u}}} = (0.9,0,0,0.1)^{}$ is a boundary steady state for system given by $\Lambda_1$, while the system given by $\Lambda_\mathrm{mb}$ evolves from $\bar{{\boldsymbol{u}}}$ to the MLE.
• Figure 4: The least negative nonzero eigenvalues for systems constructed using increasing sequences of spanning sets $\Lambda_i$, which increase towards either (solid) the unique Markov basis $\Lambda_{10} = \Lambda_{\mathrm{mb}}$, or (dashed) a random spanning set. The models considered were (a) the binary 4-cycle model (\ref{['ex:4cycle']}), and (b) the independence model (\ref{['ex:independence']}) with $r_1 = 10$ and $r_2 = 10$.
• Figure 5: All nonzero eigenvalues of the systems constructed using $\Lambda_i$, which increases towards the Markov basis $\Lambda_{10}$ for (a) the binary 4-cycle model (\ref{['ex:4cycle']}), and (b) the independence model (\ref{['ex:independence']}) with $r_1 = 10$ and $r_2 = 10$. As $\Lambda_i$ increases towards the Markov basis, all nonzero eigenvalues tend to become more negative.

Theorems & Definitions (46)

• Definition 2.1
• Definition 2.2
• Lemma 2.3
• Example 2.4
• Theorem 2.5: Birch's Theorem
• Example 2.6
• Example 2.7
• Example 2.8
• Definition 3.1
• Definition 3.2