Computation with finite stochastic chemical reaction networks

Abstract

A highly desired part of the synthetic biology toolbox is an embedded chemical microcontroller, capable of autonomously following a logic program specified by a set of instructions, and interacting with its cellular environment. Strategies for incorporating logic in aqueous chemistry have focused primarily on implementing components, such as logic gates, that are composed into larger circuits, with each logic gate in the circuit corresponding to one or more molecular species. With this paradigm, designing and producing new molecular species is necessary to perform larger computations. An alternative approach begins by noticing that chemical systems on the small scale are fundamentally discrete and stochastic. In particular, the exact molecular counts of each molecular species present, is an intrinsically available form of information. This might appear to be a very weak form of information, perhaps quite difficult for computations to utilize. Indeed, it has been shown that error-free Turing universal computation is impossible in this setting. Nevertheless, we show a design of a chemical computer that achieves fast and reliable Turing-universal computation using molecular counts. Our scheme uses only a small number of different molecular species to do computation of arbitrary complexity. The total probability of error of the computation can be made arbitrarily small (but not zero) by adjusting the initial molecular counts of certain species. While physical implementations would be difficult, these results demonstrate that molecular counts can be a useful form of information for small molecular systems such as those operating within cellular environments.

A. Proof details

1.1 A.1 Clock analysis

The following three lemmas refer to the Markov chain in Fig. 3. We use p _i (t) to indicate the probability of being in state i at time t. CDF stands for cumulative distribution function.

Fig. 3

Continuous-time Markov chain for Lemmas 5–7. States i = 1,…,l indicate the identity of the currently present clock species C ₁,…,C _l . Transition to state 0 represents reaction dec ₂ for the RM simulation or the state transition initiation reaction of the CTM simulation

Lemma 5

Suppose the process starts in state l. Then $\forall t,\,p_1(t) \leq (1-p_0(t))\mu$ where $\mu=1/(1+\frac{r}{f}+(\frac{r}{f})^2+\cdots+(\frac{r}{f})^{l-1}).$

Proof

Consider the Markov chain restricted to states 1,…,l. We can prove that the invariance p _i+1(t)/p _i (t) ≥ r/f (for i = 1,…,l−1) is maintained at all times through the following argument. Letting ϕ _i (t) = p _i+1(t)/p _i (t), we can show dϕ _i (t)/dt ≥ 0 when ϕ _i (t) = r/f and $\forall i^{\prime}, \phi_{i^{\prime}}(t) \geq r/f,$ which implies that for no i can ϕ _i (t) fall below r/f if it starts above. This is done by showing that dp _i (t)/dt = p _i+1(t) f + p _i − 1(t) r−(r + f) p _i (t) ≤ 0 since ϕ _i (t) = r/f and ϕ _i − 1(t) ≥ r/f, and dp _i+1(t)/dt = p _i+2(t) f + p _i (t) r−(r + f) p _i+1(t) ≥ 0 since ϕ _i (t) = r/f and ϕ _i+1(t) ≥ r/f (the p _i − 1 or the p _i+2 terms are zero for the boundary cases).

Now p _i (t) = ϕ _i − 1(t) ϕ _i − 2(t) …ϕ₁(t) p ₁(t). Thus $\sum_i p_i(t)=1$ implies $p_1(t)=1/(1+\phi_1+\phi_2 \phi_1+\cdots+\phi_{l-1}\phi_{l-2}\cdots \phi_1)\leq 1/(1+{{r}\over{f}}+(\frac{r}{f})^2+\cdots+(\frac{r}{f})^{l-1}).$ This is a bound on the probability of being in state 1 given that we haven’t reached state 0 in the full chain of Fig. 3. Thus multiplying by 1−p ₀(t) gives us the desired result.$\square$

Lemma 6

Suppose for some μ we have $\forall t, \,p_1(t) \leq (1-p_0(t))\mu.$ Let T be a random variable describing the time until absorption at state 0. Then $\Pr[T < t]\leq 1-e^{-\lambda t}$ for λ = f μ (i.e. our CDF is bounded by the CDF for an exponential random variable with rate λ = f μ).

Proof

The result follows from the fact that dp ₀(t)/dt = p ₁(t) f ≤  (1−p ₀(t)) μ f. $\square$

Lemma 7

Starting at state l, the expected time to absorb at state 0 is $O((\frac{r}{f})^{l-1}/f)$ assuming sufficiently large r/f.

Proof

The expected number of transitions to reach state 0 starting in state i is $d_i = \frac{2 p q \left(\left(q/p\right)^l-\left(q/p\right)^{l-i}\right)}{(1-2 p)^2}-\frac{i}{q-p},$ where $p = \frac{f}{f+r}$ is the probability of transitioning to a state to the left and q = 1−p is the probability of transitioning to the state to the right. This expression is obtained by solving the recurrence relation d _i  = p d _i − 1 + q d _i+1 + 1 (0 > i > l) with boundary conditions d ₀ = 0, d _l  = d _l − 1 + 1. Thus $d_l < \frac{2 p q \left(q/p\right)^{l}}{(1-2 p)^2} = \frac{2 (r/f)^{l+1}}{(r/f - 1)^2}.$ This implies that for r/f larger than some constant, $d_l = O((\frac{r}{f})^{l-1}).$ Since the expected duration of any transition is no more than 1/f, the desired bound is obtained.

By the above lemmas, the time for the clock to “tick” can be effectively thought of as an exponential random variable with rate $\lambda = f/(1+\frac{r}{f} + (\frac{r}{f})^2 + \cdots + (\frac{r} {f})^{l-1}) = \Uptheta(\frac{f}{(r/f)^{l-1}}).$ Lemma 6 shows that the CDF of the tick is bounded by the CDF of this exponential random variable. Further, Lemma 7 shows that the expected time for the tick is bounded by (the order of) expected time of this exponential random variable. Note that Lemma 6 is true no matter how long the clock has already been running (a “memoryless” property). For our clock construction (Fig. 1b), we set λ by changing #A and #A* which define the forward and reverse rates f and r. Specifically, we have $\lambda=\Uptheta \left(\frac{k {\#A^{\ast}}^{l}}{v {\#A}^{l-1}}\right).$

1.2 A.2 Time/space-bounded RM simulation

Lemma 8

For the finite RM simulation, the probability of error per step is O((1/#A)^l − 1). Further, the expected time per step is bounded by O((#A)^l − 1 v/k).

Proof

Consider the point in time when the RM simulation enters a state in which it should decrement a non-empty register. If the time until dec ₂ occurs were an exponential random variable with rate λ then the probability of error per step would be bounded by λ/(k/v + λ). (We are making the worst case assumption that there is exactly one register molecule; otherwise, the error is even smaller.) The time until dec ₂ is not exponentially distributed, but by Sect. A.1, it can be bounded by an exponential random variable with rate $\lambda=O(\frac{k}{v{\#A}^{l-1}})$ (#A* = 1 for the RM construction). Note that the clock may have been running for a while since the last dec operation (while the RM performs inc operations for example); however, this amount of time is irrelevant by the memoryless property established in Sect. A.1. Thus the probability of error per step is bounded by λ/(k/v + λ) = O((1/#A)^l − 1). The expected time per RM step is bounded by the expected time for dec ₂ which is O((#A)^l − 1 v/k) by Sect. A.1. $\square$

The above lemma implies that we can use ${\#A}=\Uptheta((t/\delta)^{1/(l-1)})$ resulting in the expected time for the whole computation of $O(\frac{vt^2}{k \delta})$ and the total probability of error being bounded by δ.

1.3 A.3 Time/space-bounded CTM simulation

In the following lemmas, we say a reaction completely finishes when it happens enough times that one of the reactants is used up.

Lemma 9

Starting with $\Uptheta(m)$ molecules of X and $\Uptheta(m)$ molecules of Y, the expected time for the reaction X + Y →Y to completely finish is $O(\frac{v}{k m} \hbox{log}\,{m}).$ The variance of the completion time is $O((\frac{v}{k m})^2).$

Proof

When there are q molecules of X remaining, the waiting time until next reaction is an exponential random variable with rate $\Uptheta(kqm/v)$ and therefore mean $O(\frac{v}{kqm}).$ Each waiting time is independent. Thus the total expected time is $\sum_{q=1}^{\Uptheta(m)}O(\frac{v}{kqm}) = O(\frac{v} {km}\hbox{log}\,{m}).$ ⁸ The variance of each waiting time is $O((\frac{v}{kqm})^2).$ Thus the total variance is $\sum_{q=1}^{\Uptheta(m)}O((\frac{v}{kqm})^2) = O((\frac{v}{k m})^2).$ $\square$

Lemma 10

Starting with $\Uptheta(m)$ molecules of X and $\Uptheta(m)$ molecules of Y such that ${\Updelta}={\#Y}-{\#X}=\Upomega(m)$ the expected time for the reaction $X + Y \rightarrow \emptyset$ to completely finish is $O(\frac{v}{k m} \hbox{log}\,{m}).$ The variance of the completion time is $O((\frac{v}{k m})^2).$

Proof

This case can be proven by reducing to Lemma 9 with initial amounts $\#Y^{\prime} =\Updelta$ and #X′ = #X. $\square$

Lemma 11

Starting with $\Uptheta(m)$ molecules of X and 1 molecule of Y, the expected time for the reaction X + Y →2Y to completely finish is $O(\frac{v}{k m} \hbox{log}\,{m}).$ The variance of the completion time is $O((\frac{v}{k m})^2).$

Proof

Consider splitting the process into two halves, with the first part bringing the amount of X to half its initial value and the second part using up the remainder. The time-reverse of the first part, as well as the second part, can both be bounded by processes covered by Lemma 9. (Assume that #X is fixed at its minimal value for part one, and assume #Y is fixed at its minimal value for part two. The variance can only decrease.)$\square$

Lemma 12

Some $\lambda = \Uptheta(\frac{k \varepsilon^{3/2} 3^{s_{ct}}}{v s_{ct}})$ attains error at most ɛ per microstep of the CTM simulation.

Proof

Using the above lemmas with $m = 3^{s_{ct}-1},$ by Chebyshev’s inequality,⁹ with probability at least 1−ɛ/2 all reactions finish before some time $t_{f} = \Uptheta(\frac{v}{km} (\hbox{log}(m) + 1/\sqrt{\varepsilon})) = O\left(\frac{v\,{\hbox{log}}\,{m}}{k m \varepsilon^{1/2}}\right).$ Now we set λ such that the probability that the clock ticks before time t _f is smaller than ɛ/2 (for a total probability of error ɛ). Since the time until the clock ticks is bounded by the CDF of an exponential random variable with rate λ (Sect. A.1), it is enough that $\lambda < \frac{\varepsilon}{2 t_f}$ and so we can choose some $\lambda = \Uptheta \left(\frac{\varepsilon^{3/2} k m}{v\,{\hbox{log}}\,{m}}\right).$

Lemma 13

Any TM with a two-way infinite tape using at most s _tm space and t _tm time can be converted to a CTM using s _ct  = 2s _tm space and t _ct  = O(t _tm s _tm ) time. If $\Updelta$ extra bits of padding on the CTM tape is used, then $t_{ct} =O(t_{tm} (s_{tm}+\Updelta))$ time is required.

Proof

(sketch, see Neary and Woods 2005) Two bits of the CTM are used to represent a bit of the TM tape. The extra bit is used to store a TM head position marker. To move in the direction corresponding to moving the CTM head clockwise (the easy direction) is trivial. To move in the opposite direction, we use the temporary marker to record the current head position and then move each tape symbol clockwise by one position. Thus, a single TM operation in the worst case corresponds to O(s) CTM operations. $\square$

In order to simulate t _tm steps of a TM that uses s _tm bits of space on a CTM using $\Updelta$ bits of padding requires $t_{ct} = O(t_{tm} (s_{tm}+\Updelta))$ CTM steps and a circular tape of size $s_{ct} =2s_{tm}+\Updelta$ (Lemma 13). Recall that in our CTM simulation, there are four microsteps corresponding to a single CTM operation, which is asymptotically still O(t _ct ). Thus, in order for the total error to be at most δ, we need the error per CTM microstep to be $\varepsilon = O(\frac{\delta}{t_{tm} (s_{tm}+\Updelta)}).$ Setting the parameters of the clock module (#A, #A*) to attain the largest λ satisfying Lemma 12, the expected time per microstep is $O\left(\frac{v s_{ct}}{k 3^{s_{ct}} \varepsilon^{3/2}}) = O(\frac{v (s_{tm}+\Updelta)^{5/2} t_{tm}^{3/2}}{k 3^{2s_{tm}+\Updelta} \delta^{3/2}}\right).$ This can be done, for example, by setting ${\#A^{\ast}}^l = \Uptheta \left(\frac{3^{s_{ct}}}{s_{ct}}\right)$ and ${\#A}^{l-1} = \Uptheta \left(\frac{1}{\varepsilon^{3/2}}\right).$ Since there are total $O(t_{tm}(s_{tm}+\Updelta))$ CTM microsteps, the total expected time is $O\left(\frac{v (s_{tm}+\Updelta)^{7/2} t_{tm}^{5/2}}{k 3^{(2s_{tm}+\Updelta)} \delta^{3/2}}\right).$

How large is the total molecular count? If we keep δ constant while increasing the complexity of the computation being performed, and setting #A* and #A as suggested above, we have that the total molecular count is $\Uptheta(m+\#A)$ where $m = 3^{2 s_{tm}+\Updelta}.$ Now m increases at least exponentially with $s_{tm}+\Updelta,$ while #A increases at most polynomially. Further, m increases at least quadratically with t _tm (for any reasonable algorithm $2^{s_{tm}} \geq t_{tm}$) while #A increases at most as a polynomial of degree $(3/2)\frac{1}{l-1} < 2.$ Thus m will dominate.

1.4 A.4 Unbounded RM simulation

After i dec ₂ steps, we have #A = i ₀ + i where i ₀ is the initial number of A’s. The error probability for the next step is O(1/#A ²) = O(1/(i ₀ + i)²) by Lemma 8 when l = 3. The total probability of error over an unbounded number of steps is $O(\sum_{i=0}^\infty 1/(i_0+i)^2).$ To make sure this is smaller than δ we start out with $i_0 =\Uptheta(1/\delta)$ molecules of A.¹⁰

Now what is the total expected time for t steps? By Lemma 8 the expected time for the next step after i dec ₂ steps is $O(\#A^2 v/k)=O((i_0+i)^2 v/k).$ Since each step at most increases the total molecular count by 1, after t total steps v is not larger than O(i ₀ + t + s ₀), where s ₀ is the sum of the initial values of all the registers. Thus the expected time for the tth step is bounded by $O((i_0+i)^2 (i_0+t+s_0)/k) =O((1/{\delta}+t)^2 (1/{\delta}+t+ s_0)/k)$ and so the expected total time for t steps is $O(t (1/{\delta+}t)^2 (1/{\delta}+t+s_0)/k).$

1.5 A.5 Unbounded CTM simulation

We want to follow a similar strategy as in the RM simulation (Sect. A.4) and want the error probability on the ith CTM step to be bounded by ${\varepsilon}=1/(\Uptheta(1/\delta)+i)^2$ such that the total error probability after arbitrarily many steps is bounded by δ. By Lemma 12, we can attain per step error probability (taking the union bound over the 4 microsteps in a step) bounded by this ɛ when we choose a small enough $\lambda = \Uptheta(\frac{k \varepsilon^{3/2} 3^{s_{ct}}}{v s_{ct}}) = \Uptheta(\frac{k 3^{s_{ct}}}{v (1/\delta + i)^3 s_{ct}}),$ where s _ct is the current CTM tape size. Recall that λ is set by #A and #A* such that $\lambda = \Uptheta(\frac{k {\#A^{\ast}}^l}{v \#A^{l-1}})$ (Sect. A.1). It is not hard to see that we can achieve the desired λ using clock Markov chain length l = 5, and appropriate $\#A = \Uptheta((i_0+i) 3^{s_{ct}})$ and $\#A^{\ast} = \Uptheta(3^{s_{ct}}),$ for appropriate $i_0 = \Uptheta(1/\delta + {s_{ct}}_0),$ where ${s_{ct}}_0$ is the initial size of the tape. These values of #A and #A* can be attained if the SCRN triples the amount of A and A* whenever extending the tape and increases #A by an appropriate amount $\Uptheta(3^{s_{ct}})$ on every step.

How fast is the simulation with these parameters? From Sect. A.1 we know that the expected time per microstep is $O(1/\lambda) = O(\frac{v (1/\delta+{s_{ct}}_0+i)^4}{k 3^{s_{ct}}}).$ Since the total molecular count is asymptotically $O(\#A) = O((1/\delta + {s_{ct}}_0+i) 3^{s_{ct}}),$ this expected time is $O((1/\delta + {s_{ct}}_0+ i)^5/k).$ However, unlike in the bounded time/space simulations and the unbounded RM simulation, this expected time is conditional on all the previous microsteps being correct because if a microstep is incorrect, A and A* may increase by an incorrect amount (for example reactions tripling #A akin to $A \rightarrow A^\dagger$ and $A^\dagger \rightarrow 3A$ can drive A arbitrarily high if the catalyst state species for both reactions are erroneously present simultaneously). Nonetheless, the expected duration of a microstep conditional on the entire simulation being correct is at most a factor of 1/(1−δ) larger than this.¹¹ Since we can assume δ will always be bounded above by a constant less than one, the expected duration of a microstep conditional on the entire simulation being correct is still $O((1/{\delta}+{s_{ct}}_0+i)^5/k).$ By Lemma 13, this yields total expected time to simulate t _tm steps of a TM using at most s _tm space and with initial input of size ${s_{tm}}_0$ is $O((1/\delta + {s_{tm}}_0 + t_{tm} s_{tm})^5 t_{tm} s_{tm}/k)$ assuming the entire simulation is correct.

1.6 A.6 Decidability of reachability

We reduce the reachability question in SCRNs to the reachability question in Vector Addition Systems (VAS), a model of asynchronous parallel processes developed by Karp and Miller (1969). In the VAS model, we consider walks through a p dimensional integer lattice, where each step must be one of a finite set of vectors in ${\mathbb{N}}^p,$ and each point in the walk must have no negative coordinates. It is known that the following reachability question is decidable: given points x and y, is there a walk that reaches some point ${\mathbf y}^{\prime} \geq {\mathbf y}$ from x (Karp and Miller 1969). The correspondence between VASs and SCRNs is straightforward (Cook 2005). First consider chemical reactions in which no species occurs both as a reactant and as a product (i.e. reactions that have no catalysts). When such a reaction $\alpha=\langle {\mathbf l},{\mathbf r},k \rangle$ occurs, the state of the SCRN changes by addition of the vector $-{\mathbf l}+{\mathbf r}.$ Thus the trajectory of states is a walk through ${\mathbb{N}}^p$ wherein each step is any of a finite number of reactions, subject to the constraint requiring that the number of molecules of each species remain non-negative. Karp and Miller’s decidability results for VASs then directly imply that our reachability question of whether we ever enter a state greater than or equal to some target state is decidable for catalyst-free SCRNs. The restriction to catalyst-free reactions is easily lifted: each catalytic reaction can be replaced by two new reactions involving a new molecular species after which all reachability questions (not involving the new species) are identical for the catalyst-free and the catalyst-containing networks.