Rate-independent Computation in Continuous Chemical Reaction Networks

1.1 Summary of Main Results

Our first contribution is to define a reachability relation that captures the broadest reasonable notion of “what could happen” and is of independent interest. Although concentration trajectories of mass-action kinetics (and other standard rate laws) follow smooth curves, we base the reachability relation on taking simple-to-analyze straight-line paths. Theorem 2.27 shows that this notion of reachability is exactly equivalent to satisfying the three intuitive properties described above, (1) nonnegative concentrations, (2) reactions require their reactants present, and (3) respecting causal relationships between the production of species. Thus, our reachability relation has all reasonable rate laws as special cases; i.e., if any of them can reach a state, then so can our reachability relation.

The reachability relation allows us to formally define stable computation in Definition 3.2, analogously to similar definitions of probability 1 computation in discrete systems [6, 18]. Stable computation allows us to delineate when a function cannot be computed rate-independently. Roughly, unless the system is stably computing, then an adversary can always push it “far” away from the correct output (theorem 3.3), precluding the system from being reasonably rate-independent.

For the positive direction, the CRN should converge to the correct output no matter the reaction rates. While a CRN that does not stably compute is not rate-independent (which is sufficient for negative results), the positive direction does not directly follow from stable computation for continuous systems. Indeed, we show examples of CRNs that stably compute a function by our definition, yet under standard mass-action kinetics fail to converge to the correct output; see Section 4. Instead, we capture a very strong notion of “convergence despite perturbations” in fair computation (Definition 4.3), based on generalized rate laws (so-called fair rate schedules, Definition 4.1). A CRN that fairly computes converges to the correct output as time \(t \rightarrow \infty\) under any trajectory satisfying the three intuitive conditions above, plus an additional requirement that reactions do occur when applicable. (In particular, mass-action (Corollary 4.11) satisfies these conditions, but the range of rate laws satisfying the conditions is much broader.) Luckily, stable computation and fair computation can be tightly connected, and we show that for a special class of CRNs we call feedforward (Definition 4.5) the two definitions coincide. In other words, a feedforward CRN stably computes a function if and only if it fairly computes the function (Lemmas 4.9 and 4.4). We show that all functions stably computable by CRNs are computable by feedforward CRNs (Lemmas 5.16 and 5.15), implying that the class of functions computable by CRNs under either definition—stable computation or fair computation—is identical. In other words, we can freely work with the simpler definition of stable computation, knowing that we are actually reasoning about a very general notion of rate-independence.

The above line of reasoning leads us to conclude that exactly the functions that are positive-continuous, piecewise linear (direct) or continuous, piecewise linear (dual-rail) can be rate-independently computed (Theorems 5.9 and 5.10). Positive-continuous means that the only discontinuities occur on a “face” of \(\mathbb {R}_{\ge 0}^k\)—i.e., the function may discontinuously jump only at a point where some input goes from 0 to positive. We already saw a simple example of a continuous, piecewise linear function (max function, 1(a,b,c)). 1(d,e) shows a more complex example and a CRN that computes it. 1(f,g) shows a discontinuous but positive-continuous function and a CRN that computes it. Although our work shows that the computational power of rate-independent CRNs is limited, the power of the computable class of functions should not be underestimated. For example, allowing a fixed non-zero initial concentration of non-input species (see Section 6.2), such CRNs are equivalent to ReLU neural networks—arguably the most widely used type of neural networks in machine learning [52].

1.2 Chemical Motivation

Traditionally, CRNs have been used as a descriptive language to analyze naturally occurring chemical reactions, as well as various other systems with a large number of interacting components such as gene regulatory networks and animal populations. However, CRNs also constitute a natural choice of programming language for engineering artificial systems. For example, nucleic-acid networks can be rationally designed to implement arbitrary chemical reaction networks [14, 20, 50, 51]. Thus, since in principle any CRN can be physically built, hypothetical CRNs with interesting behaviors are becoming of more than theoretical interest. One day artificial CRNs may underlie embedded controllers for biochemical, nanotechnological, or medical applications, where environments are inherently incompatible with traditional electronic controllers. However, to effectively program chemistry, we must understand the computational power at our disposal. In turn, the computer science approach to CRNs is also beginning to generate novel insights regarding natural cellular regulatory networks [15].

At the fine-grained level of detail, chemistry is discrete and stochastic. This level is typically modeled by discrete CRNs, where the state is a vector of nonnegative integers representing the counts of each species in the given reaction vessel, and reactions are modeled by a Markov jump process [29]. The continuous model is governed by a system of mass-action ordinary differential equations, which can be derived as a limiting case of the discrete model when volume and counts are large [34].⁴ While of physical primacy, the discrete model can be less suitable for reasoning about feasible chemical algorithms. Many algorithms in the discrete model rely on a single molecule (called a “leader”) to coordinate computation [8]. Whether the initial state is assumed to have a leader or the CRN is designed to eliminate all but one copy of the leader species (“leader election”), such algorithms relying on single-molecule behavior are currently infeasible, since any single molecule can get damaged or become effectively lost.

An important reason for our focus on stoichiometric computation is that algorithms relying only on stoichiometry make easier design targets. The rates of reactions are real-valued quantities that can fluctuate with reaction conditions such as temperature, while the stoichiometric coefficients are immutable whole numbers set by the nature of the reaction. Methods for physically implementing CRNs naturally yield systems with digital stoichiometry that can be set exactly [14, 50], whereas these methods often suffer from imprecise control over reaction rates [20, 51]. Further, relying on specific rate laws can be problematic: Many systems do not apparently follow mass-action rate laws, and chemists have developed an array of alternative rate laws such as Michaelis-Menten (modeling enzymes) and Hill-function kinetics (widely used for gene regulation).⁵ It is well-known that cells are not well-mixed, and many models have been developed to take space into account (e.g., reaction-diffusion [32]). Moreover, robustness of rate laws is a recurring motif in systems biology due to much evidence that biological regulatory networks tend to be robust to the form of the rate laws and the rate parameters [10]. Thus, we are interested in what computations can be understood or engineered without regard for the reaction rate laws.

1.3 Related Works

An earlier conference version of this article appeared as Reference [19]. Besides replacing a number of informal arguments with rigorous proofs, this journal version also expands and generalizes the results of the conference version. For example, we introduce new machinery for representing and manipulating trajectories as linear objects (piecewise linear paths). We also define a broad class of rate laws, formalized by Definition 2.22, which captures mass-action kinetics and all other known rate laws such as Michaelis-Menten and Hill-function kinetics, and we prove that our definition of reachability is as general as any in this class. For the constructive part, this version also generalizes Lemma 3.4 of Reference [19] (in addition to correcting its proof) by introducing feedforward CRNs and proving that correct computation in our setting implies convergence under any “reasonable” rate law (one that produces a fair schedule of rates; definition 4.1) for any feedforward CRN (lemma 4.8).

The relationship between the discrete and continuous CRN models is a complex and much studied one in the natural sciences [45]. The computational abilities of discrete CRNs have been investigated more thoroughly than of continuous CRNs and have been shown to have a surprisingly rich computational structure. Of most relevance here is the work in the discrete setting showing that the class of functions that can be computed depends strongly on whether the computation must be correct or just likely to be correct (under the usual stochastic kinetics)—which is the discrete version of the distinction between rate-independent and rate-dependent computation. While Turing universal computation is possible with an arbitrarily small, non-zero probability of error over all time [49], forbidding error altogether limits the computational power: Error-free computation by stochastic CRNs is limited to semilinear predicates and functions [6, 18]. (Intuitively, semilinear functions are expressible as a finite union of affine functions, with “simple, periodic” domains of each affine function [18].) The study of error-free computation in discrete CRNs is heavily based on the results first developed for a model of distributed computing called population protocols [6, 9]. We formally refer to our notion of rate-independent computation as stable computation in direct reference to the analogous notion in population protocols.

While our notion of rate-independent computation is the natural extension of deterministic computation in the discrete model, there are many differences between the two settings. As mentioned above, many discrete algorithms such as those that rely on a single “leader” molecule fail to work in the continuous setting, and some functions like distinguishing between even and odd molecular counts do not make sense. Broadly speaking, the proof techniques appear to require very different machinery, and the importance of stable computation itself needs substantial justification in the continuous model (as the examples shown at the beginning of Section 4 demonstrate).

Continuous CRNs have been proven to be Turing universal under mass-action rate laws [27], a consequence of the surprising computational power of polynomial ODEs [12]. In ODEs without the CRN semantics, there is no natural notion of stoichiometry and thus no notion of rate-independence analogous to ours. In chemistry, the same physical process (a reaction) is responsible for multiple monomials across multiple ODEs, which justifies these monomials being exactly the same or in fixed ratios (corresponding to obeying reaction stoichiometry). Such forced relationships do not seem natural for more general polynomial ODEs that do not correspond to chemical reactions.⁶

Our notion of reachability (definition 2.3) is intended to capture a wide diversity of possible rate laws. Generalized rate laws (extending mass-action, Michaelis-Menten, etc.) have been previously studied, although not in a computational setting. For example, certain conditions were identified on global convergence to equilibrium based on properties intuitively similar to ours [4]. A related idea in the literature, generalizing mass-action, is differential inclusion [30]. In that model, the mass-action rate constants are not fixed to be particular real numbers constant over time, but instead can vary over time within some bounded interval \([l,u]\) fixed in advance, with \(0 \lt l \le u \lt \infty\). Another related idea is the notion of a reaction system [26], which generalizes even beyond mass-action, allowing reaction rates to be an (almost) arbitrary function of species concentrations.⁷ Other generalized rate laws have been defined as well [5, 23].

Since the original publication of the conference version of this article [19], a number of works have used our framework. A key concept in capturing rate-independent computation is the reachability relation (segment-reachability, definition 2.3). Reference [16] showed that, given two states, deciding whether one is reachable from the other is solvable in polynomial time. This contrasts sharply with the hardness of the reachability problem for discrete CRNs, which, although computable [39], is not even primitive recursive [22, 37]. (These results were proven using the terminology of the equivalent models of Petri nets/vector addition systems.)

The question of deciding whether a given CRN is rate-independent was studied in Reference [23]. The work provides sufficient graphical conditions on the structure of the CRN that ensure rate-independence for the whole CRN or only for certain output species. Interestingly, the authors of Reference [23] applied this method to the Biomodels repository of curated CRNs of biological origin and found a number of CRNs that satisfy the rate-independence conditions.

An important motivation for the dual-rail representation in this work is to allow composition of rate-independent CRN modules (section 5.1). Such rate-independent modules can be composed into overall rate-independent computation simply by concatenating their chemical reactions and relabeling species (such that the output species of the first is the input species of the second, and all other species are distinct). In contrast, rate-independent composition with the direct (non-dual rail) representation, introduces an additional “superadditivity” constraint that for all input vectors \(\mathbf {x}\) and \(\mathbf {x}^{\prime }\), \(f(\mathbf {x}) + f(\mathbf {x}^{\prime }) \le f(\mathbf {x}+ \mathbf {x}^{\prime })\) [17]. Thus, for example, the non-superadditive max function (Figure 1) provably cannot be composably computed with a rate-independent CRN in the direct representation. Composable computation has also been characterized in the discrete model [31, 48].

Other input encodings have been considered besides direct and dual-rail. For example, the so-called “fractional encoding” encodes a real number between 0 and 1 as a ratio \(\frac{x_1}{x_0+x_1}\) where \(x_0,x_1\) are concentrations of two input species [44]. Other notions of chemical “rate-independence” include CRNs that work independently of the rate law as long as there is a separation into fast and slow reactions [47]. For a detailed survey on computation with CRNs (both continuous and discrete), see Reference [13].

2.1 Chemical Reaction Networks

We first explain our notation for vectors of concentrations of chemical species and then formally define chemical reaction networks.

Given a finite set F, let \(\mathbb {R}^F\) denote the set of functions \(\mathbf {c}: F \rightarrow \mathbb {R}\). We view \(\mathbf {c}\) equivalently as a vector of real numbers indexed by elements of F. Given \(x \in F\), we write \(\mathbf {c}(x)\), or sometimes \(\mathbf {c}_x\), to denote the real number indexed by x. The notation \(\mathbb {R}_{\ge 0}^F\) is defined similarly for nonnegative real vectors. Throughout this article, let \(\Lambda\) be a finite set of chemical species. Given \(S\in \Lambda\) and \(\mathbf {c}\in \mathbb {R}_{\ge 0}^\Lambda\), we refer to \(\mathbf {c}(S)\) as the concentration of S in \(\mathbf {c}\). For any \(\mathbf {c}\in \mathbb {R}_{\ge 0}^\Lambda\), let \([\mathbf {c}] = \lbrace S \in \Lambda \ |\ \mathbf {c}(S) \gt 0 \rbrace\), the set of species present in \(\mathbf {c}\) (a.k.a., the support of \(\mathbf {c}\)). We write \(\mathbf {c}\le \mathbf {c}^{\prime }\) to denote that \(\mathbf {c}(S) \le \mathbf {c}^{\prime }(S)\) for all \(S \in \Lambda\). Given \(\mathbf {c},\mathbf {c}^{\prime } \in \mathbb {R}_{\ge 0}^\Lambda\), we define the vector component-wise operations of addition \(\mathbf {c}+\mathbf {c}^{\prime }\), subtraction \(\mathbf {c}-\mathbf {c}^{\prime }\), and scalar multiplication \(x \mathbf {c}\) for \(x \in \mathbb {R}\). If \(\Delta \subset \Lambda\), we view a vector \(\mathbf {c}\in \mathbb {R}_{\ge 0}^\Delta\) equivalently as a vector \(\mathbf {c}\in \mathbb {R}_{\ge 0}^\Lambda\) by assuming \(\mathbf {c}(S)=0\) for all \(S \in \Lambda \setminus \Delta .\) For \(\Delta \subset \Lambda\), we write \(\mathbf {c}\upharpoonright \Delta\) to denote \(\mathbf {c}\) restricted to \(\Delta\); in particular, \(\mathbf {c}\upharpoonright \Delta = \mathbf {0} \iff (\forall S\in \Delta)\ \mathbf {c}(S)=0.\) (We use the convention that \(\mathbf {c}\upharpoonright \emptyset = \mathbf {0}\) for all states \(\mathbf {c}\).)

A reaction over \(\Lambda\) is a pair \(\alpha = \langle \mathbf {r},\mathbf {p}\rangle \in \mathbb {N}^\Lambda \times \mathbb {N}^\Lambda\), such that \(\mathbf {r}\ne \mathbf {p}\), specifying the stoichiometry of the reactants and products, respectively.⁸ For instance, given \(\Lambda =\lbrace A,B,C\rbrace\), the reaction \(A+2B \rightarrow A+3C\) is the pair \(\left\langle (1,2,0) , (1,0,3) \right\rangle .\) We represent reversible reactions such as \(A \mathop {\rightleftharpoons }\limits B\) as two irreversible reactions \(A \rightarrow B\) and \(B \rightarrow A\). In this article, we assume that \(\mathbf {r}\ne \mathbf {0}\), i.e., we have no reactions of the form \(\emptyset \rightarrow \ldots\).⁹ A (finite) chemical reaction network (CRN) is a pair \(\mathcal {C}=(\Lambda ,R)\), where \(\Lambda\) is a finite set of chemical species, and R is a finite set of reactions over \(\Lambda\). A state of a CRN \(\mathcal {C}=(\Lambda ,R)\) is a vector \(\mathbf {c}\in \mathbb {R}_{\ge 0}^\Lambda\). Given a state \(\mathbf {c}\) and reaction \(\alpha =\left\langle \mathbf {r} , \mathbf {p} \right\rangle\), we say that \(\alpha\) is applicable in \(\mathbf {c}\) if \([\mathbf {r}] \subseteq [\mathbf {c}]\) (i.e., \(\mathbf {c}\) contains positive concentration of all of the reactants). If no reaction is applicable in state \(\mathbf {c}\), then we say \(\mathbf {c}\) is static. We say a species S is produced in reaction \(\langle \mathbf {r},\mathbf {p}\rangle\) if \(\mathbf {r}(S) \lt \mathbf {p}(S)\), and consumed if \(\mathbf {r}(S) \gt \mathbf {p}(S)\). (Note that a catalyst, such as C in the reaction \(C+X \rightarrow C+Y\), is neither produced nor consumed.)

2.2 Segment Reachability

In the previous section, we defined the syntax of CRNs. Toward studying rate-independent computation, we now want to define the semantics of what “could happen” if reaction rates can vary arbitrarily over time. This is captured by a notion of reachability, which is the focus of this section. Intuitively, \(\mathbf {d}\) is reachable from \(\mathbf {c}\) if applying some amount of reactions to \(\mathbf {c}\) results in \(\mathbf {d}\), such that no reaction is ever applied when any of its reactants are concentration 0. Formalizing this concept is a bit tricky and constitutes one of the contributions of this article. Intuitively, we will think of reachability via straight line segments. This may appear overly limiting; after all mass-action and other rate laws trace out smooth curves. However, in this and subsequent sections, we show a number of properties of our definition that support its reasonableness.

Throughout this section, fix a CRN \(\mathcal {C}= (\Lambda ,R)\). All states \(\mathbf {c}\), and so on, are assumed to be states of \(\mathcal {C}\). We define the \(|\Lambda | \times |R|\) reaction stoichiometry matrix \(\mathbf {M}\) such that, for species \(S \in \Lambda\) and reaction \(\alpha = \langle \mathbf {r},\mathbf {p}\rangle \in R\), \(\mathbf {M}(S,\alpha) = \mathbf {p}(S) - \mathbf {r}(S)\) is the net amount of S produced by \(\alpha\) (negative if S is consumed).¹⁰ For example, if we have the reactions \(X \rightarrow Y\) and \(X + A \rightarrow 2X + 3Y\), and if the three rows correspond to X, A, and Y, in that order, then \(\begin{equation*} \mathbf {M}= \left(\begin{array}{cc} -1 & 1 \\ 0 & -1 \\ 1 & 3 \\ \end{array}\right). \end{equation*}\)

Intuitively, by a single segment, we mean running the reactions applicable at \(\mathbf {c}\) at a constant (possibly 0) rate to get from \(\mathbf {c}\) to \(\mathbf {d}\). In the definition, \(\mathbf {u}(\alpha)\) represents the flux of reaction \(\alpha \in R\).

The next definition is used in our main notion of reachability, which uses either a finite number of straight lines or infinitely many so long as they converge to a single state.

For example, suppose the reactions are \(X \rightarrow C\) and \(C + Y \rightarrow C + Z\), and we are in state \(\lbrace 0 C, 1 X, 1 Y, 0 Z\rbrace\). With straight-line segments, any state with a positive amount of Z must be reached in at least two segments: first to produce C, which allows the second reaction to occur, and then any combination of the first and second reactions. For example, \(\lbrace 0 C\), \(1 X\), \(1 Y\), \(0 Z\rbrace\) \(\rightarrow ^1\) \(\lbrace 0.1 C\), \(0.9 X\), \(1 Y\), \(0 Z\rbrace\) \(\rightarrow ^1\) \(\lbrace 1 C\), \(0 X\), \(0 Y\), \(1 Z\rbrace\). This is a simple example showing that more states are reachable with \(\rightsquigarrow\) than \(\rightarrow ^1\). Often, Definition 2.3 is used implicitly, when we make statements such as, “Run reaction 1 until X is gone, then run reaction 2 until Y is gone,” which implicitly defines two straight lines in concentration space.

Although more effort will be needed to justify its reasonableness (see Section 2.4), segment-reachability will serve as the main notion of reachability in this article.

2.3 Bound on Number of Required Line Segments in Segment Reachability

It may seem that we can never achieve the “full diversity” of states reachable with an infinite number of line segments if we use only a bounded number of line segments. However, Theorem 2.15 shows that increasing the number of straight-line segments beyond a certain point does not make any additional states reachable. Thus, using a few line segments captures all the states reachable with arbitrarily many line segments, and in fact even in the limit of infinitely many line segments.

To prove Theorem 2.15, we first develop important machinery for representing and manipulating paths under \(\rightsquigarrow\). Note that reachability is closed under addition and scaling in the sense that if \(\mathbf {c}\rightsquigarrow \mathbf {d}\) and \(\mathbf {c}^{\prime } \rightsquigarrow \mathbf {d}^{\prime }\), then \(\alpha \mathbf {c}+ \beta \mathbf {c}^{\prime } \rightsquigarrow \alpha \mathbf {d}+ \beta \mathbf {d}^{\prime }\) for all \(\alpha , \beta \in \mathbb {R}_{\ge 0}\). The following definition captures this property by defining a linear space of all paths. This machinery will also be key to proving the piecewise linearity of the computed function in Section 5.5.

Definition 2.4 allows a prepath to be essentially any sequence of vectors in the linear subspace spanned by reaction vectors. The next definition restricts the vectors with three physical constraints: Species concentrations are nonnegative, reaction fluxes are nonnegative (i.e., reactions can only go one way, turning reactants into products), and reactions cannot occur if any reactant is 0.

Definitions 2.4 and 2.5 allow an infinite sequence of reaction flux vectors (each corresponding to a straight line in the definition of 1-segment reachability). A finite number of straight lines can be specified by letting all but finitely many \(\mathbf {u}_i = \mathbf {0}\). The next definition bounds how many can be nonzero.

Intuitively, \(\Gamma _\infty\) is the space of all of the valid piecewise linear paths that the system can take starting from any given initial state and \(\Gamma _k\) (\(k \in \mathbb {N}\)) is the set of all such paths that have length at most k; thus, \(\Gamma _0 \subseteq \Gamma _1 \subseteq \ldots \subseteq \Gamma _\infty\).

The next lemma shows that if it is possible to reach from a state \(\mathbf {c}\) to several other states, each containing some species possibly distinct from each other, then it is possible to reach from \(\mathbf {c}\) to a state with all of those species present at once.

The next lemma shows that with at most a constant number of straight line segments, it is possible to reach from any state \(\mathbf {c}\) to a state containing all species possible to produce from \(\mathbf {c}\).

Recall that a set is closed if it contains all of its limit points.

Note that Lemma 2.11 is false if we replace “straight-line reachable” with “segment-reachable.” For example, consider \(X \rightarrow C\) and \(X + C \rightarrow Y + C\), where we take the initial state \(\mathbf {c}= \lbrace 1X, 0C, 0Y\rbrace\). Note that for any \(\varepsilon \gt 0\), we can reach the state \(\mathbf {d}_\varepsilon = \lbrace 0 X, \varepsilon C, (1 - \varepsilon) Y \rbrace\). However, because producing Y first requires consuming a positive amount of X to create the catalyst C, the state \(\mathbf {d}_0 = \lbrace 0X, 0C, 1 Y\rbrace\) is not segment reachable from \(\mathbf {c}\), even though \(\mathbf {d}_0 = \lim _{\varepsilon \rightarrow 0} \mathbf {d}_{\varepsilon }\).

If we have an infinite sequence of states such that \(\mathbf {c}= \mathbf {b}_0 \rightsquigarrow \mathbf {b}_1 \rightsquigarrow \mathbf {b}_2 \dots\) and \(\lim _{i \rightarrow \infty } \mathbf {b}_i = \mathbf {d}\), then this does not immediately imply that \(\mathbf {c}\rightsquigarrow ^\infty \mathbf {d}\) by definition 2.2. This is because although the endpoints of the paths \(\mathbf {b}_i \rightsquigarrow \mathbf {b}_{i+1}\) converge to \(\mathbf {d}\), the intermediate states on the paths related by \(\rightarrow ^1\) (i.e., \(\mathbf {b}_i \rightarrow ^1 \mathbf {b}_i^{\prime } \rightarrow ^1 \mathbf {b}_i^{\prime \prime } \rightarrow ^1 \dots \rightarrow ^1 \mathbf {b}_{i+1}\)) may not converge. To capture this weaker notion of convergence, we introduce the following definition, which generalizes \(\rightsquigarrow ^\infty\) by requiring only that there be a converging subsequence of states. (The weaker notion will be eventually needed to prove theorem 3.3.)

The main results of this section have this weaker notion of convergence as a precondition, which will imply that, despite appearances, \(\rightsquigarrow ^\infty\) and \(\rightsquigarrow ^\infty _{\mathrm{ss}}\) are actually equivalent.

The next lemma shows that if no more species are producible from state \(\mathbf {c}\) than are already present in \(\mathbf {c}\), then any state \(\mathbf {d}\) that is \(\rightsquigarrow ^\infty _{\mathrm{ss}}\) reachable from \(\mathbf {c}\) is reachable via a single straight line segment.

Finally, the previous lemmas can be combined to show that at most a constant number of straight line segments (depending on the CRN) are required to reach from any state to any other reachable state. In fact, this holds even for states that are only \(\rightsquigarrow ^\infty _{\mathrm{ss}}\) reachable.

Note that theorem 2.14 immediately implies that \(\rightsquigarrow ^\infty\) and \(\rightsquigarrow ^\infty _{\mathrm{ss}}\) are the same relation, since \(\mathbf {c}\rightsquigarrow ^\infty _{\mathrm{ss}}\mathbf {d}\) implies \(\mathbf {c}\rightsquigarrow ^{m+1} \mathbf {d}\) implies \(\mathbf {c}\rightsquigarrow ^{\infty } \mathbf {d}\).

Although the full power of Theorem 2.14 is useful later in Theorem E.1, the most important consequence of Theorem 2.14 is the following result, which we will use repeatedly:

Corollary 2.15. will allow us to assume without loss of generality that there are a constant number of line segments between any two states, simplifying many arguments. For example, it is not obvious that the relation \(\rightsquigarrow ^\infty\) is transitive, since one cannot concatenate two infinite sequences. However, since two finite sequences of segments can be concatenated, the following corollary is immediate:

The goal of the reachability relation is to capture “what could happen” in chemical reaction networks independently of rates. Thus, it is natural to satisfy several properties: The relation should be reflexive (true for \(\rightsquigarrow\), since \(\mathbf {x}\rightsquigarrow ^0 \mathbf {x}\)) and transitive (Corollary 2.16). Further, the relation should be additive in the intuitive sense that the presence of additional molecules cannot entirely prevent reactions from happening (although in a kinetic model it could effectively slow down reactions due to competition); formally, if \(\mathbf {x}\rightsquigarrow \mathbf {y}\), then \(\mathbf {x}+ \mathbf {c}\rightsquigarrow \mathbf {y}+ \mathbf {c}\) for any state \(\mathbf {c}\). Additivity is a crucial property of the more standard notion of discrete CRN reachability, used for example in many cases to prove impossibility results for those systems [1, 2, 7, 11, 25]. We also employ additivity of \(\rightsquigarrow\) for impossibility results, for example, in the proofs of Lemmas 5.30 and 5.31.¹¹

While satisfying these properties is a good start for justifying the reasonableness of \(\rightsquigarrow\), it is natural to wonder whether there are some “reasonable” rate laws that segment-reachability fails to capture, i.e., perhaps some CRN rate law would take \(\mathbf {x}\) to \(\mathbf {y}\) even though \(\mathbf {x}\not\rightsquigarrow \mathbf {y}\). In Section 2.4, we define an apparently much more general notion of reachability (Definition 2.22) that captures all commonly studied rate laws, while still respecting the fundamental semantics of reactions. We prove that our reachability relation is in fact identical to it, i.e., \(\mathbf {x}\) can reach to \(\mathbf {y}\) under this notion if and only if \(\mathbf {x}\rightsquigarrow \mathbf {y}\) (Theorem 2.27 and Lemma 2.26).

2.4 Generality of Segment Reachability

In this section, we justify that our notion of reachability via straight lines actually corresponds to the most general notion of “being able to get from one state to another,” restricted only by the non-negativity of concentrations, reaction stoichiometry and the need for catalysts—as long as we maintain the causal relationships between the production of species. This notion of reachability admits “time-varying” rate laws where reactions occur according to some arbitrary schedule, which captures situations such as solutions that are not well-mixed, or where physical parameters, such as temperature, change in some arbitrary way. The main result of this section, Theorem 2.27, formalizes this idea and justifies calling segment-reachability simply “reachability” in the rest of this article.

We begin with a review of mass-action kinetics, the most commonly used rate law in chemistry and show the (physically and intuitively obvious but mathematically subtle) features that make it consistent with segment reachability. We then generalize rate-law trajectories to arbitrary “valid rate schedules” and prove that these are exactly captured by segment-reachability.

A CRN with positive rate constants assigned to each reaction defines a mass-action ODE (ordinary differential equation) system with a variable for each species, which represents the time-varying concentration of that species. We follow the convention of upper-case species names and lower-case concentration variables. Each reaction contributes one term to the ODEs for each species produced or consumed in it. The term from reaction \(\alpha\) appearing in the ODE for x is the product of: the rate constant, the reactant concentrations, and the net stoichiometry of X in \(\alpha\) (i.e., the net amount of X produced by \(\alpha\), negative if consumed). For example, the CRN \(\begin{align*} X + X &\mathop {\rightarrow }\limits ^{k_1} C \\ C + X &\mathop {\rightarrow }\limits ^{k_2} C + Y \end{align*}\) corresponds to ODEs: (2.1) \(\begin{align} dx/dt &= -2 k_1 x^2 - k_2 c x, \end{align}\) (2.2) \(\begin{align} dc/dt &= k_1 x^2, \end{align}\) (2.3) \(\begin{align} dy/dt &= k_2 c x, \end{align}\) where \(k_1, k_2\) are the rate constants of the two reactions.

Given a CRN \(C = (\Lambda , R)\), let \(\mathbf {A}(\mathbf {d}): \mathbb {R}_{\ge 0}^\Lambda \rightarrow \mathbb {R}_{\ge 0}^R\) be the rates of all the reactions in state \(\mathbf {d}\) as given by the mass-action ODEs. Given an assignment of (strictly) positive rate constants, and an initial state \(\mathbf {c}\), the mass-action trajectory is a function \({\boldsymbol \rho }: [0,t_\mathrm{max}) \rightarrow \mathbb {R}_{\ge 0}^\Lambda\), where \(t_\mathrm{max} \in \mathbb {R}_{\ge 0}\cup \lbrace \infty \rbrace\), such that \({\boldsymbol \rho }\) is the solution to \(d {\boldsymbol \rho }/dt = \mathbf {M}\cdot \mathbf {A}({\boldsymbol \rho }(t))\) with \({\boldsymbol \rho }(0) = \mathbf {c}\), where \(t_\mathrm{max}\) is the maximum time, typically \(\infty\), for which the solution is defined on all of \([0,t_\mathrm{max})\). Although beyond the scope of this article, mass-action ODEs are locally Lipschitz, so a CRN admits exactly one mass-action trajectory for a fixed collection of rate constants and initial state \(\mathbf {c}\). Note that for some CRNs (e.g., \(2X \rightarrow 3X\)), the solution of the ODEs goes to infinite concentration in finite time,¹² and for such CRNs, \(t_\mathrm{max}\) is finite.

To prove Theorem 2.27, we need to introduce the notion of a siphon from the Petri net literature. This notion will be used, as well, to prove negative results in Section 5.5.

The following lemma, due to Angeli, De Leenheer, and Sontag [5], shows that this is equivalent to the notion that “the absence of \(\Omega\) is forward-invariant” under mass-action: If all species in \(\Omega\) are absent, then they can never again be produced (under mass-action).¹⁴ For the sake of completeness, we give a self-contained proof in Appendix A.

We show that the same holds true for segment-reachability. Due to the discrete nature of segment-reachability, the proof is more straightforward than that of Lemma 2.19. It follows the same essential structure one would use to prove this in the discrete CRN model: If the siphon \(\Omega\) is absent, then no reaction with a reactant in \(\Omega\) can be the next reaction to fire, so by the siphon property, no species in \(\Omega\) is produced in the next step.

Recall that \(\mathsf {P}(\mathbf {c})\) represents the set of species producible from state \(\mathbf {c}\). The next lemma shows that the set of species that cannot ever be produced from a given state is a siphon.

The main result of this section is Theorem 2.27, which justifies that our (seemingly limited) notion of reachability via straight lines is actually quite general. To state the theorem, we define a very general notion of “reasonable rate laws,” which are essentially schedules of rates to assign to reactions over time. All known rate laws, such as mass-action, Michaelis-Menten, Hill function kinetics, as well as our own nondeterministic notion of adversarial rates following straight lines (segment reachability, Definition 2.3), obey this definition. (We justify this below explicitly for mass-action and Definition 2.3, but it is straightforward to verify in the other cases.)

Recall that R is the set of all reactions in some CRN, and \(\Lambda\) is the set of its species.

We note that because \(\mathbf {f}\) is Lebesgue integrable, by Reference [43, Theorem 6.11], \(\mathbf {F}_\alpha (t)\) is locally absolutely continuous.

Although Definition 2.22 explicitly constrains the states \({\boldsymbol \rho }(t)\) to be non-negative, the non-negativity of \({\boldsymbol \rho }(t)\) actually follows from condition (2).¹⁶

Conditions (2) and (3) may appear redundant, but in fact each can be obeyed while the other is violated.

For example, consider the reaction \(\alpha : X \rightarrow 2X\), starting in the state \(\lbrace 0 X\rbrace\), with invalid rate schedule \(\mathbf {f}_\alpha (t) = t\), with trajectory \({\boldsymbol \rho }_X(t) = t^2 / 2\). Since the rate \(\mathbf {f}_\alpha (0)\) is 0, this vacuously satisfies (2) at time 0, and, since \({\boldsymbol \rho }_X(t) \gt 0\) for \(t \gt 0\) (X is present at all positive times), (2) is also satisfied for positive times. However, (3) is violated, since \(\lbrace X\rbrace\) is a siphon absent at time 0 but present at future times. This example also demonstrates why condition (3) is required to satisfy our intuitive understanding of reasonable rate laws respecting “causality of production” among species: With only the reaction \(X \rightarrow 2X\), the only way to produce more X is already to have some X.

To see the other case, take reactions \(\alpha : X \rightarrow C\) and \(\beta : C+X \rightarrow C+Y\), starting in state \(\lbrace 1 X, 1 Y, 0 C\rbrace\). Consider the invalid rate schedule \(\mathbf {f}_\alpha (t) = 0\) for all t, \(\mathbf {f}_\beta (t) = 1\) for \(0 \le t \le 1/2\), and \(\mathbf {f}_\beta (t)=0\) for \(t \gt 1/2\), i.e., run only \(\beta\), until X is half gone. This violates (2), since \(\beta\) occurs without its reactant C present. However, the only set of species absent along this trajectory is \(\lbrace C\rbrace\), which is not a siphon, since reaction \(\alpha\) has C as a product but not a reactant, so (3) is vacuously satisfied.

The next lemma shows that the most commonly used rate law, mass-action, gives a valid rate schedule and trajectory according to Definition 2.22. Recall that \(\mathbf {A}({\boldsymbol \rho }(t)): \mathbb {R}_{\ge 0}^\Lambda \rightarrow \mathbb {R}_{\ge 0}^R\) represents the rates of all the reactions in state \({\boldsymbol \rho }(t)\) as given by the mass-action ODEs. For instance, for our mass-action example at the beginning of this section, the function \(\mathbf {f}(t) = \mathbf {A}({\boldsymbol \rho }(t))\) corresponds to the ODEs of Equations (2.1)–(2.3) when written as \(d{\boldsymbol \rho }/dt = \mathbf {M}\cdot \mathbf {f}(t)\).

We say that a rate schedule \(\mathbf {f}: \mathbb {R}_{\ge 0}\rightarrow \mathbb {R}_{\ge 0}^R\) is finite if there is \(t_0 \ge 0\) such that \(\mathbf {f}(t) = \mathbf {0}\) for all \(t \ge t_0\), i.e., reactions eventually stop occurring. The following observation is straightforward to verify, showing that the concatenation of two valid rate schedules, with the first finite, is also a valid rate schedule.

The next lemma shows essentially that our definition of segment-reachability creates a valid rate schedule.

Finally, we have the main result of this section, which shows that segment reachability is as general as any valid rate schedule.

Recall mass-action trajectories correspond to valid rate schedules by Lemma 2.24. Thus, Theorem 2.27 implies the following corollary, which intuitively says that if a state is reachable via a mass-action trajectory (even in the limit of infinite time), then it is segment-reachable.¹⁷

4.1 Feedforward CRNs

As we saw before, for general CRCs it is possible that the output-stable state is always reachable but the mass-action trajectory does not converge to it, for instance, the reactions \(X + X \rightarrow Y + Y\) and \(Y + X \rightarrow X + X\) discussed at the start of Section 4. Thus, stable computation does not necessarily imply that the system will eventually produce the correct output. Since the mass-action trajectory defines a fair rate schedule (lemma 4.2), the above example shows that some CRCs stably compute a function but do not fairly compute it; i.e., the converse of lemma 4.4 does not hold for all CRCs.

In contrast to the above example, the feedforward property defined in this section allows us to bridge the definition of stable computation, defined in terms of reachability (what could happen), to convergence (what will happen), defined in terms of fair rate schedules.

Recall that a reaction \(\alpha =\langle \mathbf {r},\mathbf {p}\rangle\) produces a species S if \(\mathbf {r}(S) \lt \mathbf {p}(S)\) and consumes S if \(\mathbf {r}(S) \gt \mathbf {p}(S)\). We say a CRN is feedforward if the species can be ordered so every reaction that produces a species also consumes another species earlier in the ordering. Formally:

Intuitively, we want to avoid situations, as in the example above, where the output-stable state is always reachable but the trajectory does not converge to it. This can happen if the trajectory does not converge or converges to a dynamic equilibrium where reactions balance each other. In contrast, the total flux of reactions in a feedforward CRN must be bounded, because there cannot be a complete “cycle” among the species that balance consumption with production.

We start with the following simple observation: General CRNs can have reactions such as \(A+B \rightarrow A+2B+C\) that do not consume any reactant, and such CRNs clearly have infinite total flux. Luckily, by our definition of a CRN, any reaction \(\alpha\) must either produce or consume some species, and if \(\alpha\) produces a species, the feedforward condition guarantees that \(\alpha\) consumes some other species.

The following lemma will used in formalizing the idea that a feedforward CRN cannot converge to a dynamic equilibrium (like the example beginning this section). General CRNs can have reactions that undo each other’s effect (for instance, \(A \rightarrow B\) and \(B \rightarrow A\)). For such CRNs, we cannot bound total flux as a function of the change in species concentration—indeed, concentrations might remain constant but the two reactions canceling each other can have arbitrarily large flux—allowing for a dynamic equilibrium. In contrast, for feedforward CRNs, the following lemma shows that total flux can be bounded by the change in state:

Finally, we are ready to prove the main lemma about feedforward CRNs, a consequence of which (Lemma 4.9) establishes the connection between stable computation and convergence to the output-stable state for feedforward CRNs.

Proof.

First, define \(\begin{equation*} K = 1 + \max _{\alpha \in R} \sum _{S_i \in \Lambda } \max (0, \mathbf {M}(S_i,\alpha)). \end{equation*}\) In particular, K is larger than the sums of the positive entries in the columns of \(\mathbf {M}\). Intuitively, we will assign a “mass” to each species, and K will be the ratio of the masses assigned to consecutive species in the feedforward ordering. By making K sufficiently large, we can guarantee that running any reaction has the effect of decreasing the total mass \(V(\mathbf {x})\) of any state \(\mathbf {x}\). Formally, define \(V: \mathbb {R}_{\ge 0}^\Lambda \rightarrow \mathbb {R}_{\ge 0}\) as \(\begin{equation*} V(\mathbf {x}) = \sum _{i =1}^n \frac{\mathbf {x}(S_i)}{K^i}. \end{equation*}\)

Recall \(\mathbf {f}_\alpha (t)\) is the rate of reaction \(\alpha\) at time t. Note that whenever \({\boldsymbol \rho }(t)\) is differentiable²⁰ \(\begin{align*} \frac{d}{dt} V({\boldsymbol \rho }(t)) &= \sum _{i =1}^n \frac{d}{dt}\frac{({\boldsymbol \rho }(t))(S_i)}{K^i} \\ &= \sum _{i =1}^n \frac{1}{K^i} \sum _{\alpha \in R} \mathbf {M}(S_i,\alpha) \mathbf {f}_\alpha (t) \\ &= \sum _{\alpha \in R} \sum _{i =1}^n \frac{1}{K^i} \mathbf {M}(S_i,\alpha) \mathbf {f}_\alpha (t) \\ &= \sum _{\alpha \in R} C_\alpha \mathbf {f}_\alpha (t), \end{align*}\) where \(C_\alpha = \sum _{i =1}^n \frac{1}{K^i} \mathbf {M}(S_i,\alpha)\). For any fixed \(\alpha\), let \(i_0\) be the smallest i such that \(\mathbf {M}(S_i, \alpha) \ne 0\), i.e., \(S_{i_0}\) is the first species in the feedforward ordering that is produced or consumed by \(\alpha\). By the feedforward condition, \(\mathbf {M}(S_{i_0}, \alpha) \le -1\), i.e., \(S_{i_0}\) is consumed. As a result, \(\begin{align*} C_\alpha =&\, \frac{1}{K^{i_0}}\left(\mathbf {M}(S_{i_0}, \alpha) + \frac{1}{K} \sum _{i = i_0 + 1}^n \frac{\mathbf {M}(S_i,\alpha)}{K^{i - i_0 - 1}}\right) \\ &\lt \frac{1}{K^{i_0}}\left(\mathbf {M}(S_{i_0}, \alpha) + 1\right) \\ &\le 0. \end{align*}\)

Now, we show that the total flux through all of the reactions is finite. To fix notation, let us write \(\mathbf {F}_\alpha\) for the total flux through reaction \(\alpha\) as \(t \rightarrow \infty\), i.e., \(\begin{equation*} \mathbf {F}_\alpha = \int _0^\infty \mathbf {f}_\alpha (t) dt. \end{equation*}\) Since every \(C_\alpha\) is negative, for any fixed reaction \(\alpha _0\), \(\begin{align*} \frac{d}{dt} V({\boldsymbol \rho }(t)) &= \sum _{\alpha \in R} C_\alpha \mathbf {f}_\alpha (t) \le C_{\alpha _0} \mathbf {f}_{\alpha _0}(t) \end{align*}\) so by integrating²¹ both sides \(\begin{align*} \lim _{t \rightarrow \infty } V({\boldsymbol \rho }(t)) - V(\mathbf {c}) &\le C_{\alpha _0} \int _0^\infty \mathbf {f}_{\alpha _0}(t) dt = C_{\alpha _0}\mathbf {F}_{\alpha _0}. \end{align*}\) If \(\mathbf {F}_{\alpha _0}\) were infinite, then because \(C_{\alpha _0} \lt 0\), we see that \(V({\boldsymbol \rho }(t))\) would be unbounded below. Since \({\boldsymbol \rho }(t)\) always remains in \(\mathbb {R}_{\ge 0}^\Lambda\) where \(V({\boldsymbol \rho }(t))\)\(\ge 0\), we conclude that \(\mathbf {F}_{\alpha _0}\) must be finite.

If we write \(\mathbf {v}_\alpha \in \mathbb {R}^\Lambda\) for the \(\alpha\)-column of \(\mathbf {M}\) (i.e., \(\mathbf {v}_\alpha (S) = \mathbf {M}(S, \alpha)\)), then \(\begin{align*} \int _0^\infty \left\Vert \frac{d{\boldsymbol \rho }(t)}{dt}\right\Vert dt =&\, \int _0^\infty \left\Vert \sum _{\alpha \in R} \mathbf {v}_\alpha \mathbf {f}_\alpha (t)\right\Vert dt \\ &\le \int _0^\infty \sum _{\alpha \in R} \mathbf {f}_\alpha (t)|\mathbf {v}_\alpha | dt \\ =&\, \sum _{\alpha \in R}|\mathbf {v}_\alpha | \int _0^\infty \mathbf {f}_\alpha (t) dt = \sum _{\alpha \in R}|\mathbf {v}_\alpha | \mathbf {F}_\alpha \\ &\lt \infty . \end{align*}\) In other words, \({\boldsymbol \rho }\) has finite length.

Let us now show that \({\boldsymbol \rho }\) having finite length implies that it converges. For any states \(\mathbf {x},\mathbf {z}\) define \(d(\mathbf {x},\mathbf {z}) = \Vert \mathbf {x}- \mathbf {z}\Vert\) as the Euclidean distance from \(\mathbf {x}\) to \(\mathbf {z}\). It suffices to show

Indeed, taking this as given, let \(\mathbf {x}_k = {\boldsymbol \rho }(k)\). Then the sequence \(\left(\mathbf {x}_k\right)_{k = 1}^\infty\) is a Cauchy sequence in \(\mathbb {R}^\Lambda\), so it must converge to some state \(\mathbf {y}\). Moreover, we actually know that \({\boldsymbol \rho }(t)\) converges to some state \(\mathbf {y}\) as \(t \rightarrow \infty\): For any \(\epsilon \gt 0\) there is some N such that for all \(n \gt N\), we know \(d(\mathbf {x}_n, \mathbf {y}) \lt \epsilon /2\) and there is some M such that for all \(t, s \gt M\), we know \(d({\boldsymbol \rho }(t), {\boldsymbol \rho }(s)) \lt \epsilon /2\). In particular, taking m to be an integer larger than N and M, we see that \(\begin{equation*} d({\boldsymbol \rho }(t), \mathbf {y}) \le d({\boldsymbol \rho }(t), \mathbf {x}_m) + d(\mathbf {x}_m, \mathbf {y}) \lt \epsilon \end{equation*}\) for any \(t \gt M\). This shows that, taking (*) for granted, \({\boldsymbol \rho }\) must converge to \(\mathbf {y}\).

Let us now prove (*). Proceed by contradiction, so suppose that there is some \(\epsilon \gt 0\) such that for all \(M \in \mathbb {R}_{\ge 0}\), there exists some \(t \gt s \gt M\) such that \(d({\boldsymbol \rho }(t), {\boldsymbol \rho }(s)) \ge \epsilon\). Then take \(M_1 = 0\), and label the points we get by this assumption \(t_1, s_1\). Then for any \(n \gt 1\), take \(M_n = t_{n-1}\) and label the next pair of points \(t_n, s_n\). Then \(\begin{align*} \int _0^\infty \left\Vert \frac{d{\boldsymbol \rho }(t)}{dt}\right\Vert dt &\ge \sum _{i = 1}^\infty \int _{s_i}^{t_i}\left\Vert \frac{d{\boldsymbol \rho }(t)}{dt}\right\Vert dt \\ &\ge \sum _{i = 1}^\infty d({\boldsymbol \rho }(t_i), {\boldsymbol \rho }(s_i)) \\ &\ge \sum _{i = 1}^\infty \epsilon = \infty , \end{align*}\) where the second inequality uses the absolute continuity of \({\boldsymbol \rho }\). This gives a contradiction, establishing (*), which proves that \({\boldsymbol \rho }(t)\) converges to some state \(\mathbf {y}\) as \(t \rightarrow \infty\).

Now let us assume that \(\mathbf {f}\) is fair and establish that \(\mathbf {y}\) is a static state. If not, then there is some reaction \(\alpha\) applicable at \(\mathbf {y}\), so \(\mathbf {H}_\alpha (\mathbf {y}) = C \gt 0\). By the continuity of \(\mathbf {H}_\alpha\), there is some \(\delta\) so \(\mathbf {H}_\alpha (\mathbf {x}) \gt C/2\) for all \(\mathbf {x}\in \mathbb {R}_{\ge 0}^\Lambda\) with \(\Vert \mathbf {x}- \mathbf {y}\Vert \lt \delta\). Since \({\boldsymbol \rho }\) converges to \(\mathbf {y}\), this implies that \(\mathbf {H}_\alpha ({\boldsymbol \rho }(t)) \gt C/2\) for all t greater than some \(t_0\).

For any measurable \(U \subseteq \mathbb {R}_{\ge 0}\), let \(\mu (U) \in \mathbb {R}_{\ge 0}\cup \lbrace \infty \rbrace\) denote the measure of U. Letting \(T_\alpha\) be the subset of \(\mathbb {R}_{\ge 0}\) as in Definition 4.1 where \(\mathbf {f}_\alpha (t) \ge \mathbf {H}_\alpha ({\boldsymbol \rho }(t))\) for \(t \in T_\alpha\) and \(\mu (T_\alpha) = \infty\), we see that \(\mathbf {f}_\alpha (t) \gt C/2\) for every \(t \in T_\alpha \cap [t_0, \infty)\). Also note that \(\mu (T_\alpha \cap [t_0, \infty)) = \infty\). By the contrapositive direction of Lemma 4.7, we see that this implies that \(\Vert \frac{d{\boldsymbol \rho }(t)}{dt}\Vert \gt C/2K_\alpha\) for every \(t \in T_\alpha \cap [t_0, \infty)\). As a result, \(\begin{equation*} \int _0^\infty \left\Vert \frac{d{\boldsymbol \rho }(t)}{dt}\right\Vert dt \ge \int _{T_\alpha \cap [t_0, \infty)} \left\Vert \frac{d{\boldsymbol \rho }(t)}{dt}\right\Vert dt \ge \mu (T_\alpha \cap [t_0, \infty))\frac{C}{2K_\alpha } = \infty , \end{equation*}\) contradicting the finite length of \({\boldsymbol \rho }\).□

By applying Lemma 4.8 to a feedforward CRN that also stably computes a function, we obtain the following result, which states that such a CRN will reach the correct output under any fair rate schedule, from any state reachable from \(\mathbf {x}\). This is almost a converse to lemma 4.4; however, note that unlike in lemma 4.4, the CRC is required to be feedforward.

We point out that the feedforward property is not a necessary condition for stable computation to coincide with fair computation. For example, consider the non-feedforward CRC with reactions \(X + X \mathop {\rightarrow }\limits R + Y\) and \(R + R \mathop {\rightarrow }\limits X\). This CRC stably computes \(f(x) = (2/3) x\) [52], and it can also be shown that it fairly computes f. The key property it shares with feedforward CRCs is Lemma 4.7: A large flux through its reactions implies a large change in state.

Since mass-action yields a valid rate schedule (Lemma 2.24) that is fair (Lemma 4.2) and defined at all times for feedforward CRNs (Lemma 4.10), the following corollary is immediate. The corollary states that for a CRC stably computing a function f, the CRC will also converge to the correct output under mass-action kinetics, no matter the positive rate constants, and even if an adversary can first “steer” the CRC to some reachable state before letting mass-action kinetics take over.

5.1 Dual-rail Representations

The direct concentration-to-value mapping articulated in definition 3.2 is a straightforward way to represent non-negative input and output values. However, there are two reasons why an alternative, albeit more complex, encoding may be preferred. First, since concentrations cannot be negative, computable functions are restricted to non-negative domain and range. Second, as explained below, the direct output encoding frustrates the composition of smaller CRC modules into larger CRCs. The dual-rail representation we introduce in this section is a natural way to solve both problems.

A natural way to represent a (possibly negative) real value in chemistry is to encode it as the difference of two concentrations. Formally, let \(f:\mathbb {R}^k \rightarrow \mathbb {R}\) be a function. A function \(\hat{f}:\mathbb {R}_{\ge 0}^{2k} \rightarrow \mathbb {R}_{\ge 0}^2\) is a dual-rail representation of f if, for all \(\mathbf {x}^+,\mathbf {x}^- \in \mathbb {R}_{\ge 0}^k\), if \((y^+,y^-) = \hat{f}(\mathbf {x}^+,\mathbf {x}^-)\), then \(f(\mathbf {x}^+ - \mathbf {x}^-) = y^+ - y^-\). In other words, \(\hat{f}\) represents f as the difference of its two outputs \(y^+\) and \(y^-\), and it works for any input pair \((\mathbf {x}^+,\mathbf {x}^-)\) whose difference is the input value to f. We can define a CRC to stably compute such a function in the same manner as in Section 3, but having 2k input species \(\Sigma = \lbrace X_1^+,X_1^-,X_2^+,X_2^-,\ldots ,X_k^+,X_k^-\rbrace\) and two output species \(\Gamma = \lbrace Y^+,Y^-\rbrace\).

This definition implies that, for all \(\mathbf {x}= (\mathbf {x}^+, \mathbf {x}^-) \in \mathbb {R}_{\ge 0}^{2k}\), for all \(\mathbf {c}\) such that \(\mathbf {x}\rightsquigarrow \mathbf {c}\), there exists an output-stable state \(\mathbf {o}\) such that \(\mathbf {c}\rightsquigarrow \mathbf {o}\) and \(\mathbf {o}(Y^+) - \mathbf {o}(Y^-) = f(\mathbf {x}^+ - \mathbf {x}^-)\). Note that a single function has an infinite number of dual-rail representations; we require only that a CRC exists to compute one of them to say that the function is stably dual-computable by a CRC.

Besides making negative values chemically representable, we will see that the dual-rail representation plays a key role in allowing the composition of smaller CRC modules into a larger CRC. A key concept to enable such composition is output-obliviousness:

To recognize the problem with composition of the direct representation (definition 3.2), define the composition of two CRCs as the CRC that has the union of their reactions, relabeling the output species of the upstream CRN to be the input species of the downstream one [17, Definition 16]. Then with the direct output representation, output-obliviousness is necessary for composability but provably restricts computational power [17]: (1) Composing two stably computing CRCs stably computes the function composition if and only if the upstream CRC is output-oblivious (except in trivial cases). Intuitively, the downstream CRC can interfere with the upstream computation by prematurely consuming the output. (2) An output-oblivious CRC can only stably compute “superadditive” functions. For example, any CRC stably computing the function \(f(x_1,x_2) = x_1 - x_2\) must necessarily be able to consume its output species, since a state with more than the desired amount of output is reachable by additivity (e.g., state with \(x_1\) amount of Y).

Our results imply that the dual-rail representation allows composition without sacrificing computational power. In particular, our dual-rail constructions are all output-oblivious and thus composable by concatenation, and our negative results apply to dual-rail CRCs whether or not they are output-oblivious. To see roughly why the dual-rail representation helps with composition, consider the \(f(x_1,x_2) = x_1 - x_2\) function above. We can now compute this function with an output-oblivious (and therefore composable) CRC by producing \(Y^-\) to decrease the value of the output, without consuming any output species.

Recall fair computation (Definition 4.3). We define dual-rail fair computation analogously.

Since dual-rail computation is defined by a CRC that directly computes a dual-rail representation function, Lemmas 4.4 and 4.9 imply the analogous results for dual-rail:

5.2 Statement of Main Results

Below, we summarize our results about the computational power of stable computation, which we prove in subsequent subsections. First, we formally define the relevant classes of functions that will correspond to direct and dual-rail stable computation:

We note that rational linearity has the equivalent characterization that f is linear and maps rational inputs \(\mathbf {x} \in \mathbb {Q}^n\) to rational outputs.

Definition 5.8.

A function \(f: \mathbb {R}_{\ge 0}^k \rightarrow \mathbb {R}_{\ge 0}\) is positive-continuous if, for all \(U \subseteq \lbrace 1,\ldots ,k\rbrace\), f is continuous on the domain

In other words, f is continuous on any subset \(D \subset \mathbb {R}_{\ge 0}^k\) that does not have any coordinate \(i \in \lbrace 1,\ldots ,k\rbrace\) that takes both zero and positive values in D.

The following theorems are the main results of this article, exactly characterizing the functions stably computable with direct and dual-rail representation of inputs and outputs, and showing the equivalence between fair computation and stable computation. Furthermore, although non-feedforward CRCs do exist to compute functions in this class, the theorem shows that feedforward CRCs suffice to compute all such functions.

The following is the dual-rail analog of Theorem 5.9.

5.3 Positive Result: Continuous Piecewise Rational Linear Functions are Dual-rail Computable

Definition 5.7 does not stipulate how complex the “boundaries” between the linear pieces of a piecewise rational linear function can be. The boundaries can even be irrational in some sense, e.g., the function \(f(x_1,x_2) = 0\) if \(x_1 \gt \sqrt {2} \cdot x_2\) and \(f(x_1,x_2) = x_1+x_2\) otherwise. However, if we additionally require that f be continuous, then the following theorem of Ovchinnikov [42, Theorem 2.1] states that f has a particularly clean form, conducive to computation by CRCs.

Note that as a special case, the above result applies when f is continuous piecewise rational linear. The above theorem as stated slightly generalizes the result due to Ovchinnikov [42] (by not requiring D to be closed), although the proof technique is essentially the same. For completeness, we provide the proof in Appendix B.

We use the theorem above to dual-compute continuous piecewise rational linear functions by composing CRC modules for rational linear functions, min, and max. These modules are developed in the following three lemmas:

Proof.

To stably compute a dual-rail representation of max, observe that it is equivalent to computing the min function with the roles of the “plus” and “minus” species reversed (which negates the value represented in dual-rail), because \(\max (x_1,x_2) = - \min (-x_1,-x_2)\). In other words, use the reactions

□

Note that Lemma 5.15 applies for arbitrary convex domains \(D \subseteq \mathbb {R}^k\), which will be useful in proving Lemma 5.16, where we take D to be a strict convex subset of \(\mathbb {R}_{\ge 0}^k\) in which no coordinate takes both 0 and positive values in D. Since \(\mathbb {R}^k\) is itself convex, it also establishes the “(3) implies (4)” implication in the proof of Theorem 5.10.

5.4 Positive Result: Positive-continuous Piecewise Rational Linear Functions Are Directly Computable

The following lemma is the direct stable computation analog of Lemma 5.15. Intuitively, it is proven by using Lemma 5.15 to stably dual-compute \(2^k\) different continuous piecewise rational linear functions in parallel, one for each possible choice of which input species \(X_1,\ldots ,X_k\) are 0. A separate computation determines which inputs are positive and selects the appropriate output. Note that positive inputs may be discovered “piecemeal,” so the system must be robust to a continual updating of the decision. However, there is a monotonicity to this process that will make consistent updating possible: Once an input species is discovered to be present, we know for sure it is. (Whereas a species that appears to be absent simply may have not yet reacted.)

Proof.

The CRC will have input species \(X_1,\ldots ,X_k\) and output species \(Y^+\). (While it will be helpful to think of a \(Y^+\) and \(Y^-\) species, and during the computation the output will be encoded in their difference, the output of the CRC is only the \(Y^+\) species as per direct computability.)

Let \(f:\mathbb {R}_{\ge 0}^k \rightarrow \mathbb {R}_{\ge 0}\) be a positive-continuous piecewise linear function. Since it is positive-continuous, there exist \(2^k\) domains

one for each subset \(U \subseteq \lbrace 1,\ldots ,k\rbrace\), such that \(f \upharpoonright D_U\) is continuous. Define \(f_U = f \upharpoonright D_U\). Since \(D_U\) is convex, by Lemma 5.15 there is a CRC \(\mathcal {C}_U\) computing a dual-rail representation \(\hat{f}_U:\mathbb {R}_{\ge 0}^k\times \mathbb {R}_{\ge 0}^k\rightarrow \mathbb {R}\times \mathbb {R}\) of \(f_U\). By letting the initial concentration of the “minus” version of the ith input species \(X_i^-\) be 0, we convert \(\mathcal {C}_U\) into a CRC that directly computes an output dual-rail representation of \(f_U\).

The intuition of the proof is as follows: The case \(U=\emptyset\) is trivial, as we will have no reactions of the form \(\emptyset \rightarrow A\) for any species A, so if no species are initially present, then no species (including \(Y^+\)) will ever be produced; this is correct, since any linear function f obeys \(f(\mathbf {0})=0\). For each non-empty U, we compute \(f_U\) independently in parallel by CRC \(\mathcal {C}_{U}\), modifying each reaction producing \(Y^+\) to produce an equivalent amount of species \(Y_U\), which is specific to U. For each such U there are inactive and active “indicator” species \(J_U\) and \(I_U\). In parallel, there are reactions that will activate indicator species \(I_U\) (i.e., convert \(J_U\) to \(I_U\)) if and only if all species \(X_i\) are present initially for each \(i \in U\). These \(I_U\) species will then counteract the effect of any CRC computing \(f_{U^{\prime }}\) for \(U^{\prime } \subsetneq U\) by catalytically converting all \(Y_{U^{\prime }}^+\) to \(Y^-\) and all \(Y_{U^{\prime }}^-\) to \(Y^+\). If U is the complete set of indices of non-zero inputs, then only CRCs computing \(f_{U^{\prime }}\) for subsets \(U^{\prime } \subsetneq U\) have produced any amount of \(Y^+\), so eventually all of these will be counteracted by \(I_U\).

Formally, construct the CRC as follows. Let \(l = 2^{k-1}-1\). For each \(i \in \lbrace 1,\ldots ,k\rbrace\), add the reaction \(X_i \rightarrow I_{\lbrace i\rbrace } + J_{U_1} + X_i^{U_1} + J_{U_2} + X_i^{U_2} + \cdots + J_{U_{l}} + X_i^{U_{l}},\) where \(U_1,U_2,\ldots ,U_{l}\) are all subsets of \(\lbrace 1,\ldots ,k\rbrace\) that are strict supersets of \(\lbrace i \rbrace\). The extra versions \(X_i^{U_1}, \ldots , X_i^{U_l}\) of \(X_i\) are used as inputs to the parallel computation of each \(f_U\). We generate the inactive indicator species from the input species in this manner, because the CRC is not allowed to start with anything other than the input.

The indicator species are activated as follows: For each nonempty \(U,U^{\prime } \subseteq \lbrace 1,\ldots ,k\rbrace\) such that \(U \ne U^{\prime }\), add the reaction \(I_U + I_{U^{\prime }} + J_{U \cup U^{\prime }} \rightarrow I_U + I_{U^{\prime }} + I_{U \cup U^{\prime }}.\)

For each nonempty \(U \subseteq \lbrace 1,\ldots ,k\rbrace\), let \(\mathcal {C}_U\) be the CRC computing an output dual-rail representation of \(f_U\) (i.e., dual rail on the output). Modify \(\mathcal {C}_U\) as follows: Rename the output species of \(\mathcal {C}_U\) to \(Y^+\) and \(Y^-\), i.e., all parallel CRCs share the same output species. For each reaction producing the output species \(Y^+\), add the product \(Y^+_U\) (which is a species-specific to \(\mathcal {C}_U\)) with the same net stoichimetry. Similarly, for each reaction producing the output species \(Y^-\), add the product \(Y^-_U\) with the same net stoichimetry. For instance, replace the reaction \(A + B \rightarrow Y^+\) by the reaction \(A + B \rightarrow Y^+ + Y^+_U\), and replace the reaction \(A + Y^+ \rightarrow B + 4 Y^+\) by the reaction \(A + Y^+ \rightarrow B + 4 Y^+ + 3 Y^+_U\). Therefore, the eventual amount of \(Y_U^+\) is equal to the total amount of \(Y^+\) produced by \(\mathcal {C}_U\), and similarly for \(Y_U^-\) and \(Y^-\). For each \(U^{\prime } \subset U\), add the reactions \(I_U + Y^+_{U^{\prime }} \rightarrow I_U + Y^-\), \(I_U + Y^-_{U^{\prime }} \rightarrow I_U + Y^+.\) Also, for each reaction in \(\mathcal {C}_U\), add \(I_U\) as a catalyst. This ensures that \(\mathcal {C}_U\) cannot execute any reactions (and therefore cannot produce any amount of \(Y^+\) or \(Y^-\)) unless all species \(X_i\) for \(i \in U\) are present.

We observe that the dual-rail CRC described above (with output species \(Y^+\) and \(Y^-\)) is output-oblivious, as it involves the output species only as products of reactions. Further, the output value is non-negative for any input. Thus, we can convert the dual-rail representation to the direct one (definition 3.2) with a single output species \(Y^+\) by adding the reaction \(Y^+ + Y^- \rightarrow \emptyset\).

To complete the proof, since the CRC of Lemma 5.15 is feedforward, we can confirm by inspection our modifications preserve the feedforward property as well. In particular, one should order species within the CRN dual-computing \(f_U\) before any species for supersets of U and after any species for subsets of U.□

5.5 Negative Result: Directly Computable Functions Are Positive-continuous Piecewise Rational Linear

5.5.2 Linearity Restricted to Inputs Draining a Siphon.

This section aims to prove that the function computed by \(\mathcal {C}\), when restricted to inputs that can drain a particular stabilizing siphon, is linear. This is the first step establishing the “piecewise linear” portion of the “positive-continuous piecewise rational linear” claims in Theorem 5.9.

Recall that \(\Psi\), defined in Definition 2.4, is the space of all prepaths (i.e., the vector space in which paths live), and \(\Gamma _\infty\), defined in Definition 2.5, is the space of all paths (allowing infinitely many line segments). We now define a map \(\mathbf {o}\) that intuitively sends a path to the final state that it reaches.

Lemma 5.19.

The map \(\mathbf {o}: \Psi \rightarrow \mathbb {R}^\Lambda\) sending \(\begin{equation*} {\boldsymbol \gamma }\mapsto \lim _{n \rightarrow \infty } \mathbf {x}_n({\boldsymbol \gamma }) \end{equation*}\) is linear (recall \(\mathbf {x}_n({\boldsymbol \gamma })\), defined in Definition 2.4, is the state reached after traversing n segments along \({\boldsymbol \gamma }\)).

Proof.

To check that \(\mathbf {o}\) is a linear function, note that for any \({\boldsymbol \gamma }_0, {\boldsymbol \gamma }_1 \in \Psi\) and \(\lambda \in \mathbb {R}\), we have

□

Definition 5.20.

Let \(\Omega\) be a stabilizing siphon. Define \(\Gamma (\Omega)\) to consist of paths \({\boldsymbol \gamma }\in \Gamma _\infty\) such that \(\mathbf {o}({\boldsymbol \gamma })\upharpoonright \Omega = \mathbf {0}\).

These are the paths that converge to a state where a given stabilizing siphon is drained.

Lemma 5.21.

\(\Gamma (\Omega)\) is convex for each stabilizing siphon \(\Omega\).

Proof.

Suppose that \({\boldsymbol \gamma }_0\) and \({\boldsymbol \gamma }_1\) are in \(\Gamma (\Omega)\) and let \({\boldsymbol \gamma }_\lambda = (1-\lambda){\boldsymbol \gamma }_0 + \lambda {\boldsymbol \gamma }_1\). By Lemma 2.7, we know that \({\boldsymbol \gamma }_\lambda\) is in \(\Gamma _\infty\). Moreover, \(\begin{equation*} \mathbf {o}({\boldsymbol \gamma }_\lambda) = (1 - \lambda)\mathbf {o}({\boldsymbol \gamma }_0) + \lambda \mathbf {o}({\boldsymbol \gamma }_1). \end{equation*}\) Since \(\Omega\) is drained at both \(\mathbf {o}({\boldsymbol \gamma }_0)\) and \(\mathbf {o}({\boldsymbol \gamma }_1)\), we conclude that it must also be drained at \(\mathbf {o}({\boldsymbol \gamma }_\lambda)\), so \({\boldsymbol \gamma }_\lambda \in \Gamma (\Omega)\).□

Definition 5.22.

Let

denote those input states from which the siphon \(\Omega\) is drainable.

Lemma 5.23.

Let \(f:\mathbb {R}_{\ge 0}^k \rightarrow \mathbb {R}_{\ge 0}\) be stably computed by a CRC \(\mathcal {C}=(\Lambda ,R,\Sigma ,\lbrace Y\rbrace)\). Let \(\Omega\) be a stabilizing siphon. Then f restricted to \(\Sigma (\Omega)\) is a linear function.

Proof.

Recall \(\Gamma (\Omega)\) from Definition 5.20, the map \(\mathbf {x}_0\) from Definition 2.4, and the map \(\mathbf {o}\) from Lemma 5.19. First project \(\Gamma (\Omega)\) to \(\mathbb {R}_{\ge 0}^\Lambda \times \mathbb {R}_{\ge 0}^\Lambda\) by the map \({\boldsymbol \gamma }\mapsto (\mathbf {x}_0({\boldsymbol \gamma }), \mathbf {o}({\boldsymbol \gamma }))\). Let \(G \subseteq \mathbb {R}^{k+1}\) be the further projection to the \((k+1)\)-dimensional subspace corresponding only to the input species \(X_1,\ldots ,X_k\) and output species Y. G is the graph of the function \(y = f(\mathbf {x})\) restricted to inputs \(\mathbf {x}\in \Sigma (\Omega)\). Since G is the image of a convex set under a linear transformation, it is also convex. We claim that G must be a subset of a k-dimensional hyperplane.

For the sake of contradiction, suppose not. Then there are \(k+1\) non-coplanar points in G. Since G is convex, it contains the entire \((k+1)\)-dimensional convex hull H of these points. Since H is a \((k+1)\)-dimensional convex polytope, it contains two different values of y corresponding to the same value of \(\mathbf {x}\), contradicting the fact that only a single y value exists in all output-stable states reachable from \(\mathbf {x}\). This establishes the claim that G must be a subset of a k-dimensional hyperplane.

Since the graph of f is a subset of a k-dimensional hyperplane, f is an affine function. Since there are no reactions of the form \(\emptyset \rightarrow \ldots\), Y cannot be produced from the initial state \(\mathbf {x}=\mathbf {0}\) (nor can any other species), so \(f(\mathbf {0})=0\). Therefore, this hyperplane passes through the origin, so it defines a linear function.□

In Lemma 5.23, the reason that we restrict attention to a single output siphon \(\Omega\) is that if different output siphons are drained, then different linear functions may be computed by the CRC. For example, \(X_1+X_2\rightarrow Y\) computes \(f(x_1,x_2) = x_1\) if siphon \(\lbrace X_2\rbrace\) is drained and \(f(x_1,x_2) = x_2\) if siphon \(\lbrace X_1\rbrace\) is drained. If a CRC stably computes, then an output-stable state is reachable from any input state. Thus, by Lemma 5.18, from every input state some stabilizing siphon is drainable, and the following corollary is immediate:

Corollary 5.24.

Let \(f:\mathbb {R}_{\ge 0}^k \rightarrow \mathbb {R}_{\ge 0}\) be stably computed by a CRC. Then f is piecewise linear.

5.5.3 Positive-continuity.

Ideally, to prove that the function stably computed by a CRC is positive-continuous, we would like to prove the following: For any stabilizing siphon \(\Omega\), the set \(\Sigma (\Omega)\) of input states that can drain \(\Omega\) is closed relative to the positive orthant. If that were true, then we could use a fundamental topological result that if a function is piecewise continuous with finitely many pieces (e.g., piecewise linear), and if the domain defining each piece is closed (with agreement between pieces on intersecting domains), then the whole function is continuous. However, the above statement is not true in general. Consider the following counterexample: \(\begin{align*} X_1 &\rightarrow C \\ X_1 + X_2 + C &\rightarrow C + Y. \end{align*}\) If initially \(\mathbf {i}(X_1) \gt \mathbf {i}(X_2)\), then the stabilizing siphon \(\lbrace X_2\rbrace\) is drainable, by producing \((\mathbf {i}(X_1)\, -\, \mathbf {i}(X_2)) / 2\) of C via the first reaction (leaving an excess of \(X_1\) over \(X_2\) still), then running the second reaction until \(X_2\) is gone to produce Y. Because \(X_2\) can only be consumed if C is produced, which requires consuming a positive amount of \(X_1\), the set of inputs from which \(\lbrace X_2\rbrace\) can be drained is the non-closed set \(\lbrace \mathbf {i}\mid \mathbf {i}(X_1) \gt \mathbf {i}(X_2)\rbrace\). Note that the above CRC does not stably compute anything because, starting from a state with \(\mathbf {i}(X_1) \gt \mathbf {i}(X_2)\), the first reaction could run until \(X_2\) exceeds \(X_1\) before starting the second reaction, which would imply that the amount of Y produced depends on how much reaction 1 happens. It is still unclear whether such counterexamples exist for CRCs stably computing some function that is not identically 0.

Instead of relying on \(\Sigma (\Omega)\) being closed, we must make a more careful argument. In lieu of working directly with the sets \(\Sigma (\Omega)\), we consider “shifted” sets \(\tilde{\Sigma }_{(\mathbf {y}, \mathbf {z})}(\Omega)\) (see Definition 5.28 below). Each \(\tilde{\Sigma }_{(\mathbf {y}, \mathbf {z})}(\Omega)\) is (possibly strictly) contained in the original \(\Sigma (\Omega)\), but they still cover the set of inputs. Crucially, we are able to show that the shifted sets \(\tilde{\Sigma }_{(\mathbf {y}, \mathbf {z})}(\Omega)\) are closed, allowing us to apply the argument at the start of this section to prove that every function stably computed by a CRC is positive-continuous.

Definition 5.25.

Let

denote those states from which siphon \(\Omega\) is drainable via a single straight line segment.

Lemma 5.26.

Let \(\Omega\) be a siphon. Let \(\mathbf {a}_1,\mathbf {a}_2,\ldots \in X(\Omega)\) be a convergent sequence of states, where \(\mathbf {a}=\lim _{i\rightarrow \infty } \mathbf {a}_i\). Suppose \([\mathbf {a}] = \Lambda\). Then \(\mathbf {a}\in X(\Omega).\)

Proof.

Consider the set P of \({\boldsymbol \gamma }= (\mathbf {x}_0, \mathbf {u}) \in \mathbb {R}^\Lambda _{\ge 0} \times \mathbb {R}^R_{\ge 0}\) such that \(\mathbf {o}({\boldsymbol \gamma }) \in \mathbb {R}^\Lambda _{\ge 0}\) and \(\mathbf {o}({\boldsymbol \gamma }) \upharpoonright \Omega = \mathbf {0}\). Note that reactions occurring with positive flux in \(\mathbf {u}\) might not be applicable at \(\mathbf {x}_0\). P is cut out by a system of non-strict linear inequalities (in other words, it is a polyhedron). By Reference [53], \(\mathbf {x}_0(P)\) is also a polyhedron, and is in particular closed.

Note that \(X(\Omega) \subseteq \mathbf {x}_0(P)\) and \(\mathbf {x}_0(P) \cap \mathbb {R}_{\gt 0}^\Lambda = X(\Omega) \cap \mathbb {R}_{\gt 0}^\Lambda\). The first relation follows, since if \(\mathbf {x} \in X(\Omega)\), then the straight-line path draining \(\Omega\) produces a \({\boldsymbol \gamma }\in P\) such that \(\mathbf {x}_0({\boldsymbol \gamma }) = \mathbf {x}\). The second relation holds, since if every species is present in \(\mathbf {x}_0\) then every reaction is applicable at \(\mathbf {x}_0\), so any of the points \({\boldsymbol \gamma }\in P\) that project to \(\mathbf {x}_0\) are valid paths that drain \(\Omega\).

Since \(\mathbf {a}= \lim _{i\rightarrow \infty } \mathbf {a}_i\), we have that for all but finitely many i, \(\mathbf {a}_i \in X(\Omega) \cap \mathbb {R}_{\gt 0}^\Lambda\). As a result, these \(\mathbf {a}_i\) are in \(\mathbf {x}_0(P)\), so \(\mathbf {a}\) in in \(\mathbf {x}_0(P)\), too, since the set is closed. Since \(\mathbf {a}\in \mathbf {x}_0(P) \cap \mathbb {R}_{\gt 0}^\Lambda\), we conclude that \(\mathbf {a}\in X(\Omega).\)□

Note that the hypothesis \([\mathbf {a}]=\Lambda\) is necessary. Otherwise, consider the reactions \(A\rightarrow C\), \(A+B\rightarrow \emptyset\), and \(F+C\rightarrow C\), with \(\mathbf {a}_i(C)=0\), \(\mathbf {a}_i(F)=1\), \(\mathbf {a}_i(B)=1\), and \(\mathbf {a}_i(A)\) approaching 1 from above as \(i \rightarrow \infty\) (whence \(C \not\in [\mathbf {a}]\)). Then the siphon \(\Omega =\lbrace A,B,F\rbrace\) is drainable from each \(\mathbf {a}_i\) by running \(A \rightarrow C\) until A and B have the same concentration, then running the other two reactions to completion. However, \(\mathbf {a}(A)=\mathbf {a}(B)\), so running any amount of reaction \(A \rightarrow C\) prevents reaction \(A+B\rightarrow \emptyset\) from draining B. Therefore, \(\mathbf {a}\not\in X(\Omega)\) but \(\mathbf {a}_i \in X(\Omega)\) for all i.

Definition 5.27.

A pair \((\mathbf {y}, \mathbf {z}) \in \mathbb {R}_{\ge 0}^\Lambda \times \mathbb {R}_{\ge 0}^\Lambda\) of states is a full input pair if \([\mathbf {y}] = \Sigma\), \([\mathbf {z}] = \Lambda\), and \(\mathbf {y}\rightsquigarrow \mathbf {z}\).

Definition 5.28.

If \((\mathbf {y}, \mathbf {z})\) is a full input pair and \(\Omega\) is a stabilizing siphon, then define

Intuitively, in the CRC at the beginning of the section, the obstacle to the set \(\Sigma (\Omega)\) being closed for the siphon \(\Omega = \left\lbrace X_2\right\rbrace\) was that not all species were present initially: To drain \(\Omega\), you first need to produce some (arbitrarily small) amount of C, which leads to the requirement \(\mathbf {i}(X_1) \gt \mathbf {i}(X_2)\). To fix this problem, \(\tilde{\Sigma }_{(\mathbf {y}, \mathbf {z})}(\Omega)\) considers only the states that can drain \(\Omega\) after having been “perturbed” into a state where all species are present, where the full input pair \((\mathbf {y}, \mathbf {z})\) specifies how to perform this perturbation.

To see how this works in the example CRC, let \((\mathbf {y}, \mathbf {z})\) be the full input pair where \(\mathbf {y}= 1X_1, 1X_2, 0C, 0Y\) and \(\mathbf {z}= \left\lbrace .25X_1, .75X_2, .5C, .25Y\right\rbrace\). Then for an input state \(\mathbf {x}= \left\lbrace aX_1, bX_2, 0C, 0Y\right\rbrace\), \(\lambda \mathbf {y}\lt \mathbf {x}\) when \(\lambda \lt \min (a,b)\). As a result, \(\mathbf {x}\in \tilde{\Sigma }_{(\mathbf {y}, \mathbf {z})}(\Omega)\) if \(\mathbf {x}- \lambda \mathbf {y}+ \lambda \mathbf {z}= (a - .75\lambda)X_1, (b - .25\lambda)X_2, .5\lambda C, .25\lambda Y\) can drain \(X_2\) for all \(0 \lt \lambda \lt \min (a,b)\). This happens when \(b - .25\lambda \le a - .75\lambda\), so \(b + .5\lambda \le a\). Taking the limit as \(\lambda \rightarrow 0\), we see that \(b \le a\), so \(\min (a,b) = b\), and taking the limit as \(\lambda \rightarrow \min (a,b) = b\), we see that \(1.5 b \le a\). Since \(1.5b \le a\) is also a sufficient condition for \(\mathbf {x}\) to be in \(\tilde{\Sigma }_{(\mathbf {y}, \mathbf {z})}(\Omega)\), we see that \(\tilde{\Sigma }_{(\mathbf {y}, \mathbf {z})}(\Omega)\) is closed and contained in \(\Sigma (\Omega)\). The next two lemmas show that these properties hold in general.

Lemma 5.29.

For any full input pair \((\mathbf {y}, \mathbf {z})\) and any stabilizing siphon \(\Omega\), \(\tilde{\Sigma }_{(\mathbf {y}, \mathbf {z})}(\Omega)\) is closed relative to \(\mathbb {R}_{\gt 0}^\Sigma\).

Proof.

Let \(\mathbf {x}\) be a state such that \([\mathbf {x}] = \Sigma\) and let \(\lbrace \mathbf {x}_i\rbrace\) be a sequence such that \(\mathbf {x}= \lim _{i \rightarrow \infty } \mathbf {x}_i\) and \(\mathbf {x}_i \in \tilde{\Sigma }_{(\mathbf {y}, \mathbf {z})}(\Omega)\). For any \(\lambda \gt 0\) such that \(\lambda \mathbf {y}\lt \mathbf {x}\), by throwing out finitely many terms in the sequence \(\lbrace \mathbf {x}_i\rbrace\), we can guarantee that \(\lambda \mathbf {y}\lt \mathbf {x}_i\) for all i, too. Since \(\mathbf {x}_i \in \tilde{\Sigma }_{(\mathbf {y}, \mathbf {z})}(\Omega)\), we know that \(\mathbf {x}_i - \lambda \mathbf {y}+ \lambda \mathbf {z}\in X(\Omega)\) for all i. Since \(X(\Omega)\) is closed, we see that \(\begin{equation*} \mathbf {x}-\lambda \mathbf {y}+ \lambda \mathbf {z}= \left(\lim _{i \rightarrow \infty } \mathbf {x}_i\right) -\lambda \mathbf {y}+ \lambda \mathbf {z}= \lim _{i \rightarrow \infty } (\mathbf {x}_i -\lambda \mathbf {y}+ \lambda \mathbf {z}) \end{equation*}\) is also in \(X(\Omega)\). Since this is true for every \(\lambda\) such that \(\lambda \mathbf {y}\lt \mathbf {x}\), we conclude that \(\mathbf {x}\in \tilde{\Sigma }_{\mathbf {y}, \mathbf {z}}(\Omega)\). Since this is true for any \(\mathbf {x}\) in \(\mathbb {R}_{\gt 0}^\Sigma\), we conclude that \(\tilde{\Sigma }_{\mathbf {y}, \mathbf {z}}(\Omega)\) is closed relative to \(\mathbb {R}_{\gt 0}^\Sigma\).□

The following lemma is almost immediate from the definition. The only possible concern one might have is that an input state \(\mathbf {x}\) is contained in \(\tilde{\Sigma }_{(\mathbf {y}, \mathbf {z})}(\Omega)\) “vacuously”—in other words, that there simply does not exist a \(\lambda \gt 0\) such that \(\lambda \mathbf {y}\lt \mathbf {x}\).

Lemma 5.30.

For any full input pair \((\mathbf {y}, \mathbf {z})\) and any stabilizing siphon \(\Omega\), \(\tilde{\Sigma }_{(\mathbf {y}, \mathbf {z})}(\Omega) \subseteq \Sigma (\Omega)\).

Proof.

Let \(\mathbf {x}\) be in \(\tilde{\Sigma }_{(\mathbf {y}, \mathbf {z})}\). By definition, \([\mathbf {y}] = \Sigma\), so the following is a well-defined real number: \(\begin{equation*} \lambda _0 = \min _{S \in \Sigma } \left\lbrace \frac{\mathbf {x}(S)}{\mathbf {y}(S)}\right\rbrace . \end{equation*}\) In other words, \(\lambda _0\) is the number so \(\lambda \mathbf {y}\lt \mathbf {x}\) if and only if \(\lambda \lt \lambda _0\). Because \([\mathbf {x}] = \Sigma\), we know that \(\lambda _0 \gt 0\). Also \((\lambda _0/2) \mathbf {y}\rightsquigarrow (\lambda _0/2) \mathbf {z}\) because \((\mathbf {y}, \mathbf {z})\) is a full input pair. Then, by additivity of \(\rightsquigarrow\), \(\mathbf {x}\rightsquigarrow \mathbf {x}- (\lambda _0/2) \mathbf {y}+ (\lambda _0/2) \mathbf {z}\rightarrow ^1 \mathbf {o}\), where \(\Omega\) is drained at \(\mathbf {o}\). Because \(\rightsquigarrow\) is transitive (Corollary 2.16), we conclude that \(\mathbf {x}\rightsquigarrow \mathbf {o}\), so \(\mathbf {x}\in \Sigma (\Omega)\).□

The above two lemmas show that \(\tilde{\Sigma }_{(\mathbf {y}, \mathbf {z})}(\Omega)\) is topologically better behaved than \(\Sigma (\Omega)\), although possibly smaller. However, for these sets to be useful for analyzing the behavior of a CRC, we need to show that they are not “too small”—in particular, the next lemma shows that under the assumption that a CRC stably computes a function then for any fixed choice of a full input pair \((\mathbf {y}, \mathbf {z})\), the sets \(\tilde{\Sigma }_{(\mathbf {y}, \mathbf {z})}(\Omega)\) are still big enough to cover all of \(\mathbb {R}_{\gt 0}^\Sigma\).

Lemma 5.31.

If \(\mathcal {C}\) is a stably computing CRC, then for any fixed full input pair \((\mathbf {y}, \mathbf {z})\), as \(\Omega\) varies among all of the stabilizing siphons, the sets \(\tilde{\Sigma }_{(\mathbf {y}, \mathbf {z})}(\Omega)\) cover \(\mathbb {R}_{\gt 0}^\Sigma\).

Proof.

Let \(\mathbf {x}\) be an input state such that \([\mathbf {x}] = \Sigma\) and let \(\lambda _0\) be as in the proof of Lemma 5.30. We know that \(\lambda _0\mathbf {y}\rightsquigarrow \lambda _0 \mathbf {z}\) because \((\mathbf {y}, \mathbf {z})\) is a full input pair. By additivity of \(\rightsquigarrow\),

so, since \(\mathcal {C}\) is a stably-computing CRC there must be some output-stable state \(\mathbf {o}\) such that \(\mathbf {x}- \lambda _0 \mathbf {y}+ \lambda _0 \mathbf {z}\rightsquigarrow \mathbf {o}\). By Lemma 5.18, there is an stabilizing siphon \(\Omega\) so \(\Omega\) is drained at \(\mathbf {o}\).

For any \(\lambda \gt 0\) such that \(\lambda \mathbf {y}\lt \mathbf {x}\), we know that \(\lambda \lt \lambda _0\), so

Since \(\mathbf {x}- \lambda _0 \mathbf {y}+ \lambda _0 \mathbf {z}\rightsquigarrow \mathbf {o}\), we see that \(\mathbf {x}- \lambda \mathbf {y}+ \lambda \mathbf {z}\rightsquigarrow \mathbf {o}\), and, since \([\mathbf {z}] = \Lambda\), by Lemma 2.13, we conclude that \(\mathbf {x}- \lambda \mathbf {y}+ \lambda \mathbf {z}\rightarrow ^1 \mathbf {o}\). Since this is true for any \(\lambda \gt 0\) such that \(\lambda \mathbf {y}\lt \mathbf {x}\), we conclude that \(\mathbf {x}\in \tilde{\Sigma }_{(\mathbf {y}, \mathbf {z})}(\Omega)\). This shows that every \(\mathbf {x}\) with \([\mathbf {x}] = \Sigma\) is in \(\tilde{\Sigma }_{(\mathbf {y}, \mathbf {z})}(\Omega)\) for some \(\Omega\), as desired.□

We now use the above technical machinery to prove the following result, which is almost the full negative result for direct computation, but leaves out the constraint that f is rational linear. Rationality is shown in Section 5.5.4 below. Recall the positive-continuous functions from Definition 5.8.

Lemma 5.32.

Let \(f:\mathbb {R}_{\ge 0}^k \rightarrow \mathbb {R}_{\ge 0}\) be stably computed by a CRC. Then f is positive-continuous and piecewise linear.

Proof.

Piecewise linearity follows from Corollary 5.24. For positive continuity, we proceed as follows: Let \(U \subseteq \lbrace 1,\ldots ,k\rbrace\), let \(\mathbf {x}\in D_U\) (where \(D_U\) is as defined in Definition 5.8), and let \(\mathbf {x}_1\), \(\mathbf {x}_2\), \(\ldots \in D_U\) be an infinite sequence of points such that \(\lim _{i\rightarrow \infty } \mathbf {x}_i = \mathbf {x}\). It suffices to show that \(\lim _{i\rightarrow \infty } f(\mathbf {x}_i) = f(\mathbf {x})\)—i.e., that f is continuous on \(D_U\). We take \(\mathbf {x}_i\) and \(\mathbf {x}\) equivalently to represent an initial state of the CRC giving the concentrations of species in \(\Sigma =\lbrace X_1,\ldots ,X_k\rbrace\).

In analyzing the behavior of the CRC on states in \(D_U\), it will help us to consider the functionally equivalent CRC in which we remove species that are not producible from any state in \(D_U\). For the purposes of this proof, we consider this reduced CRC, and let \(\Lambda\) be the corresponding reduced set of species.

Let \((\mathbf {y}, \mathbf {z})\) be some full input pair. Then, as \(\Omega\) varies among the stabilizing siphons, \(\tilde{\Sigma }_{(\mathbf {y}, \mathbf {z})}(\Omega)\) gives a finite collection of closed sets covering \(\mathbb {R}_{\gt 0}^\Sigma\) by Lemmas 5.31 and 5.29. Since \(\tilde{\Sigma }_{(\mathbf {y}, \mathbf {z})}(\Omega) \subseteq \Sigma (\Omega)\) by Lemma 5.30 and, since f is linear (and therefore continuous) on \(\Sigma (\Omega)\) by Lemma 5.23, we see that f is continuous on each of the closed sets in this covering. By Reference [40], if a topological space is a union of finitely many closed sets and \(f_i\) are continuous function on each closed set that agree on overlaps, then they combine to give a continuous function. From this result, we conclude that f is continuous on \(D_U\), as desired.□

5.6 Negative Result: Dual-rail Computable Functions Are Continuous Piecewise Rational Linear

The following result, a dual-rail analog of Lemma 5.36, is not necessary for the proof of the main result of this section (Lemma 5.40), but it may be of independent interest.

The following is our main negative result for dual-rail CRC’s, a dual-rail analog of Lemma 5.38:

6.1 A Game-theoretic Formulation of Rate-independent Computation

In this section, we use our notion of a valid rate schedule (Definition 2.22) to propose a framework for studying rate-independent computation. Intuitively, we consider a model where the CRC experiences intermittent “shocks” where the system behaves erratically and the user of the CRC loses some amount of control of the rates of its reactions. We then say that the CRC rate-independently computes a function if, despite these shocks, the CRC still converges to the correct output.

In principle, one could consider different versions of the above model, depending on how much control the user of the CRC is expected to have over the kinetics of the system and how severe the shocks to the system are expected to be. Our goal in this section is to introduce a framework that is general enough to accommodate these various situations.

We will formalize the situation presented above by describing it as an infinite game, played by two players, which we will call the Demon and the Chemist. Here, the Chemist represents the user of the CRC, who is trying to perform some rate-independent computation, and the Demon represents (per tradition) the forces of nature, human failing, and so on, which create the intermittent shocks to the system. In this game, the Demon and the Chemist take turns building a valid rate schedule. The Chemist wins the game if the amount of the output species converges to the correct value as time goes to infinity. If the Chemist has a winning strategy for the game associated with a given CRC, then we say that the CRC rate-independently computes the desired function.

To model the different levels of control that the Chemist has over the system, and the different levels of severity of the shocks that the Demon can induce in the system, we can consider games where the Demon and the Chemist are only allowed to play moves from some restricted class of valid rate schedules. In this article, we have been interested in functions that are rate-independently computable in a very strong sense, i.e., where the shocks are arbitrarily severe—so the Demon is allowed to play any valid rate schedule (Strong Demon).

Once we have developed the game-theoretic model of rate-independent computation more explicitly, we will be able to use our main result Theorem 5.9 to deduce the following claim: The class of functions that one can compute rate-independently with a Strong Demon is effectively insensitive to the amount of control that the Chemist has over the system. In particular, the class of functions that are rate-independently computable by Strong Demon, Strong Chemist games, where both the Demon and the Chemist can play any valid rate law, is the same as any Strong Demon complexity class, where the Chemist is restricted to playing only a subset of valid rate laws.

To formalize the above discussion, we can use the notion of an infinite game [41].

Intuitively, one should imagine that two players, player I and player II, are taking turns playing moves from the set X. Together they form a sequence \((x_1, x_2, \ldots) \in X^\mathbb {N}\) and player I wins if this sequence is in A. In most games of interest, the set of moves that a player can make is restricted by the state of the game. The definition above captures this by using a payoff set such that a player will lose instantly if they play a move outside of some “legal” set of moves. For any previously played sequence, we can thus let moveset \(M(x_1, \ldots , x_n)\) be the legal set of moves for the two players (depending on whether n is odd or even) and assume that the payoff set A is defined relative to these sets of allowed moves accordingly (i.e., the first time that an illegal move is made, the other party wins). An infinite sequence of moves is consistent with the moveset if every move by both players is legal.

For our context, fix a CRC \(\mathcal {C}\) and an initial state \(\mathbf {x}\in \mathbb {R}_{\ge 0}^\Sigma\). The set of possible moves X is the set of all rate schedules that are zero at times \(t \gt 1\). Such moves (i.e., rate schedules) can be naturally concatenated into a longer rate schedule: The concatenation \(\mathbf {f}_1 \circ \cdots \circ \mathbf {f}_n\) is the rate schedule that follows \(\mathbf {f}_i\) for time \(i-1 \le t \lt i\) and zero elsewhere. Our game captures the idea that the Demon wins by making the constructed rate schedule not converge to the desired output, where the Demon is player I and the Chemist is player II. We allow the Demon to play any valid rate schedule, but the Chemist may have varying power as we discuss below.

The most unrestricted moveset contains all possible valid rate schedules. We call this case the Strong Demon Strong Chemist game.

We can also consider games where the players are restricted to only playing a subset of the valid rate schedules. We are particularly interested in games where the Demon is allowed to play any valid rate schedule as above, but the Chemist is restricted to only playing rate schedules from a particular restricted moveset M. It turns out that as long as the chemist has some fair rate schedule to play, the computational power is the same as with a Strong Chemist.

We first define the validity and fairness restrictions on movesets. Then given a restricted moveset for the Chemist, we define the corresponding game similar to definition 6.3.

It may seem that depending on what moves the Chemist is allowed to play, different classes of functions are computable (i.e., the Chemist has a winning strategy). However, the following theorem shows that the class of functions is invariant to the class of movesets that the Chemist has available, as long as there is always at least one fair move:

6.2 Non-zero Initial Context

Throughout this article, we have assumed that the only species allowed to be present at the start of the computation are the input species. Instead, one could consider a model where certain non-input species \(Z_1 \ldots Z_n\), called the initial context [18], have a fixed, nonzero rational concentration at the start of the computation. In this setting, we can clearly compute more functions than in the setting without initial context: For instance, we can easily compute \(f(x_1, \ldots , x_k) = C\) for some nonzero constant C, which is impossible without initial context because f is affine but not linear.

In fact, we can dual-rail compute any continuous piecewise rational affine function \(f: \mathbb {R}^k \rightarrow \mathbb {R}\), i.e., any function that is a rational linear function plus a rational constant: \(f(\mathbf {x}) = \mathbf {a} \cdot \mathbf {x} + b\). To see this, first note that we can compute any rational affine function by using the initial context to offset the value at \(f(\mathbf {0}) = b\). In fact, we can simply let output species \(Y^+\) and \(Y^-\) be the initial context, with \(\mathbf {x}(Y^+) = b\) initially if \(b\gt 0\) and \(\mathbf {x}(Y^-) = -b\) otherwise. Similar machinery to the proof of Lemma 5.15 can be used to extend this to continuous piecewise rational affine functions: By Theorem 5.11, any continuous piecewise rational affine function can be represented in max-min form, and then our construction from Section 5.3 shows that we can compute our given function f. In the direct computation setting, by a construction like the one in Section 5.4, we can compute any positive-continuous piecewise rational affine function \(f: \mathbb {R}^k_{\ge 0} \rightarrow \mathbb {R}_{\ge 0}\).

It also turns out that, even with initial context, we cannot compute any more functions than these. To see this, note that without loss of generality, we can assume that there is only one initial context species Z with initial concentration 1, since we can modify any CRC with initial context to include reactions that convert Z into \(Z_1,\ldots ,Z_n\) with appropriate concentrations.²³ Now let \(g(x_1, \ldots , x_k, z)\) be the value that the CRC computes when the input species have initial values \(x_1,\ldots , x_k\) and the initial context species has value z. A priori, g is only well-defined when \(z = 1\), but because paths remain valid after scaling, we know that \(\begin{equation*} g(x_1, \ldots , x_k, z) = z \cdot g(x_1/z, \ldots , x_k/z, 1) = z \cdot f(x_1/z, \ldots , x_k/z) \end{equation*}\) for any value of \(z \gt 0\). This shows that g is well-defined on \(D_U\) (recall Definition 5.8) for every U that contains Z, so we can apply the results of Sections 5.5 and 5.6 to characterize g on these domains. Using the fact that \(f(x_1, \ldots , x_k) = g(x_1, \ldots , x_k, 1)\) gives us the desired result.

Since every continuous function on a compact domain is uniformly continuous, it can be uniformly approximated by continuous piecewise rational affine functions. This shows that we can use rate-independent CRNs to approximate continuous functions. Note that for the negative argument above to work, it was important that all of the initial concentrations of the initial context species were rational. For practical purposes, this assumption is not at all restrictive, but it might be of theoretical interest to know what other functions can be computed if the initial concentrations are allowed to be arbitrary real numbers.