A reaction network scheme which implements inference and learning for Hidden Markov Models

TL;DR

This paper introduces the Chemical Baum-Welch Algorithm, a novel reaction network scheme that enables molecular systems to learn Hidden Markov Model parameters through biochemical reactions, mimicking traditional algorithms.

Contribution

It presents a biochemical reaction network implementation of the Baum-Welch algorithm for HMMs, with proven convergence properties and simulation validation.

Findings

Reaction network converges exponentially fast in expectation and maximization steps.

Fixed points of the network correspond to those of the Baum-Welch algorithm.

Simulations show the network learns parameters identical to classical HMM algorithms.

Abstract

With a view towards molecular communication systems and molecular multi-agent systems, we propose the Chemical Baum-Welch Algorithm, a novel reaction network scheme that learns parameters for Hidden Markov Models (HMMs). Each reaction in our scheme changes only one molecule of one species to one molecule of another. The reverse change is also accessible but via a different set of enzymes, in a design reminiscent of futile cycles in biochemical pathways. We show that every fixed point of the Baum-Welch algorithm for HMMs is a fixed point of our reaction network scheme, and every positive fixed point of our scheme is a fixed point of the Baum-Welch algorithm. We prove that the "Expectation" step and the "Maximization" step of our reaction network separately converge exponentially fast. We simulate mass-action kinetics for our network on an example sequence, and show that it learns the same parameters for the HMM as the Baum-Welch algorithm.

Tables

Table 1. Table 1: Baum-Welch Reaction Network : The steps and reactions are for all g , h ∈ H , w ∈ V formulae-sequence 𝑔 ℎ 𝐻 𝑤 𝑉 g,h\in H,w\in V and l = 1 , … , L − 1 𝑙 1 … 𝐿 1 l=1,\dots,L-1 . Notice there are some null reactions of the form ∑ a i ​ X i → ∑ a i ​ X i → subscript 𝑎 𝑖 subscript 𝑋 𝑖 subscript 𝑎 𝑖 subscript 𝑋 𝑖 \sum a_{i}X_{i}\to\sum a_{i}X_{i} . As these null reaction have no effect on the dynamics, they can be ignored.

Baum-Welch Algorithm	Baum-Welch Reaction Network
${\alpha }_{1h}={\pi }_h{\psi }_{h{v}_1}$	${\alpha }_{1h}+{\pi }_{h^{\ast }}+{\psi }_{h^{\ast }w}+E_{1w}$ $\text{ce}->\left[\right]{\alpha }_{1h^{\ast }}+{\pi }_{h^{\ast }}+{\psi }_{h^{\ast }w}+E_{1w}$ ${\alpha }_{1h^{\ast }}+{\pi }_h+{\psi }_{hw}+E_{1w}$ $\text{ce}->\left[\right]{\alpha }_{1h}+{\pi }_h+{\psi }_{hw}+E_{1w}$
${\alpha }_{l+1,h}=\sum _{g\in H}{\alpha }_{lg}{\theta }_{gh}{\psi }_{h{v}_{l+1}}$	$\begin{array}{c}{\alpha }_{l+1,h}+{\alpha }_{lg}+{\theta }_{g{h}^{\ast }}+{\psi }_{h^{\ast }w}+E_{l+1,w}\text{ce}->\left[\right]\\ {\alpha }_{l+1,h^{\ast }}+{\alpha }_{lg}+{\theta }_{g{h}^{\ast }}+{\psi }_{h^{\ast }w}+E_{l+1,w}\\ {\alpha }_{l+1,h^{\ast }}+{\alpha }_{lg}+{\theta }_{gh}+{\psi }_{hw}+E_{l+1,w}\text{ce}->\left[\right]\\ {\alpha }_{l+1,h}+{\alpha }_{lg}+{\theta }_{gh}+{\psi }_{hw}+E_{l+1,w}\end{array}$
${\beta }_{Lh}=1$	$\begin{array}{c}{\beta }_{lh}+{\beta }_{l+1,g}+{\theta }_{h^{\ast }g}+{\psi }_{gw}+E_{l+1,w}\text{ce}->\left[\right]\\ {\beta }_{l{h}^{\ast }}+{\beta }_{l+1,g}+{\theta }_{h^{\ast }g}+{\psi }_{gw}+E_{l+1,w}\\ {\beta }_{l{h}^{\ast }}+{\beta }_{l+1,g}+{\theta }_{hg}+{\psi }_{gw}+E_{l+1,w}\text{ce}->\left[\right]\\ {\beta }_{lh}+{\beta }_{l+1,g}+{\theta }_{hg}+{\psi }_{gw}+E_{l+1,w}\end{array}$
${\beta }_{lh}=\sum _{g\in H}{\theta }_{hg}{\psi }_{g{v}_{l+1}}{\beta }_{l+1,g}$
${\gamma }_{lh}=\frac{{\alpha }_{lh}{\beta }_{lh}}{{\sum }_{g\in H}{\alpha }_{lg}{\beta }_{lg}}$	$\begin{array}{c}{\gamma }_{lh}+{\alpha }_{l{h}^{\ast }}+{\beta }_{l{h}^{\ast }}\text{ce}->\left[\right]{\gamma }_{l{h}^{\ast }}+{\alpha }_{l{h}^{\ast }}+{\beta }_{l{h}^{\ast }}\\ {\gamma }_{l{h}^{\ast }}+{\alpha }_{lh}+{\beta }_{lh}\text{ce}->\left[\right]{\gamma }_{lh}+{\alpha }_{lh}+{\beta }_{lh}\end{array}$
${\xi }_{lgh}=\frac{{\alpha }_{lg}{\theta }_{gh}{\psi }_{h{v}_{l+1}}{\beta }_{l+1,h}}{{\sum }_{f\in H}{\alpha }_{lf}{\beta }_{lf}}$
${\theta }_{gh}\leftarrow \frac{{\sum }_{l=1}^{L-1}{\xi }_{lgh}}{{\sum }_{l=1}^{L-1}{\sum }_{f\in H}{\xi }_{lgf}}$	$\begin{array}{c}{\theta }_{gh}+{\xi }_{lg{h}^{\ast }}\text{ce}->\left[\right]{\theta }_{g{h}^{\ast }}+{\xi }_{lg{h}^{\ast }}\\ {\theta }_{g{h}^{\ast }}+{\xi }_{lgh}\text{ce}->\left[\right]{\theta }_{gh}+{\xi }_{lgh}\end{array}$
${\psi }_{hw}\leftarrow \frac{{\sum }_{l=1}^L{\gamma }_{lh}{\delta }_{w,v_l}}{{\sum }_{l=1}^L{\gamma }_{lh}}$	$\begin{array}{c}{\psi }_{hw}+{\gamma }_{lh}+E_{l{v}^{\ast }}\text{ce}->\left[\right]{\psi }_{h{v}^{\ast }}+{\gamma }_{lh}+E_{l{v}^{\ast }}\\ {\psi }_{h{v}^{\ast }}+{\gamma }_{lh}+E_{lw}\text{ce}->\left[\right]{\psi }_{hw}+{\gamma }_{lh}+E_{lw}\end{array}$

Equations

$${\gamma }_{lh}=\frac{{\alpha }_{lh}{\beta }_{lh}}{{\sum }_{g\in H}{\alpha }_{lg}{\beta }_{lg}}$$

$${\xi }_{lgh}=\frac{{\alpha }_{lg}{\theta }_{gh}{\psi }_{h{v}_{l+1}}{\beta }_{l+1,h}}{{\sum }_{f\in H}{\alpha }_{lf}{\beta }_{lf}}$$

$${\theta }_{gh}\leftarrow \frac{{\sum }_{l=1}^{L-1}{\xi }_{lgh}}{{\sum }_{l=1}^{L-1}{\sum }_{f\in H}{\xi }_{lgf}},$$

$$\sum _{X\in S}{y}_{X}X\to \sum _{X\in S}{y}_X^{\prime }X$$

$${\dot{x}}_i(t)=\sum _{a\to b\in R}{k}_{a\to b}(b_i-a_i)\prod _{j\in S}{x}_j(t{)}^{a_j}$$

$$\begin{array}{c}{\alpha }_{1h}+{\pi }_{h^{\ast }}+{\psi }_{h^{\ast }{v}_1}\text{ce}->[]{\alpha }_{1h^{\ast }}+{\pi }_{h^{\ast }}+{\psi }_{h^{\ast }{v}_1}\\ {\alpha }_{1h^{\ast }}+{\pi }_h+{\psi }_{h{v}_1}\text{ce}->[]{\alpha }_{1h}+{\pi }_h+{\psi }_{h{v}_1}\end{array}$$

$$\begin{array}{c}{\alpha }_{lh}+{\alpha }_{l-1,g}+{\theta }_{g{h}^{\ast }}+{\psi }_{h^{\ast }{v}_l}\text{ce}->[]{\alpha }_{l{h}^{\ast }}+{\alpha }_{l-1,g}+{\theta }_{g{h}^{\ast }}+{\psi }_{h^{\ast }{v}_l}\\ {\alpha }_{l{h}^{\ast }}+{\alpha }_{l-1,g}+{\theta }_{gh}+{\psi }_{h{v}_l}\text{ce}->[]{\alpha }_{lh}+{\alpha }_{l-1,g}+{\theta }_{gh}+{\psi }_{h{v}_l}\end{array}$$

$${\alpha }_{lh}\times \left(\sum _{g\in H}{\alpha }_{l-1,g}{\theta }_{g{h}^{\ast }}{\psi }_{h^{\ast }{v}_l}\right)={\alpha }_{l{h}^{\ast }}\times \left(\sum _{g\in H}{\alpha }_{l-1,g}{\theta }_{gh}{\psi }_{h{v}_l}\right)$$

$$\begin{array}{c}{\alpha }_{1h}+{\pi }_{h^{\ast }}+{\psi }_{h^{\ast }w}+E_{1w}\text{ce}->[]{\alpha }_{1h^{\ast }}+{\pi }_{h^{\ast }}+{\psi }_{h^{\ast }w}+E_{1w}\\ {\alpha }_{1h^{\ast }}+{\pi }_h+{\psi }_{hw}+E_{1w}\text{ce}->[]{\alpha }_{1h}+{\pi }_h+{\psi }_{hw}+E_{1w}\end{array}$$

$$\begin{array}{c}{\alpha }_{lh}+{\alpha }_{l-1,g}+{\theta }_{g{h}^{\ast }}+{\psi }_{h^{\ast }w}+E_{lw}\text{ce}->[]{\alpha }_{l{h}^{\ast }}+{\alpha }_{l-1,g}+{\theta }_{g{h}^{\ast }}+{\psi }_{h^{\ast }w}+E_{lw}\\ {\alpha }_{l{h}^{\ast }}+{\alpha }_{l-1,g}+{\theta }_{gh}+{\psi }_{hw}+E_{lw}\text{ce}->[]{\alpha }_{lh}+{\alpha }_{l-1,g}+{\theta }_{gh}+{\psi }_{hw}+E_{lw}\end{array}$$

$$\begin{array}[]{|c||c|c|c|c|c|c|c|c|}\hline\cr\text{Type}&\pi&\alpha&\beta&\theta&\psi&E&\gamma&\xi\\ \hline\cr\text{Number}&\,|H|&\,|H|L&\,|H|L&\,|H|^{2}&\,|H||V|&\,L|V|&\,|H|L&\,|H|^{2}(L-1)\\ \hline\cr\end{array}$$

$$\begin{array}[]{|c|c|c|c|}\hline\cr\text{Forward}&\text{Backward}&\text{Expectation}&\text{Maximization}\\ \hline\cr 2(|H|-1)V+&(2|H|(|H|-1)|V|)&2L(|H|-1)+&2|H|(|H|-1)(L-1)\\ 2|H|(|H|-1)(L-1)|V|&(L-1)&2(L-1)(|H|^{2}-1)&+2|H|L(|V|-1)\\ \hline\cr\end{array}$$

$$\dot{\alpha}_{1h}\big{|}_{\Phi}=k^{\alpha}_{1h}\left({\alpha^{\prime}_{1h^{*}}}{\pi_{h}\psi^{\prime}_{hv_{1}}}-{\alpha^{\prime}_{1h}}{\pi_{h^{*}}\psi^{\prime}_{h^{*}v_{1}}}\right)=0$$

$$\dot{\alpha}_{lh}\big{|}_{\Phi}=k^{\alpha}_{lh}\left({\alpha^{\prime}_{lh^{*}}}{\sum_{g\in H}\alpha^{\prime}_{l-1,g}\theta^{\prime}_{gh}\psi^{\prime}_{hv_{l}}}-{\alpha^{\prime}_{lh}}{\sum_{g\in H}\alpha^{\prime}_{l-1,g}\theta^{\prime}_{gh^{*}}\psi^{\prime}_{h^{*}v_{l}}}\right)=0$$

$$\dot{\beta}_{lh}\big{|}_{\Phi}=k^{\beta}_{lh}\left({\beta^{\prime}_{lh^{*}}}{\sum_{g\in H}\theta^{\prime}_{h^{*}g}\psi^{\prime}_{gv_{l+1}}\beta^{\prime}_{l+1,g}}-{\beta^{\prime}_{lh}}{\sum_{g\in H}\theta^{\prime}_{hg}\psi^{\prime}_{gv_{l+1}}\beta^{\prime}_{l+1,g}}\right)=0$$

$${\dot{\gamma }}_{lh}=k_{lh}^{\gamma }\left({\gamma }_{l{h}^{\ast }}^{\prime }{\alpha }_{lh}^{\prime }{\beta }_{lh}^{\prime }-{\gamma }_{lh}^{\prime }{\alpha }_{l{h}^{\ast }}^{\prime }{\beta }_{l{h}^{\ast }}^{\prime }\right)=0$$

$$\dot{\xi}_{lgh}\big{|}_{\Phi}=k^{\xi}_{lgh}\left({\xi^{\prime}_{lh^{*}h^{*}}}{\alpha^{\prime}_{lg}\theta^{\prime}_{gh}\psi^{\prime}_{hv_{l+1}}\beta^{\prime}_{l+1,h}}-{\xi^{\prime}_{lgh}}{\alpha^{\prime}_{lh^{*}}\theta^{\prime}_{h^{*}h^{*}}\psi^{\prime}_{h^{*}v_{l+1}}\beta^{\prime}_{l+1,h^{*}}}\right)=0$$

$$\dot{\theta}_{gh}\big{|}_{\Phi}=k^{\theta}_{gh}\left({\theta^{\prime}_{gh^{*}}}{\sum_{l=1}^{L-1}\xi^{\prime}_{lgh}}-{\theta^{\prime}_{gh}}{\sum_{l=1}^{L-1}\xi^{\prime}_{l}(g,h^{*})}\right)=0$$

$${\alpha }_{1h}^{\prime }=\left({\pi }_h{\psi }_{h{v}_1}^{\prime }\right)\frac{{\sum }_{f\in H}{\alpha }_{1f}^{\prime }}{{\sum }_{f\in H}{\pi }_f{\psi }_{f{v}_1}^{\prime }}\text{ for all }h\in H$$

$${\alpha }_{lh}^{\prime }=\left(\sum _{g\in H}{\alpha }_{l-1,g}^{\prime }{\theta }_{gh}^{\prime }{\psi }_{h{v}_l}^{\prime }\right)\frac{{\sum }_{f\in H}{\alpha }_{lf}^{\prime }}{{\sum }_{f,g\in H}{\alpha }_{l-1,g}^{\prime }{\theta }_{gf}^{\prime }{\psi }_{f{v}_l}^{\prime }}\text{ for all }h\in H\text{ and }l=2,\dots ,L$$

$${\beta }_{lh}^{\prime }=\left(\sum _{g\in H}{\theta }_{hg}^{\prime }{\psi }_{g{v}_{l+1}}^{\prime }{\beta }_{l+1,g}^{\prime }\right)\frac{{\sum }_{f\in H}{\beta }_{lf}^{\prime }}{{\sum }_{f,g\in H}{\theta }_{fg}^{\prime }{\psi }_{g{v}_{l+1}}^{\prime }{\beta }_{l+1,g}^{\prime }}\text{ for all }h\in H\text{ and }l=1,\dots ,L-1$$

$${\gamma }_{lh}^{\prime }=\left(\frac{{\alpha }_{lh}^{\prime }{\beta }_{lh}^{\prime }}{{\sum }_{g\in H}{\alpha }_{lg}^{\prime }{\beta }_{lg}^{\prime }}\right)\sum _{g\in H}{\gamma }_{lg}^{\prime }\text{ for all }h\in H\text{ and }l=1,2,\dots ,L-1$$

$${\xi }_{lgh}^{\prime }=\left(\frac{{\alpha }_{lg}^{\prime }{\theta }_{gh}^{\prime }{\psi }_{h{v}_{l+1}}^{\prime }{\beta }_{l+1,h}^{\prime }}{{\sum }_{e,f\in H}{\alpha }_{lf}^{\prime }{\theta }_{ef}^{\prime }{\psi }_{f{v}_{l+1}}^{\prime }{\beta }_{l+1,f}^{\prime }}\right)\sum _{e,f\in H}{\xi }_{lef}^{\prime }\text{ for all }g,h\in H\times H\text{ and }l=1,\dots ,L-1$$

$${\theta }_{gh}^{\prime }=\left(\frac{{\sum }_{l=1}^{L-1}{\xi }_{lgh}^{\prime }}{{\sum }_{f\in H}{\sum }_{l=1}^{L-1}{\xi }_{lgf}^{\prime }}\right)\sum _{f\in H}{\theta }_{gf}^{\prime }\text{ for all }g\in H\text{ and }h\in H$$

$${\psi }_{hw}^{\prime }=\left(\frac{{\sum }_{l=1}^L{\gamma }_{lh}^{\prime }{\delta }_{w,v_l}}{{\sum }_{l=1}^L{\gamma }_{lh}^{\prime }}\right)\sum _{v\in V}{\psi }_{hv}^{\prime }\text{ for all }h\in H\text{ and }w\in V$$

$$W=\left(\begin{array}{cc}{W}^{\prime }& I_r\end{array}\right)$$

$$A=\left(\begin{array}{cc}{A}_{11}& A_{12}\\ A_{21}& A_{22}\end{array}\right),$$

$$P=\left(\begin{array}{cc}{I}_{n-r}& 0_{n-r,r}\\ W^{\prime }& I_r\end{array}\right),P^{-1}=\left(\begin{array}{cc}{I}_{n-r}& 0_{n-r,r}\\ -W^{\prime }& I_r\end{array}\right),$$

$$PA{P}^{-1}$$

Full text

1 Introduction

The sophisticated behavior of living cells on short timescales is powered by biochemical reaction networks. One may say that evolution has composed the symphony of the biosphere, genetic machinery conducts the music, and reaction networks are the orchestra. Understanding the capabilities and limits of this molecular orchestra is key to understanding living systems, as well as to engineering molecular systems that are capable of sophisticated life-like behavior.

The technology of implementing abstract reaction networks with molecules is a subfield of molecular systems engineering that has witnessed rapid advances in recent times. Several researchers [1, 2, 3, 4, 5, 6, 7, 8] have proposed theoretical schemes for implementing arbitrary reaction networks with DNA oligonucleotides. There is a growing body of experimental demonstrations of such schemes [9, 10, 2, 7, 11]. A stack of tools is emerging to help automate the design process. We can now compile abstract reaction networks to a set of DNA oligonucleotides that will implement the dynamics of the network in solution [12]. We can computationally simulate the dynamics of these oligonucleotide molecular systems [13] to allow debugging prior to experimental implementation. In view of these rapid advances, the study of reaction networks from the point of view of their computational capabilities has become even more urgent.

It has long been known that reaction networks can compute any computable function [14]. The literature has several examples of reaction network schemes that have been inspired by known algorithms [15, 16, 17, 18, 19, 20, 21, 22]. Our group has previously described reaction network schemes that solve statistical problems like maximum likelihood [23], sampling a conditional distribution and inference [24], and learning from partial observations [25]. These schemes exploit the thermodynamic nature of the underlying molecular systems that will implement these reaction networks, and can be expressed in terms of variational ideas involving minimization of Helmholtz free energy [26, 27, 28].

In this paper, we consider situations where partial information about the environment is available to a cell in the form of a sequence of observations. For example, this might happen when an enzyme is acting processively on a polymer, or a molecular walker [29, 30, 31] is trying to locate its position on a grid. In situations like this, multiple observations are not independent. Such sequences can not be summarized merely by the type of the sequence [32], i.e., the number of times various symbols occur. Instead, the order of various observations carries information about state changes in the process producing the sequence. The number of sequences grows exponentially with length, and our previously proposed schemes are algorithmically inadequate. To deal with such situations requires a pithy representation of sequences, and a way of doing inference and learning directly on such representations. In Statistics and Machine Learning, this problem is solved by Hidden Markov Models (HMMs) [33].

HMMs are a widely used model in Machine Learning, powering sequence analysis applications like speech recognition [34], handwriting recognition, and bioinformatics. They are also essential components of communication systems as well as of intelligent agents trained by reinforcement learning methods. In this article, we describe a reaction network scheme which implements the Baum-Welch algorithm. The Baum-Welch algorithm is an iterative algorithm for learning HMM parameters. Reaction networks that can do such statistical analysis on sequences are likely to be an essential component of molecular communication systems, enabling cooperative behavior among a population of artificial cells. Our main contributions are:

In Section 2, we describe what the reader needs to know about HMMs and the Baum-Welch algorithm to be able to follow the subsequent constructions. No prerequisites are assumed beyond familiarity with matrices and probability distributions. 2. 2.

In Section 3, we describe a novel reaction network scheme to learn parameters for an HMM. 3. 3.

We prove in Theorem 4.1 that every fixed point of the Baum-Welch algorithm is also a fixed point of the continuous dynamics of this reaction network scheme. 4. 4.

In Theorem 4.2, we prove that every positive fixed point of the dynamics of our reaction network scheme is a fixed point of the Baum-Welch algorithm. 5. 5.

In Theorem 4.3 and Theorem 4.4, we prove that subsets of our reaction network scheme which correspond to the Expectation step and the Maximization step of the Baum-Welch algorithm both separately converge exponentially fast. 6. 6.

In Example 1, we simulate our reaction network scheme on an input sequence and show that the network dynamics is successfully able to learn the same parameters as the Baum-Welch algorithm. 7. 7.

We show in Example 2 that when the baum-welch fixed point is on the boundary then our scheme need not converge to a Baum-Welch fixed point. However, we conjecture that if there exists a positive Baum-Welch fixed point then our scheme must always converge to a Baum-Welch fixed point. In particular, there would be a positive Baum-Welch fixed point if the true HMM generating the sequence has all parameters positive, and the observed sequence is long enough.

2 Hidden Markov Models and the Baum Welch Algorithm

Fix two finite sets $P$ and $Q$. A stochastic map is a $|P|\times|Q|$ matrix $A=(a_{pq})_{|P|\times|Q|}$ such that $a_{pq}\geq 0$ for all $p\in P$ and $q\in Q$, and $\sum_{q\in Q}a_{pq}=1$ for all $p\in P$. Intuitively, stochastic maps represent conditional probability distributions.

An HMM $(H,V,\theta,\psi,\pi)$ consists of finite sets $H$ (for ‘hidden’) and $V$ (for ‘visible’), a stochastic map $\theta$ from $H$ to $H$ called the transition matrix, a stochastic map $\psi$ from $H$ to $V$ called the emission matrix, and an initial probability distribution $\pi=(\pi_{h})_{h\in H}$ on $H$, i.e., $\pi_{h}\geq 0$ for all $h\in H$ and $\sum_{h\in H}\pi_{h}=1$. See Figure 1a for an example.

Suppose a length $L$ sequence $(v_{1},v_{2},\dots,v_{L})\in V^{L}$ of visible states is observed due to a hidden sequence $(x_{1},x_{2},\dots,x_{L})\in H^{L}$. The following questions related to an HMM are commonly studied:

Likelihood: For fixed $\theta,\psi$, compute the likelihood $\Pr(v_{1},v_{2},\dots,v_{L}\mid\theta,\psi)$. This problem is solved by the forward-backward algorithm. 2. 2.

Learning: Estimate the parameters $\hat{\theta},\hat{\psi}$ which maximize the likelihood of the observed sequence $(v_{1},v_{2},\dots,v_{L})\in V^{L}$. This problem is solved by the Baum-Welch algorithm which is an Expectation-Maximization (EM) algorithm. It uses the forward-backward algorithm as a subroutine to compute the E step of EM. 3. 3.

Decoding: For fixed $\theta,\psi$, find the sequence $(\hat{h}_{1},\hat{h}_{2},\dots,\hat{h}_{l},\dots,\hat{h}_{L})\in H^{L}$ that has the highest probability of producing the given observed sequence $(v_{1},v_{2},\dots,v_{L})$. This problem is solved by the Viterbi algorithm.

The forward algorithm (Figure 1b) takes as input an HMM $(H,V,\theta,\psi,\pi)$ and a length $L$ observation sequence $(v_{1},v_{2},\dots,v_{L})\in V^{L}$ and outputs the $L\times|H|$ likelihoods $\alpha_{lh}=\Pr[v_{1},v_{2},\dots,v_{l},x_{l}=h\mid\theta,\psi]$ of observing symbols $v_{1},\ldots,v_{l}$ and being in the hidden state $h\in H$ at time $l$. It does so using the following recursion.

• •

Initialisation: $\alpha_{1h}=\pi_{h}\psi_{hv_{1}}$ for all $h\in H$,

• •

Recursion: $\alpha_{lh}=\sum_{g\in H}\alpha_{l-1,g}\theta_{gh}\psi_{hv_{l}}$ for all $h\in H$ and $l=2,\ldots,L$.

The backward algorithm (Figure 1c) takes as input an HMM $(H,V,\theta,\psi,\pi)$ and a length $L$ observation sequence $(v_{1},v_{2},\dots,v_{L})\in V^{L}$ and outputs the $L\times|H|$ conditional probabilities 

$\beta_{lh}=\Pr[v_{l+1},v_{l+2},\dots,v_{L}\mid x_{l}=h,\theta,\psi]$ of observing symbols $v_{l+1},\ldots,v_{L}$ given that the hidden state $x_{l}$ at time $l$ has label $h\in H$.

• •

Initialisation: $\beta_{Lh}=1$,  for all $h\in H$,

• •

Recursion: $\beta_{lh}=\sum_{g\in H}\theta_{hg}\psi_{gv_{l+1}}\beta_{l+1,g}$ for all $h\in H$ and $l=1,\ldots,L-1$.

The E step for the Baum-Welch algorithm takes as input an HMM $(H,V,\theta,\psi,\pi)$ and a length $L$ observation sequence $(v_{1},v_{2},\dots,v_{L})$. It outputs the $L\times H$ conditional probabilities $\gamma_{lh}=\Pr[x_{l}=h\mid\theta,\psi,v]$ of being in hidden state $h$ at time $l$ conditioned on the observed sequence $v$ by:

[TABLE]

for all $h\in H$ and $l=1,2,\dots,L-1$. It also outputs the $(L-1)\times H\times H$ probabilities $\xi_{lgh}=\Pr[x_{l}=g,x_{l+1}=h\mid\theta,\psi,v]$ of transitioning along the edge $(g,h)$ at time $l$ conditioned on the observed sequence $v$ by:

[TABLE]

for all $g,h\in H$ and $l=1,\ldots,L-1$.

Remark 1

Note that the E-step uses the forward and backward algorithms as subroutines to first compute the $\alpha$ and $\beta$ values. Further note that we don’t need the forward and backward algorithms to return the actual values $\alpha_{lh}$ and $\beta_{lh}$. To be precise, let $\alpha_{l}=(\alpha_{lh})_{h\in H}\in\mathbb{R}^{H}$ denote the vector of forward likelihoods at time $l$. Then for the E step to work, we only need the direction of $\alpha_{l}$ and not the magnitude. This is because the numerator and denominator in the updates are both linear in $\alpha_{lh}$, and the magnitude cancels out. Similarly, if $\beta_{l}=(\beta_{lh})_{h\in H}\in\mathbb{R}^{H}$ denotes the vector of backward likelihoods at time $l$ then the $E$ step only cares about the direction of $\beta_{l}$ and not the magnitude. This scale symmetry is a useful property for numerical solvers. We will also make use of this freedom when we implement a lax forward-backward algorithm using reaction networks in the next section.

The M step of the Baum-Welch algorithm takes as input the values $\gamma$ and $\xi$ that are output by the E step as reconstruction of the dynamics on hidden states, and outputs new Maximum Likelihood estimates of the parameters $\theta,\psi$ that best explain these values. The update rule turns out to be very simple. For all $g,h\in H$ and $w\in V$:

[TABLE]

where $\delta_{w,v_{l}}=\begin{cases}1\text{ if }w=v_{l}\ \\ 0\text{ otherwise}\end{cases}$ is the Dirac delta function.

Remark 2

Like in Remark 1, note that the M step does not require its inputs to be the actual values $\gamma$ and $\xi$. There is a scaling symmetry so that we only need the directions of the vectors $\gamma(h)=(\gamma_{lh})_{l=1,2,\dots,L}\in\mathbb{R}^{L}$ for all $h\in H$ and $\xi(g)=(\xi_{lgh})_{1\leq l\leq L-1,h\in H}\in\mathbb{R}^{(L-1)\times H}$ for all $g\in H$. This gives us the freedom to implement a lax E projection without affecting the M projection, and we will exploit this freedom when designing our reaction network.

The Baum-Welch algorithm (Figure 1d) is a fixed point EM computation that alternately runs the E step and the M step till the updates become small enough. It is guaranteed to converge to a fixed point $(\hat{\theta},\hat{\psi})$. However, the fixed point need not always be a global optimum.

3 Chemical Baum-Welch Algorithm

3.1 Reaction Networks

Following [25], we recall some concepts from reaction network theory [35, 36, 37, 38, 39, 24].

Fix a finite set $S$ of species. An $S$-reaction, or simply a reaction when $S$ is understood from context, is a formal chemical equation

[TABLE]

where the numbers $y_{X},y^{\prime}_{X}\in\mathbb{Z}_{\geq 0}$ are the stoichiometric coefficients of species $X$ on the reactant side and product side respectively. A reaction network is a pair $(S,R)$ where $R$ is a finite set of $S$-reactions. A reaction system is a triple $(S,R,k)$ where $(S,R)$ is a reaction network and $k\mathrel{\mathop{\mathchar 58\relax}}R\to\mathbb{R}_{>0}$ is called the rate function.

As is common when specifying reaction networks, we will find it convenient to explicitly specify only a set of chemical equations, leaving the set of species to be inferred by the reader.

Fix a reaction system $(S,R,k)$. Deterministic Mass Action Kinetics describes a system of ordinary differential equations on the concentration variables $\{x_{i}(t)\mid i\in S\}$ according to:

[TABLE]

3.2 Baum-Welch Reaction Network

In this section we will describe a reaction networks for each part of the Baum-Welch algorithm.

Fix an HMM $\mathcal{M}=(H,V,\theta,\psi,\pi)$. Pick an arbitrary hidden state $h^{*}\in H$ and an arbitrary visible state $v^{*}\in V$. Picking these states $h^{*}\in H$ and $v^{*}\in V$ is merely an artifice to break symmetry akin to selecting leaders, and our results hold independent of these choices. Also fix a length $L\in\mathbb{Z}_{>0}$ of observed sequence.

We first work out in full detail how the forward algorithm of the Baum-Welch algorithm may be translated into chemical reactions. Suppose a length $L$ sequence $(v_{1},v_{2},\dots,v_{L})\in V^{L}$ of visible states is observed. Then recall that the forward algorithm uses the following recursion:

• •

Initialisation: $\alpha_{1h}=\pi_{h}\psi_{hv_{1}}$ for all $h\in H$,

• •

Recursion: $\alpha_{lh}=\sum_{g\in H}\alpha_{l-1,g}\theta_{gh}\psi_{hv_{l}}$ for all $h\in H$ and $l=2,\ldots,L$.

Notice this implies

• •

$\alpha_{1h}\times\pi_{h^{*}}\psi_{h^{*}v_{1}}=\alpha_{1h^{*}}\times\pi_{h}\psi_{hv_{1}}$ for all $h\in H\setminus\{h^{*}\}$,

• •

$\alpha_{lh}\times\left(\sum_{g\in H}\alpha_{l-1,g}\theta_{gh^{*}}\psi_{h^{*}v_{l}}\right)=\alpha_{lh^{*}}\times\left(\sum_{g\in H}\alpha_{l-1,g}\theta_{gh}\psi_{hv_{l}}\right)$ for all $h\in H\setminus\{h^{*}\}$ and $l=2,\ldots,L$.

This prompts the use of the following reactions for the initialization step:

[TABLE]

for all $h\in H\setminus\{h^{*}\}$ and $w\in V$. By design $\alpha_{1h}\times\pi_{h^{*}}\psi_{h^{*}v_{1}}=\alpha_{1h^{*}}\times\pi_{h}\psi_{hv_{1}}$ is the balance equation for the pair of reactions corresponding to each $h\in H\setminus\{h^{*}\}$.

Similarly for the recursion step, we use the following reactions:

[TABLE]

for all $g\in H$, $h\in H\setminus\{h^{*}\},l=2,\ldots,L$ and $w\in V$. Again by design

[TABLE]

is the balance equation for the set of reactions corresponding to each $(h,l)\in H\setminus\{h^{*}\}\times\{2,\ldots,L\}$.

The above reactions depend on the observed sequence $(v_{1},v_{2},\dots,v_{L})\in V^{L}$ of visible state. And this is a problem because one would have to design different reaction networks for different observed sequences. To solve this problem we introduce the species $E_{lw}$ with $l=1,\dots,L$ and $w\in V$. Now with the $E_{lw}$ species, we use the following reactions for the forward algorithm:

[TABLE]

for all $h\in H\setminus\{h^{*}\}$ and $w\in V$.

[TABLE]

for all $g\in H$, $h\in H\setminus\{h^{*}\},l=2,\ldots,L$ and $w\in V$. The $E_{lw}$ species are to be initialized such that $E_{lw}=1$ iff $w=v_{l}$. So different observed sequences can now be processed by the same reaction network, by appropriately initializing the species $E_{lw}$.

The other parts of the Baum-Welch algorithm may be translated into chemical reactions using a similar logic. We call the resulting reaction network as the Baum-Welch Reaction Network $BW(\mathcal{M},h^{*},v^{*},L)$. It consists of four subnetworks corresponding to the four parts of the Baum-Welch algorithm, as shown in Table LABEL:tab:BWCRN.

The Baum-Welch reaction network described above has a special structure. Every reaction is a monomolecular transformation catalyzed by a set of species. The reverse transformation is also present, catalyzed by a different set of species to give the network a “futile cycle” [40] structure. In addition, each connected component in the undirected graph representation of the network has a topology with all transformations happening to and from a central species. This prompts the following definitions.

Definition 1 (Flowers, petals, gardens)

A graph is a triple $(\text{Nodes},\text{Edges},\pi)$ where Nodes and Edges are finite sets and $\pi$ is a map from Edges to unordered pairs of elements from Nodes. A flower is a graph with a special node $n^{*}$ such that for every edge $e\in\text{Edges}$ we have $n^{*}\in\pi(e)$. A garden is a graph which is a union of disjoint flowers. A petal is a set of all edges $e$ which have the same $\pi(e)$, i.e. they are incident between the same pair of nodes.

Figure 2 shows how the Baum-Welch reaction network can be represented as a garden graph, in which species are nodes and reactions are edges.

A collection of specific rates is permissible if all reactions in a petal have the same rate. However, different petals may have different rates. We will denote the specific rate for a petal by superscripting the type of the species and subscripting its indices. For example, the specific rate for reactions in the petal for $\alpha_{lh}$ would be denoted as $k^{\alpha}_{lh}$. The notation for the remaining rate constants can be read from Figure  2.

Remark 3

The “flower” topology we have employed for the Baum-Welch reaction network is only one among several possibilities. The important thing is to achieve connectivity between different nodes that ought to be connected, ensuring that the ratio of the concentrations of the species denoted by adjacent nodes takes the value as prescribed by the Baum-Welch algorithm. Therefore, many other connection topologies can be imagined, for example a ring topology where each node is connected to two other nodes while maintaining the same connected components. In fact, there are obvious disadvantages to the star topology, with a single point of failure, whereas the ring topology appears more resilient. How network topology affects basins of attraction, rates of convergence, and emergence of spurious equilibria in our algorithm is a compelling question beyond the scope of this present work.

The Baum-Welch reaction network with a permissible choice of rate constants $k$ will define a Baum-Welch reaction system $(BW(\mathcal{M},h^{*},v^{*},L),k)$ whose deterministic mass action kinetics equations will perform a continuous time version of the Baum-Welch algorithm. We call this the Chemical Baum-Welch Algorithm, and describe it in Algorithm 1.

4 Analysis and Simulations

The Baum-Welch reaction network has number of species of each type as follows:

[TABLE]

The total number of species is $|H|+3|H|L+|H|^{2}+|H||V|+L|V|+|H|^{2}(L-1)$ which is $O(L|H|^{2}+|H||V|+L|V|)$. The number of reactions (ignoring the null reactions of the form $\sum a_{i}X_{i}\to\sum a_{i}X_{i}$) in each part is:

[TABLE]

so that the total number of reactions is $O(|H|^{2}L|V|)$.

The first theorem shows that the Chemical Baum-Welch Algorithm recovers all of the Baum-Welch equilibria.

Theorem 4.1

Every fixed point of the Baum-Welch algorithm for an HMM $\mathcal{M}=(H,V,\theta,\psi,\pi)$ is a fixed point for the corresponding Chemical Baum-Welch Algorithm with permissible rates $k$.

See Appendix G.1 for proof.

We say a vector of real numbers is positive if all its components are strictly greater than [math]. The next theorem shows that positive equilibria of the Chemical Baum-Welch Algorithm are also fixed points of the Baum-Welch algorithm.

Theorem 4.2

Every positive fixed point for the Chemical Baum-Welch Algorithm on a Baum-Welch Reaction system $(BW(\mathcal{M},h^{*},v^{*},L),k)$ with permissible rate $k$ is a fixed point for the Baum-Welch algorithm for the HMM $\mathcal{M}=(H,V,\theta,\psi,\pi)$.

See Appendix G.1 for proof.

The Baum-Welch algorithm is an iterative algorithm. The E step and the M step are run iteratively. In contrast, the Chemical Baum-Welch algorithm is a generalized EM algorithm [41] where all the reactions are run at the same time in a single-pot reaction.

The next theorem shows that the $E$ step consisting of the forward network, the backward network, and the $E$ network, converges exponentially fast to the correct equilibrium if the $\theta$ and $\psi$ species are held fixed at a positive point.

Theorem 4.3

For the Baum-Welch Reaction System $(BW(\mathcal{M},h^{*},v^{*},L),k)$ with permissible rates $k$, if the concentrations of $\theta$ and $\psi$ species are held fixed at a positive point then the Forward, Backward and Expection step reaction systems on $\alpha,\beta,\gamma$ and $\psi$ species converge to equilibrium exponentially fast.

See Appendix G.2 for proof.

The next theorem shows that if the $\alpha,\beta,\gamma,\xi$ species are held fixed at a positive point, then the $M$ step consisting of reactions modifying the $\theta$ and $\psi$ species converges exponentially fast to the correct equilibrium.

Theorem 4.4

For the Baum-Welch Reaction System $(BW(\mathcal{M},h^{*},v^{*},L),k)$ with permissible rates $k$, if the concentrations of $\alpha,\beta,\gamma$ and $\xi$ species are held fixed at a positive point then the Maximization step reaction system on $\theta$ and $\psi$ converges to equilibrium exponentially fast.

See Appendix G.2 for proof.  

The following examples demonstrate the behavior of the Chemical Baum-Welch Algorithm.

Example 1

Consider an HMM $(H,V,\theta,\psi,\pi)$ with two hidden states $H=\{H_{1},H_{2}\}$ and two emitted symbols $V=\{V_{1},V_{2}\}$ where the starting probability is $\pi=(0.6,0.4)$, initial transition probability is $\theta=\begin{bmatrix}0.6&0.4\\ 0.3&0.7\end{bmatrix}$, and initial emission probability is $\psi=\begin{bmatrix}0.5&0.5\\ 0.5&0.5\end{bmatrix}$. Suppose we wish to learn $(\theta,\psi)$ for the following observed sequence:

$(V_{1},V_{1},V_{1},V_{2},V_{1},V_{1},V_{2},V_{2},V_{2},V_{1},V_{2},V_{2},V_{1},V_{1},V_{1},V_{1},V_{2},V_{2},V_{2},V_{1},V_{2},V_{1},V_{1},V_{2},V_{2})$. We initialize the corresponding reaction system by setting species $E_{l,v_{l}}=1$ and $E_{l,w}=0$ for $w\neq v_{l}$, and run the dynamics according to deterministic mass-action kinetics. Initial conditions of all species that are not mentioned are chosen to be nonzero, but otherwise at random.

For our numerical solution, we observe that the reaction network equilibrium point coincides with the Baum-Welch steady state $\hat{\theta}=\begin{bmatrix}0.5071&0.4928\\ 0.0000&1.0000\end{bmatrix}$ and $\hat{\psi}=\begin{bmatrix}1.0000&0.0000\\ 0.4854&0.5145\end{bmatrix}$ (See Figure 3).

The next example shows that the Chemical Baum-Welch algorithm can sometimes get stuck at points that are not equilibria for the Baum-Welch algorithm. This is a problem especially for very short sequences, and is probably happening because in such settings, the best model sets many parameters to zero. When many species concentrations are set to [math], many reactions get turned off, and the network gets stuck away from the desired equilibrium. We believe this will not happen if the HMM generating the observed sequence has nonzero parameters, and the observed sequence is long enough.

Example 2

Consider an HMM $(H,V,\theta,\psi,\pi)$ with two hidden states $H=\{H_{1},H_{2}\}$ and two emitted symbols $V=\{V_{1},V_{2}\}$, initial distribution $\pi=(0.6,0.4)$ is fixed, initial transition probability $\theta=\begin{bmatrix}0.6&0.4\\ 0.3&0.7\end{bmatrix}$, and initial emission probability $\psi=\begin{bmatrix}0.5&0.5\\ 0.5&0.5\end{bmatrix}$. Suppose we wish to learn $(\theta,\psi)$ for the sequence $(V_{1},V_{2},V_{1},V_{2},V_{1})$. We again simulate the corresponding reaction network and also perform the Baum-Welch algorithm for comparison.

Figure 4 shows that the reaction network equilibrium point does not coincide with the Baum-Welch steady state. Both converge to $\hat{\psi}=\begin{bmatrix}1.0000&0.0000\\ 0.0000&1.0000\end{bmatrix}$. However, the reaction network converges to $\hat{\theta}=\begin{bmatrix}0.3489&0.6510\\ 1.0000&0.0000\end{bmatrix}$ whereas the Baum-Welch algorithm converges to $\hat{\theta}=\begin{bmatrix}0.0000&1.0000\\ 1.0000&0.0000\end{bmatrix}$ which happens to be the true maximum likelihood point. Note that this does not contradict Theorem 2 because the fixed point of this system is a boundary point and not a positive fixed point.

5 Related work

In previous work, our group has shown that reaction networks can perform the Expectation Maximization (EM) algorithm for partially observed log linear statistical models [25, 42]. That algorithm also applies “out of the box” to learning HMM parameters. The problem with that algorithm is that the size of the reaction network would become exponentially large in the length of the sequence, so that even examples like Example 1 with an observation sequence of length $25$ would become impractical. In contrast, the scheme we have presented in this paper requires only a linear growth with sequence length. We have obtained the savings by exploiting the graphical structure of HMMs. This allows us to compute the likelihoods $\alpha$ in a “dynamic programming” manner, instead of having to explicitly represent each path as a separate species.

Napp and Adams [20] have shown how to compute marginals on graphical models with reaction networks. They exploit graphical structure by mimicking belief propagation. Hidden Markov Models can be viewed as a special type of graphical model where there are $2L$ random variables $X_{1},X_{2},\dots,X_{L},Y_{1},Y_{2},\dots,Y_{L}$ with the $X$ random variables taking values in $H$ and the $Y$ random variables in $V$. The $X$ random variables form a Markov chain $X_{1}\,\rule[3.44444pt]{10.0pt}{0.5pt}\,X_{2}\,\rule[3.44444pt]{10.0pt}{0.5pt}\dots\rule[3.44444pt]{10.0pt}{0.5pt}\,X_{L}$. In addition, there are $L$ edges from $X_{l}$ to $Y_{l}$ for $l=1$ to $L$ denoting observations. Specialized to HMMs, the scheme of Napp and Adams would compute the equivalent of steady state values of the $\gamma$ species, performing a version of the E step. They are able to show that true marginals are fixed points of their scheme, which is similar to our Theorem 4.1. Thus their work may be viewed as the first example of a reaction network scheme that exploits graphical structure to compute E projections. Our E step goes further by proving correctness as well as exponential convergence. Their work also raises the challenge of extending our scheme to all graphical models.

Poole et al. [43] have described Chemical Boltzmann Machines, which are reaction network schemes whose dynamics reconstructs inference in Boltzmann Machines. This inference can be viewed as a version of E projection. No scheme for learning is presented. The exact schemes presented there are exponentially large. The more realistically sized schemes are presented there without proof. In comparison, our schemes are polynomially sized, provably correct if the equilibrium is positive, and perform both inference and learning for HMMs.

Zechner et al. [11] have shown that Kalman filters can be implemented with reactions networks. Kalman filters can be thought of as a version of Hidden Markov Models with continuous hidden states [44]. It would be instructive to compare their scheme with ours, and note similarities and differences. In passing from position $l$ to position $l+1$ along the sequence, our scheme repeats the same reaction network that updates $\alpha_{l+1}$ using $\alpha_{l}$ values. It is worth examining if this can be done “in place” so that the same species can be reused, and a reaction network can be described that is not tied to the length $L$ of the sequence to be observed.

Recently Cherry et al. [10] have given a brilliant experimental demonstration of learning with DNA molecules. They have empirically demonstrated a DNA molecular system that can classify $9$ types of handwritten digits from the MNIST database. Their approach is based on the notion of “winner-takes-all” circuits due to Maass [45] which was originally a proposal for how neural networks in the brain work. Winner-take-all might also be capable of approximating HMM learning, at least in theory [46], and it is worth understanding precisely how such schemes relate to the kind of scheme we have described here. It is conceivable that our scheme could be converted to winner-take-all by getting different species in the same flower to give negative feedback to each other. This might well lead to sampling the most likely path, performing a decoding task similar to the Viterbi algorithm.

6 Discussion

We have described a one-pot one-shot reaction network implementation of the Baum-Welch algorithm. Firstly, this involves proposing a reaction system whose positive fixed points correspond to equilibria of the Baum-Welch Algorithm. Secondly, this involves establishing the conditions under which convergence of solutions to the equilibria can be accomplished. Further, from a practical perspective, it is essential to obtain an implementation that does not rely on repeated iteration of a reaction scheme but only requires to be run once.

As we observe in Remark 3, there is a whole class of reaction networks that implements the Baum-Welch algorithm. We have proposed one such network and are aware that there are other networks, potentially with more efficient dynamical and analytical properties than the one proposed here. Finding and characterizing efficient reaction networks that can do complicated statistical tasks will likely be of future concern.

We have only discussed deterministic dynamical properties of the network. However, in realistic biological contexts one might imagine that the network is implemented by relatively few molecules such that stochastic effects are significant. Consequently, the study of the Baum-Welch reaction system under stochastic mass-action kinetics is likely to be of interest.

Lastly, we have mentioned the Viterbi algorithm, but have made no attempt to describe how the maximum likelihood sequence can be recovered from our reaction network. This decoding step is likely to be of as much interest for molecular communication systems and molecular multi-agent systems as it is in more traditional domains of communications and multi-agent reinforcement learning. Because of the inherent stochasticity of reaction networks, there might even be opportunities for list decoding by sampling different paths through the hidden states that have high probability conditioned on the observations. This might give an artificial cell the ability to “imagine” different possible realities, and act assuming one of them to be the case, leading to an intrinsic stochasticity and unpredictability to the behavior.

Appendix G

Appendix G.1 Comparing points of equilibria

We will now prove Theorem 4.1. For the sake of convenience, we first recall the statement.

Theorem 4.1. Every fixed point of the Baum-Welch algorithm for an HMM $\mathcal{M}=(H,V,\theta,\psi,\pi)$ is a fixed point for the corresponding Chemical Baum-Welch Algorithm with permissible rates $k$.

Proof

Consider a point $\Phi=(\alpha^{\prime},\beta^{\prime},\gamma^{\prime},\xi^{\prime},\theta^{\prime},\psi^{\prime})$ with $\alpha^{\prime},\beta^{\prime},\gamma^{\prime}\in\mathbb{R}_{\geq 0}^{L\times H}$, $\xi^{\prime}\in\mathbb{R}_{\geq 0}^{(L-1)\times H\times H}$, $\theta^{\prime}\in\mathbb{R}_{\geq 0}^{H\times H}$ and $\psi^{\prime}\in\mathbb{R}_{\geq 0}^{H\times V}$. If $\Phi$ is a fixed point of Baum-Welch Algorithm then it must satisfy:

• •

$\alpha^{\prime}_{1h}=\pi_{h}\psi^{\prime}_{hv_{1}}$ for all $h\in H$. Then for the chemical Baum-Welch Algorithm we have

[TABLE]

for all $h\in H\setminus\{h^{*}\}$ and \dot{\alpha}_{1h^{*}}\big{|}_{\Phi}=-\sum_{h\not=h^{*}}\dot{\alpha}_{1h}\big{|}_{\Phi}=0.

• •

$\alpha^{\prime}_{lh}=\sum_{g\in H}\alpha^{\prime}_{l-1,g}\theta^{\prime}_{gh}\psi^{\prime}_{hv_{l}}$ for all $h\in H$ and $l=2,\ldots,L$ . Then for the chemical Baum-Welch Algorithm we have

[TABLE]

for all $h\in H\setminus\{h^{*}\}$ and $l=2,\ldots,L$ and \dot{\alpha}_{lh^{*}}\big{|}_{\Phi}=-\sum_{h\not=h^{*}}\dot{\alpha}_{lh}\big{|}_{\Phi}=0.

• •

$\beta^{\prime}_{lh}=\sum_{g\in H}\theta^{\prime}_{hg}\psi^{\prime}_{gv_{l+1}}\beta^{\prime}_{l+1,g}$ for all $h\in H$ and $l=1,\ldots,L-1$. Then for the chemical Baum-Welch Algorithm we have

[TABLE]

for all $h\in H\setminus\{h^{*}\}$ and $l=1,2,\dots,L-1$ and \dot{\beta}_{lh^{*}}\big{|}_{\Phi}=-\sum_{h\not=h^{*}}\dot{\beta}_{lh}\big{|}_{\Phi}=0.

• •

$\gamma^{\prime}_{l}(h)=\frac{\alpha^{\prime}_{lh}\beta^{\prime}_{lh}}{\sum_{g\in H}\alpha^{\prime}_{lg}\beta^{\prime}_{lg}}$ for all $h\in H$ and $l=1,2,\dots,L-1$. Then for the chemical Baum-Welch Algorithm we have

[TABLE]

for all $h\in H\setminus\{h^{*}\}$ and $l=1,2,\dots,L-1$ and $\dot{\gamma}_{lh^{*}}=-\sum_{h\not=h^{*}}\dot{\gamma}_{lh}=0$.

• •

$\xi^{\prime}_{l}(g,h)=\frac{\alpha^{\prime}_{lg}\theta^{\prime}_{gh}\psi^{\prime}_{hv_{l+1}}\beta^{\prime}_{l+1,h}}{\sum_{f\in H}\alpha^{\prime}_{lf}\beta^{\prime}_{lf}}\text{ for all }g,h\in H\text{ and }l=1,\ldots,L-1$. Then for the chemical Baum-Welch Algorithm we have

[TABLE]

$\text{ for all }g,h\in H\times H\setminus\{(h^{*},h^{*})\}\text{ and }l=1,\ldots,L-1$ and \dot{\xi}_{lh^{*}h^{*}}\big{|}_{\Phi}=-\sum_{(g,h)\not=(h^{*},h^{*})}\dot{\xi}_{lgh}\big{|}_{\Phi}=0.

• •

$\theta^{\prime}_{gh}=\frac{\sum_{l=1}^{L-1}\xi^{\prime}_{l}(g,h)}{\sum_{l=1}^{L-1}\sum_{f\in H}\xi^{\prime}_{l}(g,f)}\text{ for all }g,h\in H$. Then for the chemical Baum-Welch Algorithm we have

[TABLE]

$\text{ for all }g\in H\text{ and }h\in H\setminus\{h^{*}\}$ and \dot{\theta}_{gh^{*}}\big{|}_{\Phi}=-\sum_{h\not=h^{*}}\dot{\theta}_{gh}\big{|}_{\Phi}=0.

• •

$\psi^{\prime}_{hw}=\frac{\sum_{l=1}^{L}\gamma^{\prime}_{l}(h)\delta_{w,v_{l}}}{\sum_{l=1}^{L}\gamma^{\prime}_{l}(h)}\text{ for all }h\in H\text{ and }w\in V$. Then for the chemical Baum-Welch Algorithm we have

[TABLE]

$\text{ for all }h\in H\text{ and }w\in V\setminus\{v^{*}\}$ and \dot{\psi}_{hw^{*}}\big{|}_{\Phi}=-\sum_{w\not=w^{*}}\dot{\psi}_{hw}\big{|}_{\Phi}=0.

So $\Phi$ is fixed point of the chemical Baum-Welch Algorithm. ∎

We will now prove Theorem 4.2. For the sake of convenience, we first recall the statement.

Theorem 4.2. Every positive fixed point for the Chemical Baum-Welch Algorithm on a Baum Welch Reaction system $(BW(\mathcal{M},h^{*},v^{*},L),k)$ with permissible rate $k$ is a fixed point for the Baum-Welch algorithm for the HMM $\mathcal{M}=(H,V,\theta,\psi,\pi)$.

Proof

Consider a positive point $\Phi=(\alpha^{\prime},\beta^{\prime},\gamma^{\prime},\xi^{\prime},\theta^{\prime},\psi^{\prime})$ with $\alpha^{\prime},\beta^{\prime},\gamma^{\prime}\in\mathbb{R}_{>0}^{L\times H}$, $\xi^{\prime}\in\mathbb{R}_{>0}^{(L-1)\times H\times H}$, $\theta^{\prime}\in\mathbb{R}_{>0}^{H\times H}$ and $\psi^{\prime}\in\mathbb{R}_{>0}^{H\times V}$. If $\Phi$ is a fixed point for the Chemical Baum-Welch Algorithm then we must have:

• •

\dot{\alpha}_{1h}\big{|}_{\Phi}=0 for all $h\in H$. This implies ${\alpha^{\prime}_{1h}}\times{\pi_{h^{*}}\psi^{\prime}_{h^{*}v_{1}}}={\alpha^{\prime}_{1h^{*}}}\times{\pi_{h}\psi^{\prime}_{hv_{1}}}\quad\text{ for all }h\in H\setminus{h^{*}}$. Since $\Phi$ is positive, this implies

[TABLE]

• •

\dot{\alpha}_{lh}\big{|}_{\Phi}=0 for all $h\in H$ and $l=2,\ldots,L$. This implies ${\alpha^{\prime}_{lh}}\times{\sum_{g\in H}\alpha^{\prime}_{l-1,g}\theta^{\prime}_{gh^{*}}\psi^{\prime}_{h^{*}v_{l}}}={\alpha^{\prime}_{lh^{*}}}\times{\sum_{g\in H}\alpha^{\prime}_{l-1,g}\theta^{\prime}_{gh}\psi^{\prime}_{hv_{l}}}\quad\text{ for all }h\in H\setminus\{h^{*}\}\text{ and }l=2,\ldots,L$. Since $\Phi$ is positive, this implies

[TABLE]

• •

\dot{\beta}_{lh}\big{|}_{\Phi}=0 for all $h\in H$ and $l=1,\ldots,L$. This implies ${\beta^{\prime}_{lh}}\times{\sum_{g\in H}\theta^{\prime}_{hg}\psi^{\prime}_{gv_{l+1}}\beta^{\prime}_{l+1,g}}={\beta^{\prime}_{lh^{*}}}\times{\sum_{g\in H}\theta^{\prime}_{h^{*}g}\psi^{\prime}_{gv_{l+1}}\beta^{\prime}_{l+1,g}}\quad\text{ for all }h\in H\setminus\{h^{*}\}\text{ and }l=1,\ldots,L-1$. Since $\Phi$ is positive, this implies

[TABLE]

• •

\dot{\gamma}_{lh}\big{|}_{\Phi}=0 for all $h\in H$ and $l=1,\ldots,L$. This implies ${\gamma^{\prime}_{l}(h)}\times{\alpha^{\prime}_{lh^{*}}\beta^{\prime}_{lh^{*}}}={\gamma^{\prime}_{l}(h^{*})}\times{\alpha^{\prime}_{lh}\beta^{\prime}_{lh}}\quad\text{ for all }h\in H\setminus\{h^{*}\}\text{ and }l=1,2,\dots,L-1$. Since $\Phi$ is positive, this implies

[TABLE]

• •

\dot{\xi}_{lgh}\big{|}_{\Phi}=0 for all $g,h\in H$ and $l=1,\ldots,L-1$. This implies ${\xi^{\prime}_{l}(g,h)}\times{\alpha^{\prime}_{lh^{*}}\theta^{\prime}_{h^{*}h^{*}}\psi^{\prime}_{h^{*}v_{l+1}}\beta^{\prime}_{l+1,h^{*}}}={\xi^{\prime}_{l}(h^{*},h^{*})}\times{\alpha^{\prime}_{lg}\theta^{\prime}_{gh}\psi^{\prime}_{hv_{l+1}}\beta^{\prime}_{l+1,h}}$ $\text{ for all }g,h\in H\times H\setminus\{(h^{*},h^{*})\}\text{ and }l=1,\ldots,L-1$. Since $\Phi$ is positive, this implies

[TABLE]

• •

\dot{\theta}_{gh}\big{|}_{\Phi}=0 for all $g,h\in H$. This implies ${\theta^{\prime}_{gh}}\times{\sum_{l=1}^{L-1}\xi^{\prime}_{l}(g,h^{*})}={\theta^{\prime}_{gh^{*}}}\times{\sum_{l=1}^{L-1}\xi^{\prime}_{l}(g,h)}\quad\text{ for all }g\in H\text{ and }h\in H\setminus\{h^{*}\}$. Since $\Phi$ is positive, this implies

[TABLE]

• •

\dot{\psi}_{hw}\big{|}_{\Phi}=0 for all $h\in H$ and $w\in V$. This implies ${\psi^{\prime}_{hw}}\times{\sum_{l=1}^{L}\gamma^{\prime}_{l}(h)\delta_{v^{*},v_{l}}}={\psi^{\prime}_{hv^{*}}}\times{\sum_{l=1}^{L}\gamma^{\prime}_{l}(h)\delta_{w,v_{l}}}\quad\text{ for all }h\in H\text{ and }w\in V\setminus\{v^{*}\}$. Since $\Phi$ is positive, $E_{lv}=\delta_{v,v_{l}}$ and $\sum_{v\in V}\delta_{v,v_{l}}=1$ this implies

[TABLE]

Because of the relaxation we get by Remark 1, the point $\Phi$ qualifies as a fixed point of the Baum-Welch algorithm. ∎

Appendix G.2 Rate of convergence analysis

In this section we will prove Theorem 4.3 and Theorem 4.4, but first we will state and prove two useful lemmas.

Lemma 1

Let $A$ be an arbitrary $n\times n$ matrix. Let $W$ be an $r\times n$ matrix comprising of $r$ linearly independent left kernel vectors of $A$ so that $WA=0_{r,n}$, where $0_{i,j}$ denotes a $i\times j$ matrix with all entries zero. Further suppose $W$ is in the row reduced form, that it is,

[TABLE]

where $I_{j}$ denotes the $j\times j$ identity matrix and $W^{\prime}$ is a $r\times(n-r)$ matrix. Let $A$ be given as

[TABLE]

where $A_{11}$ is a $(n-r)\times(n-r)$ matrix, $A_{12}$ is a $(n-r)\times r$ matrix, $A_{21}$ is a $r\times(n-r)$ matrix, and $A_{22}$ is a $r\times r$ matrix.

Then the $n-r$ eigenvalues (with multiplicity) of the matrix $A_{11}-A_{12}W^{\prime}$ are the same as the eigenvalues of $A$ except for $r$ zero eigenvalues.

Proof

Consider the $n\times n$ invertible matrix $P$ given by

[TABLE]

with determinant $\det(P)=\det(P^{-1})=1$. We have

[TABLE]

This implies that the characteristic polynomial of $A$ fulfils

[TABLE]

and the statement follows. ∎

Now we revisit the observation that every reaction in the Baum-Welch reaction network is a monomolecular transformation catalyzed by a set of species. For the purposes of our analysis, each reaction can be abstracted as a monomolecular reaction with time varying rate constants. This prompts us to consider the following monomolecular reaction system with $n$ species $X_{1},\ldots,X_{n}$ and $m$ reactions

[TABLE]

where $r_{j}\not=p_{j}$, and $r_{j},p_{j}\in\{1,\ldots,n\}$, and $k_{j}(t)$, $j=1,\ldots,m$, are mass-action reaction rate constants, possibly depending on time. We assume $k_{j}(t)>0$ for $t\geq 0$ and let $k(t)=(k_{1}(t),\ldots,k_{m}(t))$ be the vector of reaction rate constants. Furthermore, we assume there is at most one reaction $j$ such that $(r_{j},p_{j})=(r,p)\in\{1,\ldots,n\}^{2}$ and that the reaction network is strongly connected. The later means there is a reaction path from any species $X_{i}$ to any other species $X_{i^{\prime}}$. (In reaction network terms it means the network is weakly reversible and deficiency zero.)

The mass action kinetics of this reaction system is given by the ODE system

[TABLE]

Define the $n\times n$ matrix $A(t)=(a_{ii^{\prime}}(t))_{i,i^{\prime}=1,\ldots,n}$ by

[TABLE]

Then the ODE system might be written as

[TABLE]

Note that the column sums of $A(t)$ are zero, implying that $\sum_{i=1}^{n}x_{i}$ is conserved.

Lemma 2

Assume $k(t)$ for $t\geq 0$, converges exponentially fast towards $k=(k_{1},\ldots,k_{m})\in R^{m}_{>0}$ as $t\to\infty$, that is, there exists $\gamma_{1}>0$ and $K_{1}>0$ such that

[TABLE]

Let $A(t)$ be the matrix as defined in equation 1. And let $A$ be the matrix obtained with $k$ inserted for $k(t)$ in the matrix $A(t)$ that is, $A=\lim_{t\to\infty}A(t)$.

Then solutions of ODE system $\dot{x}=A(t)x$ starting at $x(0)\in\mathbb{R}^{n}_{\geq 0}$ converges exponentially fast towards the equilibrium $a\in\mathbb{R}^{n}_{>0}$ of the ODE system $\dot{x}=Ax$ starting at $x(0)\in\mathbb{R}^{n}_{\geq 0}$, that is, there exists $\gamma>0$ and $K>0$ such that

[TABLE]

Proof

We will first rephrase the ODE system such that standard theory is applicable. Let rank of $A$ be $n-r$. Let $W$ be as defined in Lemma 1, that is, $W$ be an $r\times n$ matrix comprising of $r$ linearly independent left kernel vectors of $A$ so that $WA=0$. Here since rank of $A$ is $n-r$, the rows of $W$ would form a basis for the left kernel of $A$. And as in Lemma 1, further suppose $W$ is in the row reduced form, that is,

[TABLE]

Then

[TABLE]

is a linear dynamical system with $r$ conservation laws (one for each row of $W$). Let $Wx(0)=T\in\mathbb{R}^{r}$ be the vector of conserved amounts. Let $\hat{x}=(x_{1},\ldots,x_{n-r})$ and $\tilde{x}=(x_{n-r+1},\ldots,x_{n})$. We will consider the (equivalent) dynamical system in which $r$ variables are eliminated, expressed through the conservation laws

[TABLE]

As in Lemma 1, let $A$ be given as

[TABLE]

where $A_{11}$ is a $(n-r)\times(n-r)$ matrix, $A_{12}$ is a $(n-r)\times r$ matrix, $A_{21}$ is a $r\times(n-r)$ matrix, and $A_{22}$ is a $r\times r$ matrix. This yields

[TABLE]

with $C=A_{11}-A_{12}W^{\prime}$ and $D=A_{12}$. We call this as the reduced ODE system. Note that this reduced system has only $n-r$ variables and that the conservation laws are built directly into it. This implies that the differential equation changes if $T$ is changed. The role of this construction is so that we can work only with the non-zero eigenvalues of the $A$.

As we also have $WA(t)=0$ for all $t\geq 0$, the ODE $\dot{x}=A(t)x$ can also be similarly reduced to

[TABLE]

with $C(t)=A_{11}(t)-A_{12}(t)W^{\prime}$ and $D(t)=A_{12}(t)$, where analogous to $A_{11}$, we define $A_{11}(t)$ to be the top-left $(n-r)\times(n-r)$ sub-matrix of $A(t)$ and analogous to $A_{12}$, we define $A_{12}(t)$ to be the top-right $(n-r)\times r$ sub-matrix of $A(t)$.

Now if $a$ is the equilibrium of the ODE system $\dot{x}=Ax$ starting at $x(0)$, then $\hat{a}=(a_{1},\ldots,a_{n-r})$ is an equilibrium of the reduced ODE sytem $\dot{\hat{x}}=C\hat{x}+DT$ starting at $\hat{x}(0)=(x_{1}(0),\ldots,x_{n-r}(0))$. Suppose we are able to prove that solutions of reduced ODE $\dot{\hat{x}}=C(t)\hat{x}+D(t)T$ starting at $\hat{x}(0)$ converges exponentially fast towards $\hat{a}$ then because of the conservation laws $\tilde{x}=T-W^{\prime}\hat{x}$, we would also have that solutions of $\dot{x}=A(t)x$ starting at $x(0)$ converges exponentially fast towards $a$. So henceforth, we will work only with the reduced ODE systems. For notational convinience, we will drop the hats off $\hat{x}$ and $\hat{a}$ and simply refer to them as $x$ and $a$ respectively.

By subtracting and adding terms to the reduced ODE system (in equation 5), we have

[TABLE]

As $a$ is an equilibrium of the ODE system $\dot{x}=Cx+DT$, we have $Ca+DT=0$.

Define $y=x-a$. Then

[TABLE]

where it is used that $Ca+DT=0$, and $E(t)=(C(t)-C)x(t)+(D(t)-D)T$.

The solution to the above ODE system is known to be

[TABLE]

We have, using (6) and the triangle inequality,

[TABLE]

Now $A$ as defined (see equation 1 with $k$ inserted for $k(t)$) would form a Laplacian matrix over a strongly connected graph and so it follows that all the eigenvalues of $A$ are either zero or have negative real part. And using $C=A_{11}-A_{12}W^{\prime}$ and Lemma 1 it follows that all eigenvalues of $C$ have negative real part. Hence it follows that

[TABLE]

where $0<\gamma_{2}<-\Re(\lambda_{1})$ and $K_{2}>0$. Here $\lambda_{1}$ is the eigenvalue of $C$ with the largest real part.

The matrices $C(t)$ and $D(t)$ are linear in $k(t)$. And as $k(t)$ converges exponentially fast towards $k$, it follows that the matrices $C(t)$ and $D(t)$ converge exponentially fast towards $C$ and $D$ respectively. Hence it follows that

[TABLE]

where

• •

$\mathinner{\!\left\lVert(C_{0}(t)-C)\right\rVert}\mathinner{\!\left\lVert x(t)\right\rVert}\leq K_{3}e^{-\gamma_{3}t}$ for some $K_{3},\gamma_{3}>0$ as $C(t)$ converges exponentially fast towards $C$ and $x(t)$ is bounded (as in the original ODE $\sum_{i=1}^{n}x_{i}$ is conserved), and

• •

$\mathinner{\!\left\lVert(D_{0}(t)-D)\right\rVert}\mathinner{\!\left\lVert T\right\rVert}\leq K_{4}e^{-\gamma_{4}t}$ for some $K_{4},\gamma_{4}>0$ as $D(t)$ converges exponentially fast towards $D$, and

• •

$K_{5}=\frac{1}{2}\max(K_{3},K_{4})$ and $\gamma_{5}=\min(\gamma_{3},\gamma_{4})$.

Collecting all terms we have for all $t\geq 0$,

[TABLE]

by choosing $K_{0}=\max\left(K_{2}K_{5},\mathinner{\!\left\lVert y(0)\right\rVert}K_{2}\right)$ and $\gamma_{0}=\min(\gamma_{1},\gamma_{2})$. In the last line $\gamma$ is chosen such that $0<\gamma<\gamma_{0}$ and $K$ is sufficiently large. Since $y(t)=x(t)-a$ we have,

[TABLE]

as required.

∎

We will now prove Theorem 4.3. For the sake of convenience, we first recall the statement.

Theorem 4.3. For the Baum Welch Reaction System $(BW(\mathcal{M},h^{*},v^{*},L),k)$ with permissible rates $k$, if the concentrations of $\theta$ and $\psi$ species are held fixed at a positive point then the Forward, Backward and Expection step reaction systems on $\alpha,\beta,\gamma$ and $\psi$ species converge to equilibrium exponentially fast.

Proof

It follows by repeated use of Lemma 2. For $l=1$ the forward reaction network can be interpretated as the following molecular reactions:

[TABLE]

for all $h\in H\setminus\{h^{*}\}$ and $w\in V$, as they are dynamically equivalent. Here the effective rate constants ($\pi_{h^{*}}\psi_{h^{*}w}E_{1w}$ or $\pi_{h}\psi_{hw}E_{1w}$) are independent of time and so the conditions of Lemma 2 are fulfilled. Thus this portion of the reaction network converges exponentially fast.

The rest of the forward reaction network can be similarly interpretated as the following molecular reactions:

[TABLE]

for all $g\in H$, $h\in H\setminus\{h^{*}\},l=2,\ldots,L$ and $w\in V$. For layers $l=2,\ldots,L$ of the forward reaction network we observe that the effective rate constants ($\alpha_{l-1,g}\theta_{gh^{*}}\psi_{h^{*}w}E_{lw}$ or $\alpha_{l-1,g}\theta_{gh}\psi_{hw}E_{lw}$) for layer $l$ depend on time only through $\alpha_{l-1,g}$. If we suppose that the concentration of $\alpha_{l-1,g}$ converges exponentially fast, then we can use Lemma 2 to conclude that the concentration of $\alpha_{lh}$ also converges exponentially fast. Thus using Lemma 2 inductively layer by layer we conclude that forward reaction network converges exponentially fast. The backward reaction network converges exponentially fast, similarly.

For the expectation reaction network it likewise follows by induction. But here, notice if we interpreted the expectation network similarly into molecular reactions, the effective rate constants would depend on time through the products such as $\alpha_{lh}\beta_{lh}$ or $\alpha_{lh}\beta_{l+1,h}$. So to apply Lemma 2 we need the following: If $\alpha_{l}(t)$ and $\beta_{l}(t)$ converge exponentially fast towards $a_{l}$ and $b_{l}$ then the product $\alpha_{l}(t)\beta_{l}(t)$ converges exponentially fast towards $a_{l}b_{l}$ as

[TABLE]

where $K$ is some suitably large constant. We can further observe that $\alpha_{l}(t)\beta_{l}(t)$ converges exponentially fast towards $a_{l}b_{l}$ at rate $\gamma=\min(\gamma_{a},\gamma_{b})$, where $\gamma_{a}$ and $\gamma_{b}$, respectively, are the exponential convergence rates of $\alpha_{l}(t)$ and $\beta_{l}(t)$. A similar argument goes for the products of the form $\alpha_{l}(t)\beta_{l+1}(t)$. And thus the expectation reaction network, also converges exponentially fast. ∎

We will now prove Theorem 4.4. For the sake of convenience, we first recall the statement.

Theorem 4.4. For the Baum Welch Reaction System $(BW(\mathcal{M},h^{*},v^{*},L),k)$ with permissible rates $k$, if the concentrations of $\alpha,\beta,\gamma$ and $\xi$ species are held fixed at a positive point then the Maximization step reaction system on $\theta$ and $\psi$ converges to equilibrium exponentially fast.

Proof

Exponential convergence of the maximisation network follows by a similar layer by layer inductive use of Lemma 2. ∎