[SOLVED] algorithm EE 5393 UMN

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, lets 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 + 2yafast ymedium 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 xmedium 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, youll implement two full-fledged algo-rithms.(a) Euclids AlgorithmEuclids 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, Euclids 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,youll see that you eventually hit one.)Proving this is evidently difficult. Paul Erdos 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.