A Reaction Network Scheme Which Implements the EM Algorithm

A detailed algorithmic explanation is required for how a network of chemical reactions can generate the sophisticated behavior displayed by living cells. Though several previous works have shown that reaction networks are computationally universal and can in principle implement any algorithm, there is scope for constructions that map well onto biological reality, make efficient use of the computational potential of the native dynamics of reaction networks, and make contact with statistical mechanics. We describe a new reaction network scheme for solving a large class of statistical problems. Specifically we show how reaction networks can implement information projection, and consequently a generalized Expectation-Maximization algorithm, to solve maximum likelihood estimation problems in partially-observed exponential families on categorical data. Our scheme can be thought of as an algorithmic interpretation of E. T. Jaynes’s vision of statistical mechanics as statistical inference.

1. 1.

   This situation can arise because only a linear projection is observable. It can also happen because we require a rich family of probability distributions on the space of s points, but don’t want to give away the nice properties of exponential families. We can achieve both by imagining that our observation s comes from projection from a data vector x living in a higher-dimensional space, and then employ an exponential family of probability distributions on this higher-dimensional space.

Authors and Affiliations

1. India Institute of Technology Bombay, Mumbai, India

   Muppirala Viswa Virinchi, Abhishek Behera & Manoj Gopalkrishnan

Corresponding author

Correspondence to Manoj Gopalkrishnan .

Editors and Affiliations

1. University of California, Davis, CA, USA

   David Doty

2. Technical University Munich, Garching, Germany

   Hendrik Dietz

Appendix A

Proof

(Proof of Theorem 2). (1) Fix \(\alpha \in \mathbb {R}^n_{>0}\). We first prove uniqueness: suppose for contradiction that there are at least two points of intersection \(\alpha ^*_1,\alpha ^*_2\) of the polytope \((\alpha + V^\perp )\cap \mathbb {R}^n_{\ge 0}\) with the hypersurface \(e^V\). Since \(\alpha ^*_1,\alpha ^*_2\in e^V\), we have \(\log \alpha _1^* - \log \alpha _2^*\in V\). Since \(\alpha - \alpha _1^*\in V^\perp \), we have \((\alpha -\alpha _1^*)\cdot (\log \alpha _1^* - \log \alpha _2^*)=0\). Then by the Pythagorean theorem, \(D(\alpha \Vert \alpha ^*_1) = D(\alpha \Vert \alpha ^*_2) + D(\alpha ^*_2\Vert \alpha ^*_1)\) which implies \(D(\alpha \Vert \alpha ^*_1)\ge D(\alpha \Vert \alpha ^*_2)\). By a symmetric argument, \(D(\alpha \Vert \alpha ^*_2)\ge D(\alpha \Vert \alpha ^*_1)\) and we conclude \(D(\alpha \Vert \alpha ^*_2)= D(\alpha \Vert \alpha ^*_1)\). In particular, \(D(\alpha ^*_2\Vert \alpha ^*_1)=0\) which implies \(\alpha ^*_1=\alpha ^*_2\) by Note 1.

To prove that there exists at least one point of intersection, and to show (2), fix \(\beta \in e^V\). We will show that the E-Projection \(\alpha ^*\) of \(\beta \) to \((\alpha +V^\perp )\cap \mathbb {R}^n_{\ge 0}\) belongs to \(e^V\). This point \(\alpha ^*\) exists since \(D(x\Vert \beta )\) is continuous, and hence attains its minimum over the compact set \((\alpha +V^\perp )\cap \mathbb {R}^n_{\ge 0}\). Further, because \(\alpha ^*\) is an infimum, we need that \(\lim _{\lambda \rightarrow 0} \frac{d f((1-\lambda )\alpha ^* + \lambda \alpha )}{d\lambda }=0\). That is, \((\alpha - \alpha ^*)\log \frac{\alpha ^*}{\beta }=0\), which implies that \(\alpha ^*\in e^V\) since \(\alpha \) could have been replaced by any other arbitrary point of \((\alpha +V^\perp )\cap \mathbb {R}^n_{>0}\).

(3) now follows because \(\alpha ^*\in e^V\) implies \(D(\alpha \Vert \alpha ^*) + D(\alpha ^*\Vert \beta ) = D(\alpha \Vert \beta )\) for all \(\beta \in e^V\), hence \(\alpha ^*\) is the M-Projection of \(\alpha \) to \(e^V\).

Proof

(Proof of Theorem 6).

1. (1)

   From the chain rule, \(\dot{D}(x(t) \Vert y_A\circ \theta (t)) = ( \nabla _x D\cdot \dot{x} + \nabla _\theta D\cdot \dot{\theta })|_{(x(t),\theta (t))}\). From Theorem 4, the first term is nonpositive with equality iff x(t) is the E-Projection of \(y(\theta (t))\) onto \(H_{x_0}\). From Theorem 5, the second term is nonpositive with equality iff \(y(\theta (t))\) is the M-Projection of x(t) onto \(y_A(\mathbb {R}^m)\). Hence \(dD(x(t) \Vert y_A\circ \theta (t))/dt \le 0\) with equality iff both x(t) is the E-Projection of \(y_A\circ \theta (t)\) to \(H_{x_0}\) and \(y_A\circ \theta (t)\) is the M-Projection of x(t) to \(y_A(\mathbb {R}^m)\).

2. (2)

   Since \(D(x(t) \Vert y_A\circ \theta (t))\) has a lower bound and \(dD(x(t) \Vert y_A\circ \theta (t))/dt \le 0\), eventually \(dD(x(t) \Vert y_A\circ \theta (t))/dt = 0\) at which point, by the above argument, both the E-Projection and M-Projection subnetworks are stationary, so that \(\dot{x}=0\) and \(\dot{\theta }=0\). Hence the limit \((\hat{x},\hat{\theta })=\lim _{t\rightarrow \infty }(x(t), \theta (t))\) exists.

3. (3)

   follows since \(\nabla _\theta D(x\Vert y_A(\theta ))|_{\hat{x},\hat{\theta }} = \dot{\theta }(t)/\theta (t)|_{\hat{x},\hat{\theta }}=0\) when \(\hat{\theta }\in \mathbb {R}^\varTheta _{>0}\).

4. (4)

   \(\nabla _x D(x\Vert y_A(\theta ))|_{\hat{x},\hat{\theta }} = \log \left( \frac{\hat{x}}{y_A(\hat{\theta })}\right) \). By (1), the point \(\hat{x}\) is the E-Projection of \(y_A(\hat{\theta })\) to \(H_{x_0}\). Hence by Theorem 2, the point \(\hat{x}\) is the Birch point of \(x_0\) relative to the affine space \(\log y_A(\mathbb {R}^m)\), so that \((x - \hat{x})\log \left( \frac{\hat{x}}{y_A(\hat{\theta })}\right) =0\) for all \(x\in H_{x_0}\). Hence the gradient \(\nabla _x D(x\Vert y_A(\theta ))|_{\hat{x},\hat{\theta }}\) is perpendicular to \(H_{R_\mathcal {B}}\).

Example 8

Consider \(A=\left( \begin{array}{ccc}2&{}1&{}0\\ 0&{}1&{}2\end{array}\right) \) and \(\mathcal {S}=\left( \begin{array}{rrr}1&{}0&{}1\\ 1&{}1&{}1\end{array}\right) \). The vector \(\left( \begin{array}{r}1 \\ 0 \\ -1 \end{array}\right) \) spans \(\ker \mathcal {S}\). The corresponding reaction network is

Here the concentration of \(X_2\) remains invariant with time. Let c be the initial concentration of \(X_2\). If \(c < 1/3\) then the system admits two stable equilibria and one unstable equilibrium. The points \((y_1,c,y_2,\sqrt{y_1},\sqrt{y_2})\) and \((y_2,c,y_1,\sqrt{y_2},\sqrt{y_1})\) are the stable equilibria where \(y_1=\frac{1-c}{2}+\frac{\sqrt{(1-3c)(1+c)}}{2}\) and \(y_2=\frac{1-c}{2}-\frac{\sqrt{(1-3c)(1+c)}}{2}\), and \(\left( \frac{1-c}{2},c,\frac{1-c}{2},\sqrt{\frac{1}{3}},\sqrt{\frac{1}{3}}\right) \) is the unstable equilibrium. On the other hand, if \(c\ge 1/3\) then there is only one equilibrium point at \(\left( \frac{1-c}{2},c,\frac{1-c}{2},\sqrt{\frac{1}{3}},\sqrt{\frac{1}{3}}\right) \), and this point is stable.

Example 9

Consider \(A=\left( \begin{array}{ccc}2&{}1&{}0\\ 0&{}1&{}2\\ \end{array}\right) \) and \(\mathcal {S}=\left( \begin{array}{rrr}1&{}-1&{}0\\ 1&{}1&{}1\end{array}\right) \). The vector \(\left( \begin{array}{r}1 \\ 1 \\ -2 \end{array}\right) \) spans \(\ker \mathcal {S}\). The corresponding reaction network is

Here the set \(\{X_1,X_2,\theta _1\}\) is a critical siphon. If we start at the initial concentrations \(x_1=0.05,x_2=0.05,x_3=0.9,\theta _1=0.1,\theta _2=1.0\), then the system converges to \(x_1=0,x_2=0,x_3=1,\theta _1=0,\theta _2=1\), hence this system is not persistent. This provides one explanation for this data: all the outcomes were of type \(X_3\). If instead we start at \(\theta _1=0.5,\theta _2=1.0\) and the same x concentrations, then the system converges to \(x_1=x_2=x_3=1/3,\theta _1=\theta _2=1/\sqrt{c}\). This provides a different explanation for the same data: all three outcomes have occurred equally frequently.

Example 10

Boltzmann machines are a popular model in machine learning. Formally a Boltzmann machine is a graph \(G=(V,E)\), each of whose nodes can be either 1 or 0. One associates to every configuration \(s \in \{0,1\}^V\) of the Boltzmann machine an energy \(E(s)=-\sum _i b_i s_i-\sum _{ij}w_{ij}s_is_j\). The probability of the Boltzmann machine being in configuration s is given by the exponential family \(P(s;b,w)\propto \exp (-E(s))\). Boltzmann machines can be used to do inference conditioned on partial observations, and learning of the maximum likelihood values of the parameters \(b_i,w_{ij}\) can be done by a stochastic gradient descent.

Our EM scheme can be used to implement the learning rule of arbitrary Boltzmann machines in chemistry. We illustrate the construction on the 3-node Boltzmann machine with \(V=\{x_1,x_2,x_3\}\):

with biases \(b_1,b_2,b_3\) and weights \(w_{12}\) and \(w_{23}\). We will work with parameters \(\theta _i=\exp (b_i)\) and \(\theta _{ij}=\exp (w_{ij})\). The design matrix \(A=(a_{ij})_{5\times 8}\) is

and the corresponding exponential model \(y_A:\mathbb {R}^5\rightarrow \mathbb {R}^8_{>0}\) sends \(\theta = (\theta _1,\theta _2,\theta _3,\theta _{ 12},\theta _{23}) \longmapsto (\theta ^{a_{.1}},\theta ^{a_{.2}},\dots ,\theta ^{a_{.8}})\). If the node \(x_2\) is hidden then the observation matrix \(\mathcal S\) is

Our EM scheme yields the reaction network

Suppose we observe a marginal distribution (0.24, 0.04, 0.17, 0.55) on the visible nodes \(x_1,x_3\). To solve for the maximum likelihood \(\hat{\theta }\), we can initialize the system with \(X_{000}=0.24,X_{001}=0.04,X_{010}=0,X_{011}=0,X_{100}=0.17,X_{101}=0.55,X_{110}=0,X_{111}=0\) and all \(\theta \)’s initialized to 1, the system reaches steady state at \(\hat{\theta }_1=0.5176,\hat{\theta }_2=0.0018,\hat{\theta }_3=0.3881,\hat{\theta }_{12}=0.8246,\hat{\theta }_{23}=0.7969, \hat{X}_{000}=0.2391,\hat{X}_{001}=0.0389, \hat{X}_{010}=0.0009, \hat{X}_{011}=0.011, \hat{X}_{100}=0.1695, \hat{X}_{101}=0.5487, \hat{X}_{110}=0.0005, \hat{X}_{111}=0.0013\).

Reprints and permissions

© 2018 Springer Nature Switzerland AG

Policies and ethics