A Scheme for Molecular Computation of Maximum
Likelihood Estimators for Log-Linear Models

arXiv:1506.03172v2 [] 10 Jun 2016

Manoj Gopalkrishnan
Tata Institute of Fundamental Research, Mumbai 400 005, India
manoj.gopalkrishnan@gmail.com, http://www.tcs.tifr.res.in/~ manoj

Abstract. We propose a novel molecular computing scheme for statistical inference. We focus
on the much-studied statistical inference problem of computing maximum likelihood estimators for log-linear models. Our scheme takes log-linear models to reaction systems, and the
observed data to initial conditions, so that the corresponding equilibrium of each reaction
system encodes the corresponding maximum likelihood estimator. The main idea is to exploit
the coincidence between thermodynamic entropy and statistical entropy. We map a Maximum
Entropy characterization of the maximum likelihood estimator onto a Maximum Entropy
characterization of the equilibrium concentrations for the reaction system. This allows for an
efficient encoding of the problem, and reveals that reaction networks are superbly suited to
statistical inference tasks. Such a scheme may also provide a template to understanding how in
vivo biochemical signaling pathways integrate extensive information about their environment
and history.

1

Introduction

The sophisticated behavior of cells emerges from the computations that are being performed by
the underlying biochemical reaction networks. These biochemical pathways have been studied in
a “top-down” manner, by looking for recurring motifs, and signs of modularity [18]. There is also
an opportunity to study these pathways in a “bottom-up” manner by proposing primitive building
blocks which can be composed to create interesting and technologically valuable behavior. This
“bottom-up” approach connects with work in the Molecular Computation community whose goal
is to generate sophisticated behavior using DNA hybridization reactions [23,24,19,31,27,6,7,22,3,25]
and other Artificial Chemistry approaches [5,10].
We propose a new building block for molecular computation. We show that the mathematical
structure of reaction networks is particularly well adapted to compute Maximum Likelihood Estimators for log-linear models, allowing a pithy encoding of such computations by reactions. According
to [12]:
Log-linear models are arguably the most popular and important statistical models for
the analysis of categorical data; see, for example, Bishop, Fienberg and Holland (1975) [4],
Christensen (1997) [8]. These powerful models, which include as special cases graphical models [see, e.g., Lauritzen (1996) [16]] as well as many logit models [see, e.g., Agresti (2002) [1],
Bishop, Fienberg and Holland (1975) [4]], have applications in many scientific areas, ranging
from social and biological sciences, to privacy and disclosure limitation problems, medicine,
data mining, language processing and genetics. Their popularity has greatly increased in the
last decades...

In order to respond in a manner that maximizes fitness, a cell has to correctly estimate the
overall state of its environment. Receptors that sit on cell walls collect a large amount of information
about the cellular environment. Processing and integration of this spatially and temporally extensive
and diverse information is carried out in the biochemical reaction pathways. We propose that this
processing and integration may be advantageously viewed from the lens of machine learning.
Our proposal entails that schemes for statistical inference by reaction networks are of biological significance, and are deserving of as thorough and extensive a study as schemes for statistical
inference by neural networks. In particular, machine learning is not just a tool for the analysis of
biochemical data, but theoretical and technological insights from machine learning could provide a
deep and fundamental way, and perhaps “the” correct way, to think about biochemical networks.
We view the scheme we present here as a promising first step in this program of applying machine
learning insights to biochemical networks.
The problem: We illustrate the main ideas of our scheme with an example. Following [21], consider

the log-linear model (also known as toric model) described by the design matrix A = 20 11 02 .
This means that we are observing an event with three possible mutually exclusive outcomes, call
them X1 , X2 , and X3 , which represent respectively the columns of A. The rows of A represent
“hidden variables” θ1 and θ2 respectively which parametrize the statistics of the outcomes in the
following way specified by the columns of A:
P [X1 | θ1 , θ2 ] ∝ θ12
P [X2 | θ1 , θ2 ] ∝ θ1 θ2
P [X3 | θ1 , θ2 ] ∝ θ22
where the constant of proportionality normalizes the probabilities so they sum to 1. 1
Suppose several independent trials are carried out, and the outcome X1 is observed x1 ∈ (0, 1)
fraction of the time, the outcome X2 is observed x2 ∈ (0, 1 − x1 ) fraction of the time, and the
outcome X3 is observed x3 = 1 − x1 − x2 fraction of the time. We wish to find the maximum
likelihood estimator (θ̂1 , θ̂2 ) ∈ R2>0 of the parameter (θ1 , θ2 ), i.e., that value of θ which maximizes
the likelihood of the observed data.
Our contribution: We describe a scheme that takes the design matrix A to a reaction network
that solves the maximum likelihood estimation problem. In Definition 8, we describe our scheme for
every matrix A over the integers with all column sums equal. All our results hold in this generality.
– In Definition 8.5, we show how to obtain from the matrix A, a reaction network that computes
the maximum likelihood distribution. Specialized to our example, note that the kernel of the
matrix A is spanned by the vector (1, −2, 1)T . We encode this by the reversible reaction
1

−
⇀
X1 + X3 ↽
− 2X2
1

– In Theorem 5, we show that if this reversible reaction is started at initial concentrations X1 (0) =
x1 , X2 (0) = x2 , X3 (0) = x3 , and the dynamics proceeds according to the law of mass action with
all specific rates set to 1:
Ẋ1 (t) = Ẋ3 (t) = −X1 (t)X3 (t) + X22 (t),
1

Ẋ2 (t) = −2X22(t) + 2X1 (t)X3 (t)

It is more common in statistics and statistical mechanics literature to write θ1 = e−E1 and θ2 = e−E2 in
terms of “energies” E1 , E2 so that P [X2 | E1 , E2 ] ∝ e−E1 −E2 for example.

2

then the reaction reaches equilibrium (x̂1 , x̂2 , x̂3 ) where x̂1 + x̂2 + x̂3 = 1 and x̂1 ∝ θ̂12 , x̂2 ∝
θ̂1 θ̂2 , and x̂3 ∝ θ̂22 , so that (x̂1 , x̂2 , x̂3 ) represents the probability distribution over the outcomes
X1 , X2 , X3 at the maximum likelihood θ̂1 , θ̂2 .
– This part of our scheme involves only reversible reactions, and requires no catalysis (see [13,
Theorem 5.2] and Lemma 2). One difficulty with implementing such schemes has been that
empirical control over kinetics is rather poor. Exquisitely setting the specific rates of individual
reactions to desired values is very tricky, and requires a detailed understanding of molecular
dynamics. Our scheme avoids this problem since any choice of specific rates that leads to the
same equilibrium will do. Hence we can freely set the specific rates so long as the equilibrium
constants (ratio of forward and backward specific rates) have value 1. This is an equilibrium
thermodynamic condition that is much easier to ensure in vitro. This combination of reversible
reactions, no catalysis, and robustness to the values of the specific rates may make this scheme
particularly easy and efficient to implement.
– In Definition 8.2, we show how to obtain from the matrix A a reaction network that computes
the maximum likelihood estimator. Specialized to our example, we obtain the reaction network
with 5 species X1 , X2 , X3 , θ1 , θ2 and the 5 reactions:
X1 + X3 ⇋ 2X2 ,

2θ1 → 0,
θ1 + θ2 → 0,

X1 → X1 + 2θ1 ,
X 2 → X 2 + θ1 + θ2 .

The number of species equals the number of rows plus the number of columns of A. The reactions
are not uniquely determined by the problem, but become so once we choose a basis for the kernel
of A and a maximal linearly independent set of columns. Here we have chosen columns 1 and 2.
Each column of A determines a pair of irreversible reactions.
– Theorem 6 implies that if this reaction system is launched at initial concentrations X1 (0) =
x1 , X2 (0) = x2 , X3 (0) = x3 and arbitrary concentrations of θ1 (0) and θ2 (0), and the dynamics
proceeds according to the law of mass action with all specific rates set to 1:
Ẋ1 (t) = Ẋ3 (t) = −X1 (t)X3 (t) + X22 (t),

θ̇1 (t) = −2θ12 (t) + 2X1 (t) − θ1 (t)θ2 (t) + X2 (t),

Ẋ2 (t) = −2X22 (t) + 2X1 (t)X3 (t),

θ̇2 (t) = −θ1 θ2 (t) + X2 (t),

then the reaction reaches equilibrium (xˆ1 , xˆ2 , xˆ3 , θ̂1 , θ̂2 ) where (θ̂1 , θ̂2 ) is the maximum likelihood
estimator for the data frequency vector (x1 , x2 , x3 ) and (x̂1 , x̂2 , x̂3 ) represents the probability
distribution over the outcomes X1 , X2 , X3 at the maximum likelihood. We prove global convergence: our dynamical system provably converges to the desired equilibrium. Global convergence
results are known to be notoriously hard to prove in reaction network theory [14].
– A number of schemes have been proposed for translating reaction networks into DNA strand
displacement reactions [27,22,6,7]. Adapting these schemes to our setting should allow molecular
implementation of our MLE-solving reaction networks with DNA molecules.

2

Maximum Likelihood Estimation in toric models

The definitions and results in this section mostly follow [21]. Because we require a slightly stronger
statement, and Theorem 1 allows a short, easy, and insightful proof, we give the proof here for
completeness.
3

In statistics, a parametric model consists of a family of probability distributions, one for each
value of the parameters. This can be described as a map from a manifold of parameters into a
manifold of probability distributions. If this map can be described by monomials as below, then the
parametric statistical model is called a toric or log-linear model, as we now describe.
Definition 1 (Toric Model). Let m, n be positive integers. The probability simplex and its relative
interior are:
∆n := {(x1 , x2 , . . . , xn ) ∈ Rn≥0 | x1 + x2 + · · · + xn = 1}
ri(∆n ) := {(x1 , x2 , . . . , xn ) ∈ Rn>0 | x1 + x2 + · · · + xn = 1}.
P
An m × n matrix A = (aij )m×n of integer entries is a design matrix iff all its column sums i aij
a
a
a
are equal. Let aj := (a1j , a2j , . . . , amj )T be the j’th column of A. Define θaj := θ1 1j θ2 2j . . . θmmj .
an
a2
m
a1
Define the parameter space Θ := {θ ∈ R>0 | θ + θ + · · · + θ = 1}. The toric model of A is
the map
pA = (p1 , p2 , . . . , pn ) : Θ → ∆n given by pj (θ) = θaj for j = 1 to n.
We could also have defined the parameter space Θ to be all of Rm
>0 , in which case we would need
to normalize the probabilities by the partition function θa1 + θa2 + · · · + θan to make sure they add
up to 1. For our present purposes, the current approach will prove technically more direct.
Note that here pj (θ) specifies Pr[j | θ], the conditional probability of obtaining outcome j given
that the true state of the world is described by θ.
A central problem of statistical inference is the problem of parameter estimation. After performing several independent identical trials, suppose the data vector u ∈ Zn≥0 is obtained as a
record of how many times each outcome occurred. Let the norm |u|1 := u1 + u2 + · · · + un denote
the total number of trials performed. The Maximum Likelihood solution to the problem of parameter estimation finds that value of the parameter θ which maximizes the likelihood function
fu (θ) := Pr[u | θ], i.e.:
θ̂(u) := arg sup fu (θ)

(1)

θ∈Θ

is a maximum likelihood estimator or MLE for the data vector u. We will call the point p̂(u) :=
pA (θ̂(u)) a maximum likelihood distribution.
Definition 2. Let A be an m × n design matrix, and u a data vector. Then the sufficient polytope
u
is PA (u) := {p ∈ ri(∆n ) | Ap = A |u|
}.
1
The following theorem is a version of Birch’s theorem from Algebraic Statistics. It provides a
variational characterization of the maximum likelihood distribution as the unique maximum entropy
distribution in the sufficient polytope. In particular the maximum likelihood distribution always
belongs to the sufficient polytope, which justifies the name.
Theorem 1. Fix a design matrix A of size m × n.
1. If u, v ∈ Zn≥0 are nonzero data vectors such that Au/|u|1 = Av/|v|1 then they have the same
maximum likelihood estimator: θ̂(u) = θ̂(v).
2. Further if PA (u) is nonempty then
Pn
(a) There is a unique distribution p̃ ∈ PA (u) which maximizes Shannon entropy H(p) = − i=1 pi log pi
viewed as a real-valued function from the closure PA (u) of PA (u) with 0 log 0 defined as 0.
4

(b) {p̃} = PA (u) ∩ pA (Θ).
(c) p̃ = p̂(u), the Maximum Likelihood Distribution for the data vector u.
Proof. 1. Fix a data vector u. Note that fu (θ) =
Therefore the maximum likelihood estimator

|u|1 !
u1
u2
u1 !u2 !...un ! p1 (θ) p2 (θ)

. . . pn (θ)un =

|u|1 !
Au
.
u1 !u2 !...un ! θ

θ̂(u) = arg sup θAu = arg sup (θAu )1/|u|1 = arg sup θAu/|u|1
θ∈Θ

θ∈Θ

θ∈Θ

where the second equality is true because the function x 7→ xc is monotonically increasing whenever
c > 0. It follows that if v ∈ Zn≥0 is a data vector such that Au/|u|1 = Av/|v|1 then θ̂(u) = θ̂(v).
2.(a) Suppose PA (u) is nonempty. A local maximum of the restriction H|PA (u) of H to the polytope
PA (u) can not be on the boundary ∂PA (u) because for p ∈ ∂PA (u), moving in the direction of
arbitrary q ∈ PA (u) increases H, as can be shown by a simple calculation:
lim

λ→0

d
H((1 − λ)p + λq) → +∞.
dλ

Since H is a continuous function and the closure PA (u) is a compact set, H must attain its maximum
value in PA (u). Further H is a strictly concave function since its Hessian is diagonal with entries
−1/pi and hence negative definite. It follows that H|PA (u) is also strictly concave, and has a unique
local maximum at p̃ ∈ PA (u), which is also the global maximum.
(b) By concavity of H, the maximum p̃ is the unique point in PA (u) such that ∇H(p̃) is perpendicular to PA (u). We claim that q ∈ PA (u) ∩ pA (Θ) iff ∇H(q) = (−1 − log q1 , −1 − log q2 , . . . , −1 − log qn )
is perpendicular to PA (u). Since all column sums are equal, this is equivalent to requiring that log q
be in the span of the rows of A, which is true iff q ∈ pA (Θ). Hence PA (u) ∩ pA (Θ) = {p̃}.
(c) To compute the Maximum Likelihood Distribution p̂(u), we proceed as follows:
p̂(u) = pA (θ̂(u)) = pA (arg sup θAu ) = pA (arg sup θAu/|u|1 )
θ∈Θ

= pA (arg sup θAp̃ ) = arg
θ∈Θ

θ∈Θ

sup
p∈pA (Θ)

pp̃ = arg

sup

n
X

p̃i log pi = p̃

p∈pA (Θ) i=1

Pn
where the fourth equality uses Ap̃ = Au/|u|1 and the last equality follows because i=1 p̃i log pi
viewed as a function of p attains its maximum in all of ∆n , and hence in pA (Θ), at p = p̃.
This theorem already exposes the core of our idea. We will design reaction systems that maximize entropy subject to the “correct” constraints capturing the polytope PA (u). Then because the
reactions also proceed to maximize entropy, the equilibrium point of our dynamics will correspond
to the maximum likelihood distribution. Most of the technical work will go in proving convergence
of trajectories to these equilibrium points.

3

Reaction Networks

According to [20], “In building a design theory for chemistry, chemical reaction networks are usually
the most natural intermediate representation - the middle of the hourglass [11]. Many different high
5

level languages and formalisms have been and can likely be compiled to chemical reactions, and
chemical reactions themselves (as an abstract specification) can be implemented with a variety of
low level molecular mechanisms.”
In Subsection 3.1, we recall the definitions and results for reaction networks which we will need
for our main results. For a comprehensive presentation of these ideas, see [13]. In Subsection 3.2, we
prove a new result in reaction network theory. We extend a previously known global convergence
result to the case of perturbations.
3.1

Brief review of Reaction Network Theory

For
a = (ai )i∈S and b = (bi )i∈S , the notation ab will be shorthand for the formal monomial
Q vectors
bi
i∈S ai . We introduce some standard definitions.
Definition 3 (Reaction Network).
Fix a finite set S of species.
1. A reaction over S is a pair (y, y ′ ) such that y, y ′ ∈ ZS≥0 . It is usually written y → y ′ , with
reactant y and product y ′ .
2. A reaction network consists of a finite set S of species, and a finite set R of reactions.
3. A reaction network is reversible iff for every reaction y → y ′ ∈ R, the reaction y ′ → y ∈ R.
4. A reaction network is weakly reversible iff for every reaction y → y ′ ∈ R there exists a positive
integer n ∈ Z>0 and n reactions y1 → y2 , y2 → y3 , . . . , yn−1 → yn ∈ R with y1 = y ′ and yn = y.
5. The stoichiometric subspace H ⊆ RS is the subspace spanned by {y ′ − y | y → y ′ ∈ R}, and
H ⊥ is the orthogonal complement of H.
6. A siphon is a set T ⊆ S of species such that for all y → y ′ ∈ R, if there exists i ∈ T such that
yi′ > 0 then there exists j ∈ T such that yj > 0.
7. A siphon T ⊆ S is critical iff v ∈ H ⊥ ∩ RS≥0 with vi = 0 for all i ∈
/ T implies v = 0.
Definition 4. Fix a weakly reversible reaction network (S, R). The associated ideal I(S,R) ⊆ C[x]
′
where x = (xi )i∈S is the ideal generated by the binomials {xy − xy | y → y ′ ∈ R}. A reaction
network is prime iff its associated ideal is a prime ideal.
The following theorem follows from [13, Theorem 4.1, Theorem 5.2].
Theorem 2. A weakly reversible prime reaction network (S, R) has no critical siphons.
We now recall the mass-action equations which are widely employed for modeling cellular processes [29,26,28,30] in Biology.
Definition 5 (Mass Action System). A reaction system consists of a reaction network (S, R)
and a rate function k : R → R>0 . The mass-action equations for a reaction system are the
system of ordinary differential equations in concentration variables {xi (t) | i ∈ S}:
X
ẋ(t) =
(2)
ky→y′ x(t)y (y ′ − y)
y→y ′ ∈R

where x(t) represents the vector (xi (t))i∈S of concentrations at time t.
Note that ẋ(t) ∈ H, so affine translations of H are invariant under the dynamics of Equation 2.
We recall the well known notions of detailed balanced and complex balanced reaction system.
6

Definition 6. A reaction system (S, R, k) is
1. Detailed balanced iff it is reversible and there exists a point α ∈ RS>0 such that for every
y → y ′ ∈ R:
′
ky→y′ αy (y ′ − y) = ky′ →y αy (y − y ′ )
A point α ∈ RS>0 that satisfies the above condition is called a point of detailed balance.
2. Complex balanced iff there exists a point α ∈ RS>0 such that for every y ∈ ZS≥0 :
X

ky→y′ αy (y ′ − y) =

y→y ′ ∈R

X

′′

ky′′ →y αy (y − y ′′ )

y ′′ →y∈R

A point α ∈ RS>0 that satisfies the above condition is called a point of complex balance.
The following observations are well known and easy to verify.
– A complex balanced reaction system is always weakly reversible.
– If all rates ky→y′ = 1 and the network is weakly reversible then the reaction system is complex
balanced with point of complex balance (1, 1, . . . , 1) ∈ RS ; if the network is reversible then the
reaction system is also detailed balanced with point of detailed balance (1, 1, . . . , 1) ∈ RS .
– Every detailed balance point is also a complex balance point, but there are complex balanced
reversible networks that are not detailed balanced.
It is straightforward to check that every point of complex balance (respectively, detailed balance)
is a fixed point for Equation 2. The next theorem, which follows from [2, Theorem 2] and [15], states
that a converse also exists: if a reaction system is complex balanced (respectively, detailed balanced)
then every fixed point is a point of complex balance (detailed balance). Further there is a unique
fixed point in each affine translation of H, and if there are no critical siphons then the basin of
attraction for this fixed point is as large as possible, namely the intersection of the affine translation
of H with the nonnegative orthant.
Theorem 3 (Global Attractor Theorem for Complex Balanced Reaction Systems with
no critical siphons). Let (S, R, k) be a weakly reversible complex balanced reaction system with
no critical siphons and point of complex balance α. Fix a point u ∈ RS>0 . Then there exists a point
of complex balance β in (u + H) ∩ RS>0 such that for every trajectory x(t) with initial conditions
x(0) ∈ (u + H) ∩ RS≥0 , the limit limt→∞ x(t) exists and equals β. Further the function g(x) :=
Pn
i=1 xi log xi − xi − xi log αi is strictly decreasing along non-stationary trajectories and attains its
unique minimum value in (u + H) ∩ RS≥0 at β.
It is not completely trivial to show, but nevertheless true, that this theorem holds with weakly
reversible replaced by “reversible” and “complex balance” replaced by “detailed balance.” What is
to be shown is that the point of complex balance obtained in (u + H) ∩ RS≥0 by minimizing g(x)
is actually a point of detailed balance, and this follows from an examination of the form of the
d
g(x(t)) along trajectories x(t) to Equation 2.
derivative dt
3.2

A Perturbatively-Stable Global Attractor Theorem

Global attractor results usually assume that the reaction network is weakly reversible. We are going
to describe our scheme in the next section. Our scheme will employ reaction networks that are not
7

weakly reversible, yet we will prove global attractor results for them. The key idea we use is that
our reaction network can be broke into a reversible part, and an irreversible part. The reversible
part acts on, but evolves independent of, the irreversible part. So we get to use the global attractor
results “as is” on the reversible part. Further, as the reversible part approaches equilibrium, our
irreversible part behaves as a perturbation of a reversible detailed-balanced network. The closer the
reversible part gets to equilibrium, the smaller the perturbation of the irreversible part from the
dynamics of a certain reversible detailed-balanced network.
To make this proof idea work out, we will need a perturbative version of Theorem 3. The next
lemma shows that if the rates are perturbed slightly then, outside a small neighborhood of the
detailed balance point, the strict Lyapunov function g(x) from Theorem 3 continues to decrease
along non-stationary trajectories.
Lemma 1. Let (S, R, k) be a weakly reversible complex balanced reaction system with no critical
siphons and point of complex balance α. For every sufficiently small ǫ > 0 there exists δ > 0 such
d
that for all x′ outside the ǫ-neighborhood of α in (α + H) ∩ RS≥0 , the derivative dt
g(x(t))|t=0 < −δ,
′
where x(t) is a solution to the Mass-Action Equations 2 with x(0) = x .
Proof. Let Bǫ be the open ǫ ball around α in (α + H) ∩ RS≥0 , with ǫ small enough so that Bǫ
does not meet the boundary ∂RS≥0 . Consider the closed set S := (α + H) ∩ RS≥0 \ Bǫ . Define the
d
orbital derivative of g at x′ as Ok g(x′ ) := dt
g(x(t))|t=0 , where x(t) is a solution to the mass-action
′
equations 2 with x(0) = x . Define δ := inf x′ ∈S (−Ok g(x′ )). If δ ≤ 0 then since S is a closed set,
and Ok g is a continuous function, there exists a point x′ such that Ok g(x′ ) ≥ 0, which contradicts
Theorem 3.
We formalize the notion of perturbation using differential inclusions. Recall that differential
inclusions model uncertainty in dynamics in a nondeterministic way by generalizing the notion of
vector field. A differential inclusion maps every point to a subset of the tangent space at that point.
Definition 7. Let (S, R, k) be a reaction system and let δ > 0. The δ-perturbation of (S, R, k) is
S
the differential inclusion V : RS≥0 → 2R that at point x ∈ RS≥0 takes the value


 X

′
′
y ′
′
V (x) :=
ky→y
ky→y
.
′ x (y − y)
′ ∈ (ky→y ′ − δ, ky→y ′ + δ) for all y → y ∈ R


′
y→y ∈R

A trajectory of V is a tuple (I, x) where I ⊆ R is an interval and x : I → RS≥0 is a differentiable
function with ẋ(t) ∈ V (x(t)).

Theorem 4 (Perturbatively-Stable Global Attractor Theorem for Complex Balanced
Reaction Systems with no critical siphons). Let (S, R, k) be a weakly reversible complex balanced reaction system with no critical siphons. Fix a point u ∈ RS>0 . Then there exists a point of
complex balance β in (u + H) ∩ RS>0 such that:
1. For every sufficiently small ε > 0, there exists δ > 0 such that every trajectory of the form
(R≥0 , x) to the δ-perturbation of (S, R, k) with initial conditions x(0) ∈ (u + H) ∩ RS≥0 eventually
enters an ε-neighborhood of β and never leaves.
2. Consider a sequence δ1 > δ2 > · · · > 0 and a sequence 0 < t1 < t2 < . . . such that limi→∞ δi = 0
and limi→∞ ti = +∞, and a trajectory (R≥0 , x) with x(0) ∈ (u + H) ∩ RS≥0 such that ((ti , ∞), x)
is a trajectory of the δi -perturbation of (S, R, k). Then the limit lim x(t) = β.
t→∞

8

Proof (Proof sketch). 1. Fix ε > 0 such that the ε-ball Bε around β does not meet the boundary
∂RS≥0 . By Lemma 1, outside Bε , there exists δε > 0 such that the function Ok g < −δε . Since Ok g is
a continuous function of the specific rates k, a sufficiently small perturbation δ > 0 in the rates will
not change the sign of Ok g. Hence, outside Bǫ , the function g is strictly decreasing along trajectories
x(t) to Equation 2. It follows that eventually every trajectory must enter Bǫ .
2. Fix a sequence ε1 > ε2 > · · · > 0 with ε1 small enough so that the ε1 -ball around β does not
meet the boundary ∂RS≥0 and limi→∞ εi → 0. For each εi , there exists j such that δj is small enough
as per part (1) of the theorem. So every trajectory will eventually enter the ǫi neighborhood of β,
and never leave. Since this is true for every i and limi→∞ εi → 0, the result follows.

4

Main Result

The next definition makes precise our scheme, which takes a design matrix A to a reaction system
SMLE depending on A. The choice of this reaction system is not unique, but depends on two choices
of basis. We proceed in two stages. In the first stage, we construct the reaction system SMLD which
solves the problem of finding the maximum likelihood distribution. In the second stage, we add
reactions to solve for θ from the algebraic relations between the θ and X variables, obtaining SMLE .
Definition 8. Fix a design matrix A = (aij )m×n , a basis B for the free group Zn ∩ ker A, and a
maximal linearly-independent subset B ′ of the columns of A.
1. The reaction network RMLD (A, B) consists of n species X1 , X2 , . . . , Xn and for each b ∈ B, the
reversible reaction:
X
X
b j Xj ⇋
−bj Xj
j:bj >0

j:bj <0

2. The reaction system SMLD (A, B) consists of the reaction network RMLD (A, B) with an assignment of rate 1 to each reaction.
3. The reaction network RMLE (A, B, B ′ ) consists of m + n species θ1 , θ2 , . . . , θm , X1 , X2 , . . . , Xn ,
and in addition to the reactions in RMLD , the
Pfollowing reactions:
– For each column j ∈ B ′ of A, a reaction m
i=1 aij θi →
P0.
m
– For each column j ∈ B ′ of A, a reaction Xj → Xj + i=1 aij θi .
4. The reaction system SMLE (A, B, B ′ ) consists of the reaction network RMLE (A, B, B ′ ) with an
assignment of rate 1 to each reaction.
Note that by the rank-nullity theorem of linear algebra, the dimension of the kernel plus the rank
of the matrix equals the number of columns of the matrix. Hence counting the reversible reactions
as two irreversible reactions, our scheme yields a reaction system whose number of reactions is twice
the number of columns of A.
It is clear from the definition of SMLE that the reactions that come from RMLD are reversible and
evolve without being affected by the other reactions. Hence we first prove global convergence of the
reaction system SMLD to the maximum likelihood distribution. This part is fairly straightforward.
The key point is to verify that the reaction network RMLD has no critical siphons. In fact, we show
in the next lemma that RMLD is prime, which will imply “no critical siphons” by Theorem 2.
Lemma 2. Fix a design matrix A = (aij )m×n and a basis B for the free group Zn ∩ ker A. Then
the reaction network RMLD (A, B) is prime and SMLD (A, B) is detailed balanced. Consequently, the
reaction system SMLD (A, B) is globally asymptotically stable.
9

Proof. RMLD (A, B) is prime by [17, Corollary 2.15]. The idea is to look at the toric model pA as
a ring homomorphism C[x1 , x2 , . . . , xn ] → C[NA] with xj 7→ θaj . (Here NA is the affine semigroup
generated by the columns of A.) The kernel of this ring homomorphism is the associated ideal of
RMLD (A, B) by [17, Proposition 2.14], and the codomain is an integral domain, so the kernel must
be prime.
To verify that SMLD (A, B) is detailed balanced, note that the point (1, 1, . . . , 1) ∈ Rn is a point
of detailed balance since all rates are 1. Global asymptotic stability now follows from Theorem 2
and Theorem 3.
We can now obtain global convergence for SMLD .
Theorem 5 (The reaction system SMLD (A, B) computes the Maximum Likelihood Distribution). Fix a design matrix A = (aij )m×n , a basis B for the free group Zn ∩ ker A, and a
nonzero data vector u ∈ Zn≥0 . Let x(t) = (x1 (t), x2 (t), . . . , xn (t)) be a solution to the mass-action
differential equations for the reaction system SMLD (A, B) with initial conditions x(0) = u/|u|1 . Then
x(∞) := lim x(t) exists and equals the maximum likelihood distribution p̂(u).
t→∞

n
Proof. For the system SMLD
Pn (A, B), note that (x(0) + H) ∩ R>0 = PA (u/|u|1 ). By Theorem 3, x(∞)
exists, and the function i=1 xi log xi − xi −
Pxni log 1 attains its unique minimum in PA (u/|u|1 ) at
x(∞). Since the system is mass-conserving, i=1
Pnxi is constant on PA (u/|u|1 ), so this is equivalent
to the fact that Shannon entropy H(x) = − i=1 xi log xi is increasing, and attains its unique
maximum value in PA (u/|u|1 ) at x(∞). By Theorem 1, the point x(∞) must be the maximum
likelihood distribution p̂(u).

As the reversible reactions in SMLE approach closer and closer to equilibrium, we wish to absorb
the values of the X variables into reaction rates and pretend that the irreversible reactions are
reactions only in the θ variables. This has the advantage that we can treat this pretend reaction
system in the θ variables as a perturbation of a reversible, detailed balanced system. We can then
hope to employ Theorem 4 and conclude global convergence for these irreversible reactions, and
hence for SMLE .
One small technical point deserves mention. The pretend reaction system in the θ variables is not
a reaction system since the rates are not real numbers but functions of time. This will not trouble us.
We have already provisioned for this in Definition 7 by allowing perturbations of reaction systems
to be differential inclusions.
Theorem 6 (The reaction system SMLE (A, B, B ′ ) computes the Maximum Likelihood Estimator). Fix a design matrix A = (aij )m×n , a basis B for the free group Zn ∩ ker A, and a nonzero
data vector u ∈ Zn≥0 . Let x(t) = (x1 (t), x2 (t), . . . , xn (t), θ1 (t), θ2 (t), . . . , θm (t)) be a solution to the
mass-action differential equations for the reaction system SMLE (A, B, B ′ ) with initial conditions
x(0) = u/|u|1 and θ(0) = 0. Then x(∞) := limt→∞ x(t) exists and equals the maximum likelihood
distribution p̂(u), and θ(∞) := limt→∞ θ(t) exists and equals the maximum likelihood estimator θ̂(u).
Proof (Proof sketch). Fix u and let p̂ = p̂(u) and θ̂ = θ̂(u). Note that for the species X1 , X2 , . . . , Xn ,
the differential equations for SMLE (A, B) and SMLD (A, B, B ′ ) are identical, since these species
appear purely catalytically in the reactions that belong to RMLE (A, B, B ′ ) \ RMLD (A, B). Hence
x(∞) = p̂(u) follows from Theorem 5.
To see that θ(∞) = θ̂, let us first allow the X species to reach equilibrium, then treat the θ
system with replacing the X species by rate constants representing their values at equilibrium. The
10

system ΘMLE (A, B, B ′ , x(∞)) obtained in this way in only the θ species is a reaction system with
the reactions
P
– For each column j ∈ B ′ of A, a reaction m
aij θi → 0 of rate 1
i=1P
m
– For each column j ∈ B ′ of A, a reaction 0 → i=1 aij θi of rate xj (∞).
This is a reversible reaction system, and the maximum likelihood estimators θ̂ are precisely the points
of detailed balance for this system, where we are using the fact that B ′ was a maximal linearlyindependent set of the columns of A. In addition, this system has no siphonsPsince if species θi is
m
absent, and aij > 0 then θi will immediately be produced by the reaction 0 → i′ =1 ai′ j θi′ . (We are
assuming A has no 0 row. If A has a 0 row, we can ignore it anyway.) It follows from Theorem 3 that
this system is globally asymptotically stable, and every trajectory approaches a maximum likelihood
estimator θ̂.
Our actual system may be viewed as a perturbation of the system ΘMLE (A, B, B ′ , x(∞)). Consider any trajectory (x(t), θ(t)) to SMLE (A, B, B ′ ) starting at (u/|u|1 , 0). We are going to consider
the projected trajectory (R≥ , θ). We now show that it is possible to choose appropriate ti and δi so
that ((ti , ∞), θ(t)) is a trajectory of a δi -perturbation of ΘMLE (A, B, B ′ , x(∞)), for i = 1, 2, . . . .
Wait for a sufficiently large time t1 till x(t) is in a sufficiently small δ1 neighborhood of x(∞)
which it will never leave. After this time, we obtain a differential inclusion in the θ species with the
mass-action equations 2 for the reactions
P
– For each column j of A, a reaction m
aij θi → 0 of rate 1
i=1P
– For each column j of A, a reaction 0 → m
i=1 aij θi with time-varying rate lying in the interval
(xj (∞) − δ1 , xj (∞) + δ1 ).
Continuing in this way, we choose a decreasing sequence δ1 > δ2 > · · · > 0 with limi→∞ δi → 0,
and corresponding times t1 < t2 < t3 . . . with limi→∞ ti → ∞ such that after time ti , x(t) is
in a δi neighborhood of x(∞) which it will never leave. Then ((ti , ∞), θ(t)) is a trajectory of the
δi -perturbation of ΘMLE (A, B, B ′ , x(∞)). Hence θ(t) satisfies the conditions of Theorem 4. Hence
limt→∞ θ(t) = θ̂.

5

Related Work and Conclusions

The mathematical similarities of both log-linear statistics and reaction networks to toric geometry
have been pointed out before [9,17]. Craciun et al. [9] refer to the steady states of complex-balanced
reaction networks as Birch points “to highlight the parallels” with algebraic statistics. This paper
develops on these observations, and serves to flesh out this mathematical parallel into a scheme for
molecular computation.
Various building blocks for molecular computation that assume mass-action kinetics have been
proposed before. We briefly review some of these proposals.
In [19], Napp and Adams model molecular computation with mass-action kinetics, as we do here.
They propose a molecular scheme to implement message passing schemes in probabilistic graphical
models. The goal of their scheme is to convert a factor graph into a reaction network that encodes the
single-variable marginals of the joint distribution as steady state concentrations. In comparison, the
goal of our scheme is to do statistical inference and compute maximum likelihood estimators for loglinear models. Napp and Adams focus on the “forward model” task of how a given data-generating
process (a factor graph) can lead to observed data, whereas our focus is on the “backward model”
11

task of inference, going from the observed data to the data-generating process. Further our scheme
couples the deep role that MaxEnt algorithms play in Machine Learning with MaxEnt’s roots in
the Second Law of Thermodynamics whereas Napp and Adams are drawing their inspiration from
variable elimination implemented via message passing which has its roots in Boolean constraint
satisfaction problems.
Qian and Winfree [23,24] have proposed a DNA gate motif that can be composed to build large
circuits, and have experimentally demonstrated molecular computation of a Boolean circuit with
around 30 gates. In comparison, our scheme natively employs a continuous-time dynamical system
to do the computation, without a Boolean abstraction.
Taking a control theory point of view, Oishi and Klavins [20] have proposed a scheme for implementing linear input/output systems with reaction networks. Note that for a given matrix A, the
set of maximum likelihood distributions is usually not linear, but log-linear.
Daniel et al.[10] have demonstrated an in vivo implementation of feedback loops, exploiting
analogies with electronic circuits. It is possible that the success of their schemes is also related to
the toric nature of mass-action kinetics.
Buisman et al. [5] have proposed a reaction network scheme for computation of algebraic functions. The part of our scheme which reads out the maximum likelihood estimator from the maximum
likelihood distribution bears some similarity to their work.
One limitation of our present work is that the number of columns of the matrix A can become
very large, for example 2|V | for a graphical model with V nodes. Since the number of species and
number of reactions both depend on the number of columns of A, this can require an exponentially
large reaction network which may become impractical. One direction for future work is to extend our
scheme by specifying a reaction network that computes maximum likelihood for graphical models.
We have some freedom in our scheme in the choice of basis sets B and B ′ . In any chemical
implementation of this work, there might be opportunity for optimization in choice of basis.
Acknowledgements: I thank Nick S. Jones, Anne Shiu, Abhishek Behera, Ezra Miller, Thomas
Ouldridge, Gheorghe Craciun, and Bence Melykuti for useful discussions.

References
1. Alan Agresti. Categorical Data Analysis. John Wiley & Sons, 2013.
2. David Angeli, Patrick De Leenheer, and Eduardo D. Sontag. A Petri net approach to the study of
persistence in chemical reaction networks. Math. Biosci., 210(2):598–618, 2007.
3. Yaakov Benenson, Binyamin Gil, Uri Ben-Dor, Rivka Adar, and Ehud Shapiro. An autonomous molecular
computer for logical control of gene expression. Nature, 429(6990):423–429, 2004.
4. YMM Bishop, Stephen Feinberg, and Paul Holland. Discrete multivariant analysis. 1975.
5. HJ Buisman, Huub MM ten Eikelder, Peter AJ Hilbers, and Anthony ML Liekens. Computing algebraic
functions with biochemical reaction networks. Artificial life, 15(1):5–19, 2009.
6. Luca Cardelli. Strand Algebras for DNA Computing. Natural Computing, 10:407–428, 2011.
7. Luca Cardelli. Two-domain DNA strand displacement. Mathematical Structures in Computer Science,
23(02):247–271, 2013.
8. Ronald Christensen and R Christensen. Log-linear models and logistic regression, volume 168. Springer
New York, 1997.
9. Gheorghe Craciun, Alicia Dickenstein, Anne Shiu, and Bernd Sturmfels. Toric dynamical systems.
Journal of Symbolic Computation, 44(11):1551–1565, 2009. In Memoriam Karin Gatermann.

12

10. Ramiz Daniel, Jacob R Rubens, Rahul Sarpeshkar, and Timothy K Lu. Synthetic analog computation
in living cells. Nature, 497(7451):619–623, 2013.
11. John Doyle and Marie Csete. Rules of engagement. Nature, 446(7138):860–860, 2007.
12. Stephen E Fienberg, Alessandro Rinaldo, et al. Maximum likelihood estimation in log-linear models.
The Annals of Statistics, 40(2):996–1023, 2012.
13. Manoj Gopalkrishnan. Catalysis in reaction networks. Bulletin of Mathematical Biology, 73:2962–2982,
2011. 10.1007/s11538-011-9655-3.
14. Manoj Gopalkrishnan, Ezra Miller, and Anne Shiu. A geometric approach to the global attractor
conjecture. In preparation.
15. Friedrich J. M. Horn. The dynamics of open reaction systems. In Mathematical aspects of chemical and
biochemical problems and quantum chemistry, volume VIII of Proc. SIAM-AMS Sympos. Appl. Math.,
New York, 1974.
16. Steffen L Lauritzen. Graphical models. Oxford University Press, 1996.
17. Ezra Miller. Theory and applications of lattice point methods for binomial ideals. In Combinatorial
Aspects of Commutative Algebra and Algebraic Geometry, pages 99–154. Springer, 2011.
18. Ron Milo, Shai Shen-Orr, Shalev Itzkovitz, Nadav Kashtan, Dmitri Chklovskii, and Uri Alon. Network
motifs: simple building blocks of complex networks. Science, 298(5594):824–827, 2002.
19. Nils E Napp and Ryan P Adams. Message passing inference with chemical reaction networks. In
Advances in Neural Information Processing Systems, pages 2247–2255, 2013.
20. Kevin Oishi and Eric Klavins. Biomolecular implementation of linear I/O systems. Systems Biology,
IET, 5(4):252–260, 2011.
21. L. Pachter and B. Sturmfels. Algebraic Statistics for Computational Biology. Number v. 13 in Algebraic
Statistics for Computational Biology. Cambridge University Press, 2005.
22. Lulu Qian, David Soloveichik, and Erik Winfree. Efficient turing-universal computation with DNA
polymers. In DNA computing and molecular programming, pages 123–140. Springer, 2011.
23. Lulu Qian and Erik Winfree. A simple DNA gate motif for synthesizing large-scale circuits. J. R. Soc.
Interface, 2011.
24. Lulu Qian and Erik Winfree. Scaling up digital circuit computation with DNA strand displacement
cascades. Science, 332(6034):1196–1201, 2011.
25. Ehud Shapiro and Yaakov Benenson. Bringing DNA computers to life. Scientific American, 294(5):44–
51, 2006.
26. Guy Shinar and Martin Feinberg. Structural sources of robustness in biochemical reaction networks.
Science, 327(5971):1389–1391, 2010.
27. David Soloveichik, Georg Seelig, and Erik Winfree. DNA as a universal substrate for chemical kinetics.
Proceedings of the National Academy of Sciences, 107(12):5393–5398, 2010.
28. Eduardo D. Sontag. Structure and stability of certain chemical networks and applications to the kinetic
proofreading model of T-cell receptor signal transduction. IEEE Trans. Autom. Control, 46:1028–1047,
2001.
29. Matthew Thomson and Jeremy Gunawardena. Unlimited multistability in multisite phosphorylation
systems. Nature, 460(7252):274–277, 2009.
30. John J Tyson, Katherine C Chen, and Bela Novak. Sniffers, buzzers, toggles and blinkers: dynamics of
regulatory and signaling pathways in the cell. Current opinion in cell biology, 15(2):221–231, 2003.
31. Boyan Yordanov, Jongmin Kim, Rasmus L Petersen, Angelina Shudy, Vishwesh V Kulkarni, and Andrew
Phillips. Computational design of nucleic acid feedback control circuits. ACS synthetic biology, 3(8):600–
616, 2014.

13

