EE 5393 UMN
Circuits, Computation, and Biology Winter 2017
⊕ ⊕
Homework # 1
Due Thu. Feb 22, 2017
1. Analyzing Chemical Reaction Networks
The theory of reaction kinetics underpins our understanding of biological and
chemical systems. It is a simple and elegant formalism: chemical reactions
define rules according to which reactants form products; each rule fires at a
rate that is proportional to the quantities of the corresponding reactants that
are present. On the computational front, there has been a wealth of research
into efficient methods for simulating chemical reactions, ranging from ordinary
differential equations (ODEs) to stochastic simulation. On the mathematical
front, entirely new branches of theory have been developed to characterize the
dynamics of chemical reaction networks.
Most of this work is from the vantage point of analysis: a set of chemical
reaction exists, designed by nature and perhaps modified by human engineers;
the objective is to understand and characterize its behavior. Comparatively
little work has been done at a conceptual level in tackling the inverse problem
of synthesis: how can one design a set of chemical reactions that implement
specific behavior?
This homework will consider the computational power of chemical reactions
from both a deductive and an inductive point of view.
(a) A molecular system consists of a set of chemical reactions, each specifying
a rule for how types of molecules combine. For instance,
X1 + X2
k−→ X3,
specifies that one molecule of X1 combines with one molecule of X2 to
produce one molecule of X3. The rate at which the reaction fires is pro-
portional to (1) the concentrations of the participating molecular types;
and (2) a rate constant k. (This value is not constant at all; rather it
is dependent on factors such as temperature and volume; however, it is
independent of molecular quantities, and so called a “constant.”)
EE 5393, Winter ’17 2
Given a set of such reactions, we can model the behavior of the system in
two ways:
i. In a continuous sense, in terms of molecular concentrations, with
differential equations;
ii. In a discrete sense, in terms of molecular quantities, through proba-
bilistic discrete-event simulation.
Consider the reactions:
R1 : 2X1 + X2 → 4X3 k1 = 1
R2 : X1 + 2X3 → 3X2 k2 = 2
R3 : X2 + X3 → 2X1 k3 = 3
For a continuous model, let x1, x2 and x3 denote the concentrations of
X1, X2, and X3, respectively. (Recall that concentration is number of
molecules per unit volume.) The behavior of the system is described by
the following set of differential equations:
dx1
dt
= −x21x2 − 2x1x
2
3 + 6x2x3
dx2
dt
= −x21x2 + 6x1x
2
3 − 3x2x3
dx3
dt
= 4x21x2 − 2x1x
2
3 − 3x2x3
For the discrete model, let the state be S = [x1, x2, x3], where x1, x2 and
x3 denote the numbers of molecules of types X1, X2, and X3, respectively.
(Here we use actual integer quantities, not concentrations.) The firing
probabilities for R1, R2 and R3 are computed as follows:
p1(x1, x2, x3) =
1
2
x1(x1 − 1)x2
1
2
x1(x1 − 1)x2 + x1x3(x3 − 1) + 3x2x3
,
p2(x1, x2, x3) =
x1x3(x3 − 1)
1
2
x1(x1 − 1)x2 + x1x3(x3 − 1) + 3x2x3
,
p3(x1, x2, x3) =
3x2x3
1
2
x1(x1 − 1)x2 + x1x3(x3 − 1) + 3x2x3
.
Suppose that S = [3, 3, 3]. Then the firing probabilities for R1, R2 and R3
EE 5393, Winter ’17 3
are
p1(3, 3, 3) =
9
9 + 18 + 27
=
1
6
,
p2(3, 3, 3) =
18
9 + 18 + 27
=
1
3
,
p3(3, 3, 3) =
27
9 + 18 + 27
=
1
2
,
respectively.
N.B. In the continuous model, the rate of change of type is proportional
to xn where x is the concentration of a reaction and n is the coefficient. In
the discrete model, the probability is proportional to
(
x
n
)
. This is a subtle
difference. See the paper by Gillespie for an explanation.
Problem
Suppose that we define the following “outcomes”:
• C1: states S = [x1, x2, x3] with x1 > 7.
• C2: states S = [x1, x2, x3] with x2 ≥ 8.
• C3: states S = [x1, x2, x3] with x3 < 3.
Beginning from the state S = [5, 5, 5], compute (or estimate) Pr(C1),
Pr(C2), and Pr(C3).
2
EE 5393, Winter ’17 4
(b) Now, instead of “outcomes”, let’s analyze probabilities in a more general
sense.
Consider again the reactions:
R1 : 2X1 + X2 → 4X3 k1 = 1
R2 : X1 + 2X3 → 3X2 k2 = 2
R3 : X2 + X3 → 2X1 k3 = 3
Let the state be S = [x1, x2, x3], where x1, x2 and x3 denote the numbers
of molecules of types X1, X2, and X3, respectively.
Suppose that systems begins in the state S = [3, 3, 3] (with probability 1).
After one step:
• it is in state [1, 2, 7] with probability 1
6
.
• it is in state [2, 6, 1] with probability 1
3
.
• it is in state [5, 2, 2] with probability 1
2
.
So, considering the first type, X1, its discrete probability distribution after
one step is
• Pr[X1 = 1] = 16 ,
• Pr[X1 = 2] = 13 ,
• Pr[X1 = 5] = 12 ,
After many steps, the system can be in many different states, with different
quantities of X1. (Of course, some of these states may have the same
quantity of X1.) The probability distribution may look something like:
• Pr[X1 = 0] = 0.012,
• Pr[X1 = 1] = 0.025,
• Pr[X1 = 2] = 0.036,
• Pr[X1 = 3] = 0.061,
• Pr[X1 = 4] = 0.12,
• Pr[X1 = 5] = 0.19,
• Pr[X1 = 6] = 0.24,
• Pr[X1 = 7] = 0.19,
• Pr[X1 = 8] = 0.116,
• Pr[X1 = 9] = 0.010.
(Note that this is not a real calculation.)
EE 5393, Winter ’17 5
Problem
For the set of reactions above, beginning from the state [6, 6, 6] compute
(or estimate) the mean and variance for the probability distributions of
X1, X2 and X3 – each separately – after 5 steps. 2
EE 5393, Winter ’17 6
2. Synthesizing Chemical Reaction Networks
The task of synthesizing a set of chemical reactions to compute a desired func-
tion is conceptually open-ended. Like most engineering problems, there are
often many possible solutions.
In class, we saw a chemical reaction network that performs multiplication. What
follows are some other examples of chemical reaction networks that compute
simple functions. In describing the functions that the modules implement, we
add subscripts to the quantities of molecular types to denote when these quan-
tities exist: zero indicates that this is the initial quantity, whereas infinity
indicates that it is the quantity after the module has finished.
Exponentiation:
Y∞ = 2
X0
This module consumes molecules of an input type one at a time, doubling the
quantity of an output type for each. Its behavior is described by the following
pseudocode:
1 ForEach x {
2 Y = 2 * Y;
3 }.
The reactions are:
x
slow−→ a
a + y
faster−→ a + 2y′
a
fast−→ ∅
y′
medium−→ y.
Initially, Y is one and all other quantities (except X) are zero.
Logarithm:
Y∞ = log2(X0)
This module is similar to the exponentiation module, except that instead of
doubling the output, the input is forced to halve itself; each time it does so,
the output is incremented by one. Its behavior is described by the following
pseudocode:
EE 5393, Winter ’17 7
1 While Not(X==1) {
2 X = X/2;
3 Y = Y+1;
4 }.
The reactions are:
b
slow−→ a + b
a + 2x
faster−→ c + x′ + a
2c
faster−→ c
a
fast−→ ∅
x′
medium−→ x
c
medium−→ y.
Initially, B is a small but non-zero quantity and all other quantities (except X)
are zero.
Problem
Produce a chemical reaction network that computes
Z∞ = X0 log2 Y0
Demonstrate that your solution works (either mathematically, or by simulating
it continuously or stochastically).
Problem
Produce a chemical reaction network that computes
Y∞ = 2
log2 X0 (1)
(2)
(No points for noticing that Y∞ = X0. Your network must compute this as
shown!) Demonstrate that your solution works (either mathematically, or by
simulating it continuously or stochastically).
EE 5393, Winter ’17 8
3. Multiplication
(no collaboration)
Consider the following representation of real numbers. A real value x between
0 and 1 is represented as x1
x1+x2
, where x1 and x2 are positive integers.
Construct a set of chemical reactions that multiplies two real numbers a and
b represented this way, producing a resulting number c, also represented this
way. Let the quantities or concentrations of molecular types A1 and A2 rep-
resent a, those of B1 and B2 represent b, and those of C1 and C2 represent c.
Demonstrate that your solution works (either mathematically, or by simulating
it continuously or stochastically).
EE 5393, Winter ’17 9
4. Iterative Computation with Molecular Reactions
In the last couple of problems, you implemented some simple functions with
molecular reactions. In this homework, you’ll implement two full-fledged algo-
rithms.
(a) Euclid’s Algorithm
Euclid’s algorithm is an efficient method for computing the greatest com-
mon divisor (GCD) of two integers, also known as the greatest common
factor (GCF) or highest common factor (HCF). It is named after the Greek
mathematician Euclid, who described it in Books VII and X of his Ele-
ments.
In its simplest form, Euclid’s algorithm starts with a pair of positive
integers and forms a new pair that consists of the smaller number and the
difference between the larger and smaller numbers. The process repeats
until the numbers are equal. That number then is the greatest common
divisor of the original pair.
Design a set of molecular reactions that implements the procedure.
Demonstrate that your code works for x = 66 and y = 30.
(b) Collatz Procedure
The Collatz conjecture is a famous open problem in mathematics, proposed
by Lothar Collatz in 1937. Consider the following iterative procedure. For
any positive integer x,
• if x = 1 stop;
• else if x is odd, let x = 3x + 1;
• else let x = x/2.
The conjecture is that, starting with any positive integer x, the procedure
always terminates with x = 1. For instance, starting with x = 5, one
follows the sequence 16, 8, 4, 2 and 1. Starting from x = 27, one follows
the sequence 82, 41, 124, 62, 31, 94, 47, 142, 71, 214, 107, . . . (keep going,
you’ll see that you eventually hit one.)
Proving this is evidently difficult. Paul Erdös said about the con-
jecture: “Mathematics is not yet ready for such problems”. He offered a
monetary reward of $500 for its solution.
You are not asked to prove the Collatz conjecture on this homework.
Rather you are asked to design a set of molecular reactions that implements
the procedure. The input to the system is a quantity of a type X; the
system should iterate through the Collatz sequence until it hits one.
EE 5393, Winter ’17 10
Demonstrate that your code works for x = 27.