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).2EE 5393, Winter ’17 4(b) Now, instead of “outcomes”, let’s analyze probabilities in a more generalsense.Consider again the reactions:R1 : 2X1 + X2 → 4X3 k1 = 1R2 : X1 + 2X3 → 3X2 k2 = 2R3 : X2 + X3 → 2X1 k3 = 3Let the state be S = [x1, x2, x3], where x1, x2 and x3 denote the numbersof 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 16.• it is in state [2, 6, 1] with probability 13.• it is in state [5, 2, 2] with probability 12.So, considering the first type, X1, its discrete probability distribution afterone 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 differentquantities of X1. (Of course, some of these states may have the samequantity 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 5ProblemFor the set of reactions above, beginning from the state [6, 6, 6] compute(or estimate) the mean and variance for the probability distributions ofX1, X2 and X3 – each separately – after 5 steps. 2EE 5393, Winter ’17 62. Synthesizing Chemical Reaction NetworksThe task of synthesizing a set of chemical reactions to compute a desired func-tion is conceptually open-ended. Like most engineering problems, there areoften many possible solutions.In class, we saw a chemical reaction network that performs multiplication. Whatfollows are some other examples of chemical reaction networks that computesimple functions. In describing the functions that the modules implement, weadd subscripts to the quantities of molecular types to denote when these quan-tities exist: zero indicates that this is the initial quantity, whereas infinityindicates that it is the quantity after the module has finished.Exponentiation:Y∞ = 2X0This module consumes molecules of an input type one at a time, doubling thequantity of an output type for each. Its behavior is described by the followingpseudocode:1 ForEach x {2 Y = 2 * Y;3 }.The reactions are:xslow−→ aa + yfaster−→ a + 2y′afast−→ ∅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 ofdoubling 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 followingpseudocode:EE 5393, Winter ’17 71 While Not(X==1) {2 X = X/2;3 Y = Y+1;4 }.The reactions are:bslow−→ a + ba + 2xfaster−→ c + x′ + a2cfaster−→ cafast−→ ∅x′medium−→ xcmedium−→ y.Initially, B is a small but non-zero quantity and all other quantities (except X)are zero.ProblemProduce a chemical reaction network that computesZ∞ = X0 log2 Y0Demonstrate that your solution works (either mathematically, or by simulatingit continuously or stochastically).ProblemProduce a chemical reaction network that computesY∞ = 2log2 X0 (1)(2)(No points for noticing that Y∞ = X0. Your network must compute this asshown!) Demonstrate that your solution works (either mathematically, or bysimulating it continuously or stochastically).EE 5393, Winter ’17 83. Multiplication(no collaboration)Consider the following representation of real numbers. A real value x between0 and 1 is represented as x1x1+x2, where x1 and x2 are positive integers.Construct a set of chemical reactions that multiplies two real numbers a andb represented this way, producing a resulting number c, also represented thisway. 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 simulatingit continuously or stochastically).EE 5393, Winter ’17 94. Iterative Computation with Molecular ReactionsIn the last couple of problems, you implemented some simple functions withmolecular reactions. In this homework, you’ll implement two full-fledged algo-rithms.(a) Euclid’s AlgorithmEuclid’s algorithm is an efficient method for computing the greatest com-mon divisor (GCD) of two integers, also known as the greatest commonfactor (GCF) or highest common factor (HCF). It is named after the Greekmathematician 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 positiveintegers and forms a new pair that consists of the smaller number and thedifference between the larger and smaller numbers. The process repeatsuntil the numbers are equal. That number then is the greatest commondivisor 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 ProcedureThe Collatz conjecture is a famous open problem in mathematics, proposedby Lothar Collatz in 1937. Consider the following iterative procedure. Forany 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 procedurealways terminates with x = 1. For instance, starting with x = 5, onefollows the sequence 16, 8, 4, 2 and 1. Starting from x = 27, one followsthe 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 amonetary 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 implementsthe procedure. The input to the system is a quantity of a type X; thesystem should iterate through the Collatz sequence until it hits one.EE 5393, Winter ’17 10Demonstrate that your code works for x = 27.
Reviews
There are no reviews yet.