Understanding the algorithmic behaviors that are in principle realizable in a chemical system is necessary for a rigorous understanding of the design principles of biological regulatory networks. Further, advances in synthetic biology herald the time when we will be able to rationally engineer complex chemical systems and when idealized formal models will become blueprints for engineering.
Coupled chemical interactions in a well-mixed solution are commonly formalized as chemical reaction networks (CRNs). However, despite the widespread use of CRNs in the natural sciences, the range of computational behaviors exhibited by CRNs is not well understood. Here, we study the following problem: What functions f : ℝk → ℝ can be computed by a CRN, in which the CRN eventually produces the correct amount of the “output” molecule, no matter the rate at which reactions proceed? This captures a previously unexplored but very natural class of computations: For example, the reaction X1 + X2 → Y can be thought to compute the function y = min (x1, x2). Such a CRN is robust in the sense that it is correct whether its evolution is governed by the standard model of mass-action kinetics, alternatives such as Hill-function or Michaelis-Menten kinetics, or other arbitrary models of chemistry that respect the (fundamentally digital) stoichiometric constraints (what are the reactants and products?).
We develop a reachability relation based on a broad notion of “what could happen” if reaction rates can vary arbitrarily over time. Using reachability, we define stable computation analogously to probability 1 computation in distributed computing and connect it with a seemingly stronger notion of rate-independent computation based on convergence in the limit t → ∞ under a wide class of generalized rate laws. Besides the direct mapping of a concentration to a nonnegative analog value, we also consider the “dual-rail representation” that can represent negative values as the difference of two concentrations and allows the composition of CRN modules. We prove that a function is rate-independently computable if and only if it is piecewise linear (with rational coefficients) and continuous (dual-rail representation), or non-negative with discontinuities occurring only when some inputs switch from zero to positive (direct representation). The many contexts where continuous piecewise linear functions are powerful targets for implementation, combined with the systematic construction we develop for computing these functions, demonstrate the potential of rate-independent chemical computation.
1 Introduction
Understanding the dynamic behaviors that are, in principle, achievable with chemical species interacting over time is crucial for engineering of complex molecular systems capable of diverse and robust behaviors. The exploration of this space also helps to elucidate the constraints imposed upon biology by the laws of chemistry. The natural language for describing the interactions of molecular species in a well-mixed solution is that of chemical reaction networks (CRNs), i.e., finite sets of chemical reactions such as \(A + B \rightarrow A + C\). The intuitive meaning of this expression is that a unit of chemical species A reacts with a unit of chemical species B, producing a unit of a new chemical species C and regenerating a unit of A back. Typically (in mass-action kinetics) the rate with which this occurs is proportional to the product of the amounts of the reactants A and B.
Informally speaking, we can identify two sources of computational power in CRNs. First, the reaction stoichiometry transforms some specific ratios of reactants to products. For example, \(X \mathop {\rightarrow }\limits 2Y\) makes two units of Y for every unit of X. Second, in mass-action kinetics the reaction rate laws effectively perform multiplication of the reactant concentrations. In this work, we seek to disentangle the contributions of these two computational ingredients by focusing on the computational power of stoichiometry alone. Besides fundamental scientific interest, such rate-independent computation may be significantly easier to engineer than computation relying on rates (see Section 1.2). Importantly, stoichiometry is robust—not requiring the tuning of reaction conditions, nor even the assumption that the solution is well-mixed.
In the discrete model of chemical kinetics (see Section 1.2 for the distinction between the discrete and continuous models), rate-independence is formally related to probability 1 computation with passively mobile (i.e., interacting randomly) agents in distributed computing (the “population protocols” model [6, 9]; see Section 1.3). However, the continuous model of chemistry is most widely used and is more applicable for engineering chemical computation where working with bulk concentrations remains the state-of-the-art (see Section 1.2). This article formally articulates rate-independence in continuous CRNs and characterizes the computational power of stoichiometry in this model.
In the continuous setting, the amount of a species is a nonnegative real number representing its concentration (average count per unit volume).1 We characterize the class of real-valued functions computable by CRNs when reaction rates are permitted to vary arbitrarily (possibly adversarially) over time. Any computation in this setting must rely on stoichiometry alone. How can rate laws “preserve stoichiometry” while varying “arbitrarily over time”? Formally, preserving stoichiometry means that if we reach state \(\mathbf {d}\) from state \(\mathbf {c}\), then \(\mathbf {d}= \mathbf {c}+ \mathbf {M}\mathbf {u}\) for some non-negative vector \(\mathbf {u}\) of reaction fluxes, where the CRN’s stoichiometry matrix \(\mathbf {M}\) maps those fluxes to the changes in species concentrations they cause. (For example, flux 0.5 of reaction \(C+X \rightarrow C+3Y\) changes the concentrations of \(C,X,Y\), respectively, by \(0,-0.5,+1.5.\)) Subject to this constraint, the widest class of trajectories that still satisfies the intuitive meaning of the reaction semantics can be described informally as follows: (1) concentrations cannot become negative; (2) all reactants must be present when a reaction occurs (e.g., if a reaction uses a catalyst,2 then the catalyst must be present); (3) the causal relationships between the production of species is respected (e.g., if producing A requires B and producing B requires A, then neither can ever be produced if both are absent).3 This notion of “allowed trajectories” is formalized as Definition 2.22.
Fig. 1.
Fig. 1. Examples of rate-independent computation with chemical reaction networks. (a) Direct and (b) dual-rail CRNs computing the function \(f(x_1, x_2) = \max (x_1,x_2)\) plotted in (c). (d) Dual-rail CRN computing the function \(f(x_1, x_2) = x_1+\min (-x_1,x_2)-\max (x_1,x_2)\) plotted in (e). (f) Direct CRN computing the (discontinuous but still positive-continuous; see Definition 5.8) function \(f(x_1, x_2) = x_1\)if \(x_2 = 0\) and \(2 x_1\) if \(x_2 \gt 0\) plotted in (g).
The example shown in Figure 1(a) illustrates the style of computation studied here. Let \(f: \mathbb {R}_{\ge 0}^2 \rightarrow \mathbb {R}_{\ge 0}\) be the max function \(f(x_1, x_2) = \max (x_1,x_2)\) restricted to non-negative \(x_1\) and \(x_2\). The CRN of Figure 1(a) computes this function in the following sense: Inputs \(x_1\) and \(x_2\) are given as initial concentrations of input species \(X_1\) and \(X_2\). Then the CRN converges to f’s output value \(\max (x_1,x_2)\) of species Y, under a very wide interpretation of rate laws. Intuitively, the first two reactions must eventually produce \(x_1 + x_2\) of Y, and \(x_1\), \(x_2\) of \(Z_1\) and \(Z_2\), respectively. This is enforced by the stoichiometric constraint that the amount of \(Z_1\) and Y produced is equal to the amount of \(X_1\) consumed (and analogously for the second reaction). Stoichiometric constraints require the third reaction to produce the amount of K that is the minimum of the amount of \(Z_1\) and \(Z_2\) eventually produced in the first two reactions. Thus, \(\min (x_1,x_2)\) of K is eventually produced. Therefore, the fourth reaction eventually consumes \(\min (x_1,x_2)\) molecules of Y leaving \(x_1 + x_2 - \min (x_1, x_2) = \max (x_1, x_2)\) of Y behind. We can imagine an adversary pushing flux through these four reactions in any devious stratagem (i.e., arbitrary rates), yet unable to prevent the CRN from converging to the correct output, so long as applicable reactions must eventually occur.
We further consider the natural extension of such computation to handle negative real values. The example shown in Figure 1(b) computes \(f(x_1, x_2) = \max (x_1,x_2)\) (\(f: \mathbb {R}^2 \rightarrow \mathbb {R}\)), graphed in (c). To handle negative input and output values, we represent the value of each input and output by a pair of species using the so-called “dual-rail” representation. For example, in state \(\mathbf {c}\), \(x_1 = \mathbf {c}(X_1^+)-\mathbf {c}(X_1^-)\)—i.e., the difference between the concentrations of species \(X_1^+\) and \(X_1^-\). Note that when \(X_1^-\) and \(X_2^-\) are initially absent, the CRN becomes equivalent to the first three reactions of Figure 1(a) under relabeling of species. We do not need the last reaction of (a), because the output is represented as the difference of \(Y^+\) and \(Y^-\) by our convention. For the argument that the computation is correct even if \(X_1^-\) and \(X_2^-\) are initially present, we refer the reader to the proof of corollary 5.14 in section 5.3.
In addition to handling negative values, the dual-rail representation has the benefit of allowing composition. Specifically, the dual-rail representation allows a CRN to never consume its output species (e.g., rather than consuming \(Y^+\), it can produce \(Y^-\)). This monotonicity in the production of output allows directly composing CRN computations simply by concatenating CRNs and relabeling species (e.g., to make the output of one be input to the other). Since the upstream CRN never consumes its output species, the downstream CRN is free to consume them without interfering with the upstream computation. Since the class of functions computable by dual-rail CRNs ends up being invariant to whether or not, they are allowed to consume their output, and our results imply that dual-rail computation is composable without sacrificing computational power (see section 5.1).
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].4 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).5 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.6
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.7 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 Defining Reachability in Chemical Reaction Networks
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.8 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\).9 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).10 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
State \(\mathbf {d}\) is straight-line reachable (a.k.a. 1-segment reachable) from state \(\mathbf {c}\), written \(\mathbf {c}\rightarrow ^1\mathbf {d}\), if \((\exists \mathbf {u}\in \mathbb {R}_{\ge 0}^R)\ \mathbf {c}+ \mathbf {M}\mathbf {u}= \mathbf {d}\) and \(\mathbf {u}(\alpha) \gt 0\) only if reaction \(\alpha\) is applicable at \(\mathbf {c}\). In this case, write \(\mathbf {c}\rightarrow ^1_\mathbf {u}\mathbf {d}\).
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.
Definition 2.2.
Let \(k \in \mathbb {N}\cup \lbrace \infty \rbrace\). State \(\mathbf {d}\) is k-segment reachable from state \(\mathbf {c}\), written \(\mathbf {c}\rightsquigarrow ^k \mathbf {d}\), if \((\exists \mathbf {b}_0, \ldots , \mathbf {b}_{k})\ \mathbf {c}= \mathbf {b}_0 \rightarrow ^1\mathbf {b}_1 \rightarrow ^1\mathbf {b}_2 \rightarrow ^1\cdots \rightarrow ^1\mathbf {b}_{k}\), with \(\mathbf {b}_k = \mathbf {d}\) if \(k \in \mathbb {N}\), or \(\lim \nolimits _{i \rightarrow \infty } \mathbf {b}_i = \mathbf {d}\) if \(k = \infty\).
Definition 2.3.
State \(\mathbf {d}\) is segment-reachable (or simply reachable) from state \(\mathbf {c}\), written \(\mathbf {c}\rightsquigarrow \mathbf {d}\), if \((\exists k\in \mathbb {N}\cup \lbrace \infty \rbrace)\ \mathbf {c}\rightsquigarrow ^k \mathbf {d}\).
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.
Let \(\mathcal {C}= (\Lambda , R)\) be a CRN with species \(\Lambda\) and reactions R. For \(n \in \mathbb {N}\), we define a linear map \(\mathbf {x}_{n}: (\mathbb {R}^{\Lambda } \times \times _{i = 1}^{\infty } \mathbb {R}^{R}) \rightarrow \mathbb {R}^{\Lambda }\), which takes \({\boldsymbol \gamma }= (\mathbf {x}_0, \mathbf {u}_1, \mathbf {u}_2, \ldots)\) representing an initial state \(\mathbf {x}_0\) and reaction flux vectors \(\mathbf {u}_1, \mathbf {u}_2, \ldots\), and produces
which intuitively is the state reached after traversing the first n line segments. Let \(\Psi\) be the set of \({\boldsymbol \gamma }\) for which \(\lim \nolimits _{n \rightarrow \infty } \mathbf {x}_n({\boldsymbol \gamma })\) converges. We call elements of \(\Psi\)prepaths.
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.
Definition 2.5.
Let \(\Gamma _\infty\) be the subset of \(\Psi\) consisting of vectors \({\boldsymbol \gamma }= (\mathbf {x}_0, \mathbf {u}_1, \mathbf {u}_2, \ldots)\) satisfying the following conditions for all \(n\in \mathbb {N}\):
every reaction with positive flux in \(\mathbf {u}_{n+1}\) is applicable at \(\mathbf {x}_n\).
We call an element of \(\Gamma _\infty\) a piecewise linear path or sometimes just a path.
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.
Definition 2.6.
For \(k\in \mathbb {N}\), define \(\Gamma _k\) to be the subset of \(\Gamma _\infty\) consisting of all paths \({\boldsymbol \gamma }= (\mathbf {x}_0, \mathbf {u}_1, \mathbf {u}_2, \ldots)\) such that \(\mathbf {u}_i = \mathbf {0}\) for all \(i \gt k\). Say that a path is finite if it is contained in \(\Gamma _k\) for some \(k \in \mathbb {N}\).
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\).
Lemma 2.7.
For \(k \in \mathbb {N}\cup \lbrace \infty \rbrace\), \(\Gamma _k\) is convex.
Proof.
Let \({\boldsymbol \gamma }_0, {\boldsymbol \gamma }_1 \in \Gamma _\infty\) be two paths and consider \(\lambda \in (0,1)\). We need to show that \({\boldsymbol \gamma }_\lambda = (1 - \lambda){\boldsymbol \gamma }_0 + \lambda {\boldsymbol \gamma }_1\) is in \(\Gamma _\infty\). Recall that \(\mathbf {x}_n({\boldsymbol \gamma })\) is the state reached after the first n segments of path \({\boldsymbol \gamma }\). Note that for any \(n \in \mathbb {N}\),
Since \(\mathbb {R}_{\ge 0}^\Lambda\) is convex and both \(\mathbf {x}_n({\boldsymbol \gamma }_1)\) and \(\mathbf {x}_n({\boldsymbol \gamma }_2)\) are in \(\mathbb {R}_{\ge 0}^\Lambda\), we conclude that \(\mathbf {x}_n({\boldsymbol \gamma }_\lambda)\) is in \(\mathbb {R}_{\ge 0}^\Lambda\), too. Moreover, because both \(\lim _{n \rightarrow \infty } \mathbf {x}_n({\boldsymbol \gamma }_0)\) and \(\lim _{n \rightarrow \infty } \mathbf {x}_n({\boldsymbol \gamma }_1)\) converge, we see that
Below, for a path \({\boldsymbol \gamma }= (\mathbf {x}_0, \mathbf {u}_1, \mathbf {u}_2, \ldots)\), we use the notation \(\mathbf {u}_n({\boldsymbol \gamma })\) to represent the nth flux vector \(\mathbf {u}_n\) in \({\boldsymbol \gamma }.\) Since
any reaction \(\alpha\) occurs with positive flux in \(\mathbf {u}_n({\boldsymbol \gamma }_\lambda)\) only if \(\alpha\) occurs with positive flux in \(\mathbf {u}_n({\boldsymbol \gamma }_i)\) for \(i = 0\) or 1. Without loss of generality, suppose that \(\alpha\) occurs with positive flux in \(\mathbf {u}_n({\boldsymbol \gamma }_0)\). Then reaction \(\alpha\) is applicable at \(\mathbf {x}_{n-1}({\boldsymbol \gamma }_0)\), so the reactants are all present in positive concentrations in \(\mathbf {x}_{n - 1}({\boldsymbol \gamma }_0)\). This implies that they are present with positive concentrations in \(\mathbf {x}_{n - 1}({\boldsymbol \gamma }_\lambda)\) (note that we have excluded the case \(\lambda = 1\) from the outset). Therefore, reaction \(\alpha\) is applicable at \(\mathbf {x}_{n-1}({\boldsymbol \gamma }_\lambda)\). We conclude that every reaction occurring with positive flux in \(\mathbf {u}_n({\boldsymbol \gamma }_\lambda)\) is applicable at \(\mathbf {x}_{n-1}({\boldsymbol \gamma }_\lambda)\). This shows that \(\Gamma _\infty\) is convex.
To see that \(\Gamma _k\) for \(k \in \mathbb {N}\) is also convex, note that if \(\mathbf {u}_n({\boldsymbol \gamma }_0) = 0\) and \(\mathbf {u}_n({\boldsymbol \gamma }_1) = 0,\) then \(\mathbf {u}_n({\boldsymbol \gamma }_\lambda)\) will also be zero.□
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.
Lemma 2.8.
Let \(l \in \mathbb {N}\), \(k \in \mathbb {N}\cup \lbrace \infty \rbrace\), and let \(\mathbf {c},\mathbf {d}_1,\ldots ,\mathbf {d}_l\) be states such that \(\mathbf {c}\rightsquigarrow ^k \mathbf {d}_1\), \(\mathbf {c}\rightsquigarrow ^k \mathbf {d}_2\), \(\ldots\), and \(\mathbf {c}\rightsquigarrow ^k \mathbf {d}_l\). Then there exists \(\mathbf {d}\) such that \(\mathbf {c}\rightsquigarrow ^k \mathbf {d}\) and \([\mathbf {d}] = \bigcup _{i=1}^l [\mathbf {d}_i]\).
Proof.
Write \({\boldsymbol \gamma }_i\) for the path from \(\mathbf {c}\) to \(\mathbf {d}_i\); the convexity of \(\Gamma _k\) shows that the convex combination
\({\boldsymbol \gamma }\) exhibits a k-segment path from \(\mathbf {c}\) to \(\mathbf {d}\). If S is a species that is present at \(\mathbf {d}_i\) for any \(i,\) then S is also present at \(\mathbf {d}\). However, if S is present in none of the \(\mathbf {d}_i\), then S is not present in \(\mathbf {d}\). As a result, \([\mathbf {d}] = \bigcup _{i=1}^l [\mathbf {d}_i]\).□
Definition 2.9.
Given a state \(\mathbf {c}\), let \(\mathsf {P}(\mathbf {c})\) be the set of all species that are producible from \(\mathbf {c}\)—i.e., present in some state that is segment-reachable from \(\mathbf {c}\).
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}\).
Lemma 2.10.
Let m be the minimum of \(|R|\) and \(|\Lambda |\) and let \(\mathbf {c}\) be any state. Then there is a state \(\mathbf {d}\) such that \(\mathbf {c}\rightsquigarrow ^{m} \mathbf {d}\) and \([\mathbf {d}] = \mathsf {P}(\mathbf {c})\).
Proof.
Given a state \(\mathbf {c}\), let \(\mathsf {P}_{\lt \infty }(\mathbf {c})\) be the set of all species that are present in some state that is k-segment-reachable from \(\mathbf {c}\) for some \(k \lt \infty\).
We first argue that \(\mathsf {P}(\mathbf {c}) = \mathsf {P}_{\lt \infty }(\mathbf {c})\). Since clearly \(\mathsf {P}_{\lt \infty }(\mathbf {c}) \subseteq \mathsf {P}(\mathbf {c})\), it remains to show that \(\mathsf {P}(\mathbf {c}) \subseteq \mathsf {P}_{\lt \infty }(\mathbf {c})\). Let S be a species that is present in the state \(\mathbf {d}\) such that \(\mathbf {c}\rightsquigarrow ^\infty \mathbf {d}\). Then there exists a sequence \((\mathbf {b}_i)_i\) such that \(\mathbf {c}\rightarrow ^1 \mathbf {b}_1 \rightarrow ^1 \mathbf {b}_2 \rightarrow ^1 \ldots\) with \(\mathbf {d}= \lim _{i \rightarrow \infty } \mathbf {b}_i\). Because \(S \in [\mathbf {d}]\) there must be some \(i_0 \in \mathbb {N}\) where \(S \in [\mathbf {b}_{i_0}]\), and, since \(\mathbf {c}\rightsquigarrow ^{i_0} \mathbf {b}_{i_0}\), we see that \(S \in \mathsf {P}_{\lt \infty }(\mathbf {c})\). Thus, \(\mathsf {P}(\mathbf {c}) = \mathsf {P}_{\lt \infty }(\mathbf {c})\).
We show that the lemma holds for \(\mathsf {P}_{\lt \infty }(\mathbf {c})\); since \(\mathsf {P}(\mathbf {c}) = \mathsf {P}_{\lt \infty }(\mathbf {c}),\) this establishes the full lemma.
For all \(i \in \mathbb {N}\), let \(\Lambda _i\) be the set of species S such that there exists a \(\mathbf {d}\) with \(\mathbf {c}\rightsquigarrow ^i \mathbf {d}\) and \(S \in [\mathbf {d}]\). Similarly, let \(R_i\) be the set of reactions \(\alpha\) such that there exists a \(\mathbf {d}\) with \(\mathbf {c}\rightsquigarrow ^i \mathbf {d}\) and \(\alpha\) is applicable in \(\mathbf {d}\). Note that \(\Lambda _0 = [\mathbf {c}]\) and \(R_0\) is the set of reactions applicable in \(\mathbf {c}\). Also, since \(\mathbf {c}\rightsquigarrow ^i \mathbf {d}\) implies \(\mathbf {c}\rightsquigarrow ^{i + 1} \mathbf {d}\), we see that \(\Lambda _i \subseteq \Lambda _{i + 1}\) and \(R_i \subseteq R_{i + 1}\) for all i.
Now, we show that for all i there exists some \(\mathbf {x}_i\) such that \([\mathbf {x}_i] = \Lambda _i\) and \(\mathbf {c}\rightsquigarrow ^i \mathbf {x}_i\) (and therefore \(R_i\) consists of the reactions applicable at \(\mathbf {x}_i\)). To see this, for each \(S \in \Lambda _i\) let \(\mathbf {d}_S\) be a state such that \(\mathbf {c}\rightsquigarrow ^i \mathbf {d}_S\) and \(S \in [\mathbf {d}_S]\). By applying Lemma 2.8 to the set of all \(\mathbf {d}_S\), there is some \(\mathbf {d}\) such that \(\mathbf {c}\rightsquigarrow ^i \mathbf {d}\) and \([\mathbf {d}] = \Lambda _i\); this \(\mathbf {d}\) is our desired \(\mathbf {x}_i\).
Now, we will show that if \(\Lambda _i = \Lambda _{i + 1}\), then \(R_i = R_{i + 1}\) and, independently, if \(R_i = R_{i + 1}\), then \(\Lambda _{i + 1} = \Lambda _{i + 2}\) for all i. First suppose that \(\Lambda _i = \Lambda _{i + 1}\) and let \(\alpha\) be a reaction in \(R_{i + 1}\). Then there is some state \(\mathbf {d}\) such that \(\mathbf {c}\rightsquigarrow ^{i+1} \mathbf {d}\) and \(\alpha\) is applicable at \(\mathbf {d}\). Since all of the reactants of \(\alpha\) are present at \(\mathbf {d}\), they are a subset of \(\Lambda _{i + 1} = \Lambda _i\). They are therefore present at \(\mathbf {x}_i\), so \(\alpha\) is applicable at \(\mathbf {x}_i\). We conclude that \(\alpha \in R_i\) so \(R_{i + 1} = R_i\).
Now suppose that \(R_i = R_{i + 1}\) and let S be a species in \(\Lambda _{i + 2}\). Then there is some \(\mathbf {d}\) such that \(\mathbf {c}\rightsquigarrow ^{i + 2} \mathbf {d}\) and \(S \in [\mathbf {d}]\). If \(S \in [\mathbf {c}]\), then \(S \in \Lambda _{i + 1}\). Otherwise, S must be produced by some reaction \(\alpha\) in \(R_{i + 1} = R_i\), and we can apply \(\alpha\) to \(\mathbf {x}_i\) to obtain a state \(\mathbf {d}^{\prime }\) such that \(\mathbf {c}\rightsquigarrow ^i \mathbf {x}_i \rightarrow ^1 \mathbf {d}^{\prime }\) and \(S \in [\mathbf {d}^{\prime }]\). Again, we conclude that \(S \in \Lambda _{i + 1}\) so \(\Lambda _{i + 2} = \Lambda _{i + 1}\).
Combining the two statements we just proved, we see that if \(\Lambda _i = \Lambda _{i + 1}\), then \(\Lambda _i = \Lambda _j\) for all \(j \ge i\), so \(\Lambda _i = \mathsf {P}_{\lt \infty }(\mathbf {c})\). Similarly, if \(R_i = R_{i + 1}\), then \(\Lambda _{i + 1} = \mathsf {P}_{\lt \infty }(\mathbf {c})\).
If \(|\Lambda | \le |R|\), then, since \(\Lambda _0 \subseteq \Lambda _1 \subseteq \ldots\) is an increasing sequence of subsets of the finite set \(\Lambda\), it must be the case that \(\Lambda _j = \Lambda _{j + 1}\) for some \(j \le |\Lambda |\), and in this case \(\mathbf {x}_j\) gives our desired \(\mathbf {d}\). If, however, \(|R| \le |\Lambda |\), the proof is similar: First note that if \(R_0 = \emptyset\), then we are done. Otherwise, \(|R_0| \ge 1\) so, since \(R_i\) is an increasing sequence of subsets of \(R,\) there is some \(j \le |R| - 1\) such that \(R_j = R_{j + 1}\). Then \(\Lambda _{j + 1} = \mathsf {P}_{\lt \infty }(\mathbf {c})\) so \(\mathbf {x}_{j + 1}\) gives our desired \(\mathbf {d}\).□
Recall that a set is closed if it contains all of its limit points.
Lemma 2.11.
Let \(\mathbf {c}\in \mathbb {R}_{\ge 0}^\Lambda\) be any state and let \(S_\mathbf {c}\subseteq \mathbb {R}_{\ge 0}^\Lambda\) be the set of states that are straight-line reachable from \(\mathbf {c}\). Then \(S_\mathbf {c}\) is closed.
Proof.
Let \(R_\mathbf {c}\) be the set of reactions that are applicable at \(\mathbf {c}\). Then \(C = \lbrace \mathbf {u}\in \mathbb {R}_{\ge 0}^{R} | \mathbf {u}(\alpha) = 0\, \text{for}\, \alpha \notin R_\mathbf {c}\rbrace\) is a polyhedron. Then \(\mathbf {c}+ \mathbf {M}C\) is also a polyhedron (see Reference [53]), and is in particular closed. \(S_c\) is just \(\mathbb {R}_{\ge 0}^\Lambda \cap (\mathbf {c}+ \mathbf {M}C)\), and is therefore also closed.□
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.)
Definition 2.12.
State \(\mathbf {d}\) is s.s. segment reachable from state \(\mathbf {c}\), written \(\mathbf {c}\rightsquigarrow ^\infty _{\mathrm{ss}}\mathbf {d}\), if \((\exists \mathbf {b}^{\prime }_0,\mathbf {b}^{\prime }_1,\dots)\ \mathbf {c}= \mathbf {b}^{\prime }_0 \rightarrow ^1 \mathbf {b}^{\prime }_1 \rightarrow ^1 \dots\), where for some subsequence \(\mathbf {b}_0, \mathbf {b}_1, \dots\) of \(\mathbf {b}^{\prime }_0,\mathbf {b}^{\prime }_1,\dots\), \(\lim \nolimits _{i\rightarrow \infty } \mathbf {b}_i = \mathbf {d}\).
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.
Lemma 2.13.
If \([\mathbf {c}] = \mathsf {P}(\mathbf {c})\), then every state \(\mathbf {d}\) such that \(\mathbf {c}\rightsquigarrow ^\infty _{\mathrm{ss}}\mathbf {d}\) is straight-line reachable from \(\mathbf {c}\) (i.e., \(\mathbf {c}\rightsquigarrow ^\infty _{\mathrm{ss}}\mathbf {d}\) implies \(\mathbf {c}\rightarrow ^1 \mathbf {d}\)).
Proof.
First consider the finite case where \(\mathbf {c}\rightsquigarrow ^k \mathbf {d}\) for \(k \lt \infty\). Let \(\mathbf {u}_1, \ldots , \mathbf {u}_k\) be the flux vectors corresponding to the segments in the path from \(\mathbf {c}\) to \(\mathbf {d}\). Then \(\mathbf {u}= \mathbf {u}_1 + \cdots + \mathbf {u}_k\) is a vector in \(\mathbb {R}_{\ge 0}^R\). Every reaction that occurs with positive flux in \(\mathbf {u}\) has positive flux in one of the \(\mathbf {u}_i\), and thus its reactants are in \(\mathsf {P}(\mathbf {c})\), so are present in \(\mathbf {c}\) by the assumption \([\mathbf {c}] = \mathsf {P}(\mathbf {c})\). Thus, every reaction that occurs with positive flux in \(\mathbf {u}\) is applicable at \(\mathbf {c}\). The straight-line from \(\mathbf {c}\) corresponding to \(\mathbf {u}\) takes \(\mathbf {c}\) to
so \(\mathbf {c}\rightarrow _\mathbf {u}^1 \mathbf {d}\).
Now suppose that \(\mathbf {c}\rightsquigarrow ^\infty _{\mathrm{ss}}\mathbf {d}\). Then there is a sequence \(\mathbf {b}^{\prime }_1, \mathbf {b}^{\prime }_2, \dots\) of states such that \(\mathbf {c}\rightarrow ^1 \mathbf {b}^{\prime }_1 \rightarrow ^1 \mathbf {b}^{\prime }_2 \ldots\) and, for some subsequence \(\mathbf {b}_1, \mathbf {b}_2, \dots\) of \(\mathbf {b}^{\prime }_1, \mathbf {b}^{\prime }_2, \dots\), \(\mathbf {d}= \lim _{i \rightarrow \infty } \mathbf {b}_i\). For each finite \(i \ge 1\), for some finite \(j \ge i\), \(\mathbf {c}\rightsquigarrow ^j \mathbf {b}_i\). So, by the finite case shown above, \(\mathbf {c}\rightarrow ^1 \mathbf {b}_i\). Thus, \(\lbrace \mathbf {b}_1, \mathbf {b}_2, \ldots \rbrace \subseteq S_\mathbf {c}\), where \(S_\mathbf {c}\) is the set of states straight-line reachable from \(\mathbf {c}\). Since \(S_\mathbf {c}\) is closed by Lemma 2.11, it contains all its limit points, so \(\mathbf {d}\in S_\mathbf {c}\) as well, i.e., \(\mathbf {c}\rightarrow ^1 \mathbf {d}\).□
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.
Theorem 2.14.
If \(\mathbf {c}\rightsquigarrow ^\infty _{\mathrm{ss}}\mathbf {d}\), then \(\mathbf {c}\rightsquigarrow ^{m+1} \mathbf {d}\), where \(m=\min \lbrace |\Lambda |,|R|\rbrace\). Additionally, there is a constant K, depending only on the CRN, so the path from \(\mathbf {c}\) to \(\mathbf {d}\) can be chosen so the total flux of all reactions along the path is less than \(K\Vert \mathbf {d}- \mathbf {c}\Vert\).
Proof.
First, without loss of generality, we can consider the reduced CRN where we remove all of the reactions that are not used with positive flux in the given path from \(\mathbf {c}\) to \(\mathbf {d}\). By Lemma 2.10, we can find a state \(\mathbf {c}^{\prime }\) such that \(\mathbf {c}\rightsquigarrow ^m \mathbf {c}^{\prime }\) and \([\mathbf {c}^{\prime }] = \mathsf {P}(\mathbf {c})\). We now show that we can “scale down” the path \(\mathbf {c}\rightsquigarrow ^m \mathbf {c}^{\prime }\) such that no reaction occurs more than in the original path \(\mathbf {c}\rightsquigarrow \mathbf {d}\), allowing us to complete the path to \(\mathbf {d}\) using Lemma 2.13.
First consider the finite case, where \(\mathbf {c}\rightsquigarrow ^k \mathbf {d}\) for \(k \lt \infty\). We make the following general observation about finite paths: Let \({\boldsymbol \rho }\in \Gamma _l\) be any finite path with segments given by the flux vectors \(\mathbf {u}_1 \ldots \mathbf {u}_l\), and let \(F_{{\boldsymbol \rho }, \alpha }\) be the total flux through the reaction \(\alpha\) along \({\boldsymbol \rho }\), i.e.,
Let \({\boldsymbol \sigma }\in \Gamma _m\) be the path from \(\mathbf {c}\) to \(\mathbf {c}^{\prime }\), let \({\boldsymbol \sigma }_0 \in \Gamma _0 \subseteq \Gamma _m\) be the trivial path staying at \(\mathbf {c}\), and let \({\boldsymbol \gamma }\in \Gamma _k\) be the given path from \(\mathbf {c}\) to \(\mathbf {d}\). We can find some small \(\varepsilon \gt 0\) such that
for all reactions \(\alpha \in R\). As a result, letting \({\boldsymbol \gamma }^{\prime } = \varepsilon {\boldsymbol \sigma }+ (1 - \varepsilon) {\boldsymbol \sigma }_0\), we see that
Let \(\mathbf {a}\) be the state reached via \({\boldsymbol \gamma }^{\prime }\) (in particular \(\mathbf {c}\rightsquigarrow ^m \mathbf {a}\)). Since \(\mathbf {a}= \varepsilon \mathbf {c}^{\prime } + (1-\varepsilon) \mathbf {c}\) and \([\mathbf {c}^{\prime }] = \mathsf {P}(\mathbf {c})\) for the reduced CRN, we also have that \([\mathbf {a}] = \mathsf {P}(\mathbf {c})\). Thus, all reactions \(\alpha\) of the reduced CRN are applicable at \(\mathbf {a}\), and by Lemma 2.13, the final straight line from \(\mathbf {a}\) can be defined by the flux vector \(\mathbf {u}_\alpha = F_{{\boldsymbol \gamma }, \alpha } - F_{{\boldsymbol \gamma }^{\prime }, \alpha }\), so \(\mathbf {a}+ \mathbf {M}\mathbf {u}= \mathbf {d}\). This shows that \(\mathbf {a}\rightarrow ^1 \mathbf {d}\), so \(\mathbf {c}\rightsquigarrow ^{m+1} \mathbf {d}\), proving the theorem for the case of finitely many segments.
Now suppose that \(\mathbf {c}\rightsquigarrow ^\infty _{\mathrm{ss}}\mathbf {d}\), and let \({\boldsymbol \gamma }\in \Gamma _\infty\) be the path \(\mathbf {c}\rightarrow ^1_{\mathbf {u}_1} \mathbf {c}_1 \rightarrow ^1_{\mathbf {u}_2} \mathbf {c}_2 \ldots\) that starts at \(\mathbf {c}\) and has an infinite subsequence of states \(\mathbf {b}_1,\mathbf {b}_2,\dots\) converging to to \(\mathbf {d}\). Because we have assumed without loss of generality that every reaction \(\alpha\) occurs with positive flux along \({\boldsymbol \gamma }\), we know that
Let \(N = \max _{\alpha \in R} N_\alpha\) be the number of line segments required for each reaction to have had positive flux. The truncation of \({\boldsymbol \gamma }\) to a path with N segments from \(\mathbf {c}\) to \(\mathbf {c}_N\) is then a path \({\boldsymbol \gamma }_N\) where every reaction \(\alpha\) occurs with positive flux. By applying the first part of the argument, we can find a state \(\mathbf {a}\) with \([\mathbf {a}] = \mathsf {P}(\mathbf {c})\) so \(\mathbf {c}\rightsquigarrow ^m \mathbf {a}\) and \(\mathbf {a}\rightarrow ^1 \mathbf {c}_N\). But then because \(\mathbf {c}_N \rightsquigarrow ^\infty _{\mathrm{ss}}\mathbf {d}\), we see that \(\mathbf {a}\rightsquigarrow ^\infty _{\mathrm{ss}}\mathbf {d}\), so by Lemma 2.13, we see that \(\mathbf {a}\rightarrow ^1 \mathbf {d}\). As a result, we conclude that \(\mathbf {c}\rightsquigarrow ^{m + 1} \mathbf {d}\).
To see that we can bound the total flux along the path from \(\mathbf {c}\) to \(\mathbf {d}\), first note that by taking \(\varepsilon\) small enough, we can guarantee both that the total flux along \({\boldsymbol \gamma }^{\prime }\) is bounded by \(\Vert \mathbf {d}- \mathbf {c}\Vert\) and that
By the triangle inequality, this implies that \(\Vert \mathbf {a}- \mathbf {d}\Vert \lt 2\Vert \mathbf {d}- \mathbf {c}\Vert\). Now by applying Lemma D.4, we see that there’s some constant C depending only on the CRN so the flux vector \(\mathbf {u}\) of the straight line from \(\mathbf {a}\) to \(\mathbf {d}\) can chosen with \(\Vert \mathbf {u}\Vert \lt 2C \Vert \mathbf {d}- \mathbf {c}\Vert\). Taking \(K = 2C + 1\), we see that the flux along the whole path \(\mathbf {c}\rightsquigarrow ^m \mathbf {a}\rightarrow ^1 \mathbf {d}\) is bounded above by \(K\Vert \mathbf {d}- \mathbf {c}\Vert\).□
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.
If \(\mathbf {c}\rightsquigarrow \mathbf {d}\), then \(\mathbf {c}\rightsquigarrow ^{m+1} \mathbf {d}\), where \(m=\min \lbrace |\Lambda |,|R|\rbrace\).
Proof.
This follows from Theorem 2.14 and the fact that \(\mathbf {c}\rightsquigarrow \mathbf {d}\) implies \(\mathbf {c}\rightsquigarrow ^\infty _{\mathrm{ss}}\mathbf {d}\).□
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:
Corollary 2.16.
The relation \(\rightsquigarrow\) is transitive.
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.11
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*}
\begin{align} dy/dt &= k_2 c x, \end{align}
(2.3)
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,12 and for such CRNs, \(t_\mathrm{max}\) is finite.
Definition 2.17.
Fix an assignment of positive mass-action rate constants. Let \(\mathbf {c},\mathbf {d}\) be two states. We say \(\mathbf {d}\) is mass-action reachable (with respect to the rate constants) from \(\mathbf {c}\) if the associated mass action trajectory \({\boldsymbol \rho }\) obeys \({\boldsymbol \rho }(0) = \mathbf {c}\) and either \({\boldsymbol \rho }(t) = \mathbf {d}\) for some finite \(t \ge 0\) or \(\lim _{t\rightarrow \infty } {\boldsymbol \rho }(t) = \mathbf {d}\).13
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.
Definition 2.18.
Let \(\mathcal {C}=(\Lambda ,R)\) be a CRN. A siphon is a set of species \(\Omega \subseteq \Lambda\) such that, for all reactions \(\langle \mathbf {r},\mathbf {p}\rangle \in R\), \([\mathbf {p}] \cap \Omega \ne \emptyset \Rightarrow [\mathbf {r}] \cap \Omega \ne \emptyset\), i.e., every reaction that has a product in \(\Omega\) also has a reactant in \(\Omega\).
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).14 For the sake of completeness, we give a self-contained proof in Appendix A.
Lemma 2.19 ([5], Proposition 5.5).
Fix any assignment of positive mass-action rate constants. Let \(\Omega \subseteq \Lambda\) be a set of species. Then \(\Omega\) is a siphon if and only if, for any state \(\mathbf {c}\) such that \(\Omega \cap [\mathbf {c}] = \emptyset\) and any state \(\mathbf {d}\) that is mass-action reachable from \(\mathbf {c}\), \(\Omega \cap [\mathbf {d}] = \emptyset\).
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.
Lemma 2.20.
Let \(\Omega \subseteq \Lambda\) be a set of species. Then \(\Omega\) is a siphon if and only if, for any state \(\mathbf {c}\) such that \(\Omega \cap [\mathbf {c}] = \emptyset\) and any state \(\mathbf {d}\) such that \(\mathbf {c}\rightsquigarrow \mathbf {d}\), \(\Omega \cap [\mathbf {d}] = \emptyset\).
Proof.
To see the forward direction, suppose \(\Omega\) is a siphon, let \(\mathbf {c}\) be a state such that \([\mathbf {c}] \cap \Omega = \emptyset\), and let \(\mathbf {d}\) be such that \(\mathbf {c}\rightsquigarrow \mathbf {d}\). By Theorem 2.15, there is a finite path \({\boldsymbol \gamma }\) such that the nth line segment is between states \(\mathbf {x}_{n-1}({\boldsymbol \gamma })\) and \(\mathbf {x}_{n}({\boldsymbol \gamma })\), with \(\mathbf {c}= \mathbf {x}_{0}({\boldsymbol \gamma })\) and \(\mathbf {d}= \mathbf {x}_{m+1}({\boldsymbol \gamma })\). Assume inductively that \([\mathbf {x}_{n-1}({\boldsymbol \gamma })] \cap \Omega = \emptyset\); then no reaction applicable at \(\mathbf {x}_{n-1}({\boldsymbol \gamma })\) has reactants in \(\Omega\). So, by definition of siphon, no reaction applicable at \(\mathbf {x}_{n-1}({\boldsymbol \gamma })\) has products in \(\Omega\), and \([\mathbf {x}_{n}({\boldsymbol \gamma })] \cap \Omega = \emptyset\) as well. Therefore, \(\mathbf {d}\cap \Omega = \emptyset\). This shows the forward direction.
To show the reverse direction, suppose that \(\Omega\) is not a siphon. Then there is a reaction \(\alpha = \langle \mathbf {r},\mathbf {p}\rangle\) such that \([\mathbf {p}] \cap \Omega \ne \emptyset\), but \([\mathbf {r}] \cap \Omega = \emptyset\). Then from any state \(\mathbf {c}\) such that \([\mathbf {c}] = \Lambda \setminus \Omega\) (i.e., all species not in \(\Omega\) are present), all reactants of \(\alpha\) are present, so \(\alpha\) is applicable. Running \(\alpha\) produces S, hence results in a state \(\mathbf {d}\) such that \(\mathbf {c}\rightsquigarrow \mathbf {d}\) with \(\Omega \cap [\mathbf {d}] \ne \emptyset\), since \(S \in \Omega\).□
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.
Lemma 2.21.
If \(\mathbf {c}\) is any state, then \(\Omega = \Lambda \setminus \mathsf {P}(\mathbf {c})\) is a siphon.
Proof.
By Lemma 2.10 there is a state \(\mathbf {c}^{\prime }\) that is segment reachable from \(\mathbf {c}\) with all of the species in \(\mathsf {P}(\mathbf {c})\) present. If \(\Omega\) were not a siphon, then there would be a reaction \(\alpha\) that produced a species S of \(\Omega\) such that all of the reactants of \(\alpha\) would be contained in \(\Lambda \setminus \Omega = \mathsf {P}(\mathbf {c})\). This implies that \(\alpha\) would be applicable at \(\mathbf {c}^{\prime }\), so S would be in \(\mathsf {P}(\mathbf {c})\), giving a contradiction.□
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.
Definition 2.22.
A rate schedule is a function \(\mathbf {f}: \mathbb {R}_{\ge 0} \rightarrow \mathbb {R}_{\ge 0}^R\). We interpret \(\mathbf {f}_\alpha (t)\) to be the rate, or instantaneous flux, at which reaction \(\alpha\) occurs at time t. Given a state \(\mathbf {c}\in \mathbb {R}_{\ge 0}^\Lambda\), we say \(\mathbf {f}\) is a valid rate schedule starting at \(\mathbf {c}\) if:
(1)
(Total reaction fluxes are well-defined). For each \(\alpha \in R\), \(\mathbf {f}_\alpha\) is (locally Lebesgue) integrable: for each time \(t \ge 0\), \(\mathbf {F}_\alpha (t) = \int _0^t \mathbf {f}_\alpha (t) dt\) is well-defined and finite (although \(\int _0^\infty \mathbf {f}_\alpha (t) dt\) may be infinite). Let \(\mathbf {F}(t)\) be the vector in \(\mathbb {R}_{\ge 0}^R\) whose \(\alpha\) coordinate is \(\mathbf {F}_\alpha (t)\), which represents the total amount of each reaction flux that has happened by time t.15
Define the trajectory\({\boldsymbol \rho }: \mathbb {R}_{\ge 0}\rightarrow \mathbb {R}_{\ge 0}^\Lambda\) of \(\mathbf {f}\) starting at \(\mathbf {c}\) for all \(t \ge 0\) by \({\boldsymbol \rho }(t) = \mathbf {M}\cdot \mathbf {F}(t) + \mathbf {c}\), which represents the state of the CRN at time t.
(2)
(Positive-rate reactions require their reactants present). For all times \(t \ge 0\) and reactions \(\alpha \in R\), if \(\mathbf {f}_\alpha (t) \gt 0\), then \(\alpha\) is applicable in \({\boldsymbol \rho }(t)\).
(3)
(Absence of siphons is forward-invariant). For every siphon \(\Omega \subseteq \Lambda\), if \(\Omega \cap [{\boldsymbol \rho }(t)] = \emptyset\) for some time t, then \(\Omega \cap [{\boldsymbol \rho }(t^{\prime })] = \emptyset\) for all times \(t^{\prime } \ge t\).
Definition 2.23.
We say that a state \(\mathbf {d}\) is reachable from a state \(\mathbf {c}\) by a valid rate schedule if there is a valid rate schedule \(\mathbf {f}\) starting at \(\mathbf {c}\), with trajectory \({\boldsymbol \rho }\), such that either \(\mathbf {d}= {\boldsymbol \rho }(t_f)\) for some \(t_f \lt \infty\), or \(\mathbf {d}= \lim _{t \rightarrow \infty } {\boldsymbol \rho }(t)\). In the first case, we say that \(\mathbf {d}\) is reached in finite time.
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).16
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)\).
Lemma 2.24.
Suppose we fix an assignment of positive mass-action rate constants for a given CRN as well as an initial state \(\mathbf {c}\). Suppose that the associated mass action trajectory \({\boldsymbol \rho }\) is defined for all time. Then \(\mathbf {f}: \mathbb {R}_{\ge 0}\rightarrow \mathbb {R}_{\ge 0}^R\) such that \(\mathbf {f}(t) = \mathbf {A}({\boldsymbol \rho }(t))\) is a valid rate schedule whose trajectory is \({\boldsymbol \rho }\).
Proof.
First, note that because \({\boldsymbol \rho }(t)\) is a real analytic function, \(\mathbf {A}({\boldsymbol \rho }(t))\) is necessarily also real analytic, and therefore locally integrable, so condition (1) of Definition 2.22 is satisfied. Let \(\tilde{{\boldsymbol \rho }}(t)\) be the trajectory associated with the rate schedule \(\mathbf {f}\). Then, because \({\boldsymbol \rho }(t)\) is a solution to the mass-action ODEs with initial state \(\mathbf {c}\),
Because \(\mathbf {A}({\boldsymbol \rho }(t))\) can only be positive when \({\boldsymbol \rho }_X(t) \gt 0\) for all reactants X of \(\alpha\), we see that \(\mathbf {f}_\alpha (t) \gt 0\) implies that \([\tilde{{\boldsymbol \rho }}(t)] = [{\boldsymbol \rho }(t)]\) contains all of the reactants of \(\alpha\). Therefore, condition (2) of Definition 2.22 is satisfied. Finally, condition (3) of Definition 2.22 is satisfied by Lemma 2.19.□
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.
Observation 2.25.
If \(\mathbf {f}, \mathbf {g}\) are valid rate schedules, with \(\mathbf {f}\) finite such that \(\mathbf {f}(t) = \mathbf {0}\) for all \(t \ge t_0\), then \(\mathbf {h}\) defined by \(\mathbf {h}(t) = \mathbf {f}(t)\) for \(0 \le t \le t_0\) and \(\mathbf {h}(t) = \mathbf {g}(t-t_0)\) for \(t \gt t_0\), is a valid rate schedule.
The next lemma shows essentially that our definition of segment-reachability creates a valid rate schedule.
Lemma 2.26.
If \(\mathbf {c}\rightsquigarrow \mathbf {d}\), then \(\mathbf {d}\) is reachable from \(\mathbf {c}\) by a valid rate schedule.
Proof.
Since \(\mathbf {c}\rightsquigarrow \mathbf {d}\), by Theorem 2.15, we know that \(\mathbf {c}\rightsquigarrow ^{m + 1} \mathbf {d}\). Using induction and Observation 2.25, it suffices to verify that the rates defined by the straight-line reachability relation \(\mathbf {c}\rightarrow ^1_{\mathbf {u}} \mathbf {d}\) describe a valid rate schedule, since the rate schedules given by \(\rightsquigarrow ^{m+1}\) are simply concatenations of these. Define \(\mathbf {f}: \mathbb {R}_{\ge 0}\rightarrow \mathbb {R}_{\ge 0}^R\) by \(\mathbf {f}(t) = \mathbf {u}(\alpha)\) for \(0 \le t \lt 1\) and \(\mathbf {f}(t) = \mathbf {0}\) for all \(t \ge 1\). (In other words, for one unit of time, run the reactions at constant rates described by \(\mathbf {u}\).) Then \(\mathbf {f}\) is piecewise constant, and therefore integrable, so condition (1) of Definition 2.22 is satisfied. Next note that for all \(t \ge 1\), since \(\mathbf {f}_\alpha (t) = 0\) for all \(\alpha\), condition (2) of Definition 2.22 holds vacuously, and condition (3) holds because \({\boldsymbol \rho }(t) = {\boldsymbol \rho }(t^{\prime })\) for any \(t^{\prime } \gt t\). Now let \(t \lt 1\). Observe that
Since \(\mathbf {d}\) is a state (and thus non-negative on all species) and \(t\lt 1\), every species present with positive concentration in \(\mathbf {c}\) is present with positive concentration in \({\boldsymbol \rho }(t)\). Thus, all reactions applicable at \(\mathbf {c}\) are also applicable at \({\boldsymbol \rho }(t)\), so condition (2) of Definition 2.22 is also satisfied. Finally by Lemma 2.20, we see that condition (3) of Definition 2.22 is satisfied and therefore \(\mathbf {f}\) is a valid rate schedule starting at \(\mathbf {c}\). Since \({\boldsymbol \rho }(1) = \mathbf {d}\), we see that \(\mathbf {d}\) is reachable from \(\mathbf {c}\) by 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.
Theorem 2.27.
Given two states \(\mathbf {c}\) and \(\mathbf {d}\), \(\mathbf {d}\) is reachable from \(\mathbf {c}\) by a valid rate schedule if and only if \(\mathbf {c}\rightsquigarrow \mathbf {d}\).
Proof.
Lemma 2.26 establishes the reverse direction. To see the forward direction, let \(\mathbf {f}\) be the valid rate schedule from \(\mathbf {c}\) to \(\mathbf {d}\), and define \(\mathbf {F}\) and \({\boldsymbol \rho }\) for \(\mathbf {f}\) as in Definition 2.22.
First, suppose \(\mathbf {d}\) is reached in finite time \(t_f \in \mathbb {R}_{\ge 0}\), i.e., \(\mathbf {d}= {\boldsymbol \rho }(t_f)\).
We say that a reaction \(\alpha\)occurs with positive flux if \(\mathbf {F}_\alpha (t_f) \gt 0\). Let \(R_{\boldsymbol \rho }= \lbrace \alpha \in R\,|\,\mathbf {F}_\alpha (t_f) \gt 0\rbrace\) be the reactions that occur with positive flux along the trajectory \({\boldsymbol \rho }\), and let \(\Lambda _{\boldsymbol \rho }= \lbrace S \in \Lambda \mid {\boldsymbol \rho }_S(t) \gt 0 \text{ for some } 0 \le t \le t_f \rbrace\) be the species that are present with positive concentration at some point along the trajectory \({\boldsymbol \rho }\).
Consider removing species not in \(\Lambda _{\boldsymbol \rho }\) and reactions not in \(R_{\boldsymbol \rho }\). We claim that the pair \((\Lambda _{\boldsymbol \rho }, R_{\boldsymbol \rho })\) is a well-defined CRN, as defined in Section 2.1, because every reactant and product in \(R_\rho\) is in \(\Lambda _\rho\). To see why, let \(S \in \Lambda \setminus \Lambda _{\boldsymbol \rho }\) (i.e., \({\boldsymbol \rho }_S(t) = 0\) for all \(0 \le t \le t_f\)), let \(R_S\) be the reactions with S as a reactant, and let \(P_S\) be the reactions with S as a product; we must show \(R_S \cap R_{\boldsymbol \rho }= \emptyset\) and \(P_S \cap R_{\boldsymbol \rho }= \emptyset\). By Definition 2.22 part (2), no reaction in \(R_S\) has positive flux, so \(R_S \cap R_{\boldsymbol \rho }= \emptyset\). Since no reaction in \(R_S\) has positive flux, no reaction in \(P_S\) can have positive flux or else S would be produced with no reaction to consume it, contradicting \({\boldsymbol \rho }_S(t_f) = 0\), so \(P_S \cap R_{\boldsymbol \rho }= \emptyset\).
Now, we claim that every species in this reduced CRN \((\Lambda _{\boldsymbol \rho }, R_{\boldsymbol \rho })\) is segment-producible from \(\mathbf {c}\), i.e., \(\mathsf {P}(\mathbf {c}) = \Lambda _{\boldsymbol \rho }\). If not, then \(\Omega = \Lambda _{\boldsymbol \rho }\setminus \mathsf {P}(\mathbf {c})\) is non-empty. Letting \(S^*\) be some element of \(\Omega\), we know that \(S^*\) has positive concentration along \({\boldsymbol \rho }\) by our construction of the reduced CRN. However, by Lemma 2.21, \(\Omega\) is a siphon. Since \(\mathbf {c}\) is zero on \(\Omega\) and \(S^* \in \Omega\), this violates Definition 2.22 part (3).
Since every species in our reduced CRN is segment-producible from \(\mathbf {c}\), by Lemma 2.10, we can construct a state \(\mathbf {c}^{\prime }\) segment-reachable from \(\mathbf {c}\) where all of the species in the reduced CRN are present simultaneously. Since every reaction in \(R_{\boldsymbol \rho }\) is applicable at \(\mathbf {c}^{\prime }\), the remainder of the proof is similar to the proof for the finite case of Theorem 2.15: by “scaling down” the path from \(\mathbf {c}\) to \(\mathbf {c}^{\prime }\), there is a state \(\mathbf {a}\) such that \(\mathbf {c}\rightsquigarrow ^m \mathbf {a}\) (where \(m = \min \lbrace |\Lambda _{\boldsymbol \rho }|, |R_{\boldsymbol \rho }| \rbrace\)), \([\mathbf {a}] = \mathsf {P}(\mathbf {c}) = \Lambda _{\boldsymbol \rho }\), and \(\mathbf {a}\rightarrow ^1 \mathbf {d}\). Thus, \(\mathbf {c}\rightsquigarrow ^{m+1} \mathbf {d}\). This handles the case that \(\mathbf {d}\) is reached in finite time.
However, suppose that \(\mathbf {d}\) is not reached from \(\mathbf {c}\) in finite time, but instead \(\mathbf {d}= \lim _{t \rightarrow \infty } {\boldsymbol \rho }(t)\). This case is similar to the proof of the infinite case of Theorem 2.15. By definition of \(R_{\boldsymbol \rho }\), for each reaction \(\alpha \in R_{\boldsymbol \rho }\), \(\mathbf {F}_\alpha (t) \gt 0\) for some \(t \gt 0\). As a result, there is some finite \(t_\alpha\) such that \(\mathbf {F}_\alpha (t_\alpha) \gt 0.\) Let \(t_\mathrm{pos} = \max \lbrace t_\alpha \mid \alpha \in \mathbb {R}_{\boldsymbol \rho }\rbrace\), noting that \(\mathbf {F}(t_\mathrm{pos}) \gt 0\) for all \(\alpha \in R_{\boldsymbol \rho }\), i.e., each reaction has occurred by time \(t_\mathrm{pos}\). Let \(\mathbf {b}= {\boldsymbol \rho }(t_\mathrm{pos})\).
By applying the first part of the argument, we can find a state \(\mathbf {a}\) with \([\mathbf {a}] = \mathsf {P}(\mathbf {c})\) such that \(\mathbf {c}\rightsquigarrow ^m \mathbf {a}\) and \(\mathbf {a}\rightarrow ^1 \mathbf {b}\). Now let \(\mathbf {d}_t = {\boldsymbol \rho }(t)\) for every time \(t \gt t_\mathrm{pos}\). Because \({\boldsymbol \rho }\) restricted to \([t_\mathrm{pos}, t]\) gives a finite trajectory from \(\mathbf {b}\) to \(\mathbf {d}_t\), we know by the first part of the argument that \(\mathbf {b}\rightsquigarrow \mathbf {d}_t\), so \(\mathbf {a}\rightsquigarrow \mathbf {d}_t\). By Lemma 2.13, we see that \(\mathbf {a}\rightarrow ^1 \mathbf {d}_t\). Since \(\mathbf {d}= \lim _{t \rightarrow \infty } \mathbf {d}_t\), we see by Lemma 2.11 that \(\mathbf {a}\rightarrow ^1 \mathbf {d}\). As a result, we conclude that \(\mathbf {c}\rightsquigarrow \mathbf {d}\).□
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.17
Corollary 2.28.
Fix an assignment of positive mass-action rate constants for a given CRN. Let \(\mathbf {c},\mathbf {d}\) be two states such that \(\mathbf {d}\) is mass-action reachable from \(\mathbf {c}\). Then \(\mathbf {c}\rightsquigarrow \mathbf {d}\).
3 Stable Computation
We now use segment-reachability (Definition 2.3) to formalize what it means for a CRN to stably compute a function (Definition 3.2). The notion of stable computation is motivated by, and is essentially identical to, the definition of stable computation for population protocols and discrete CRNs [7, 18].
In this section, we justify stable computation by arguing for necessity: CRCs that we can reasonably call “rate-independent” must obey stable computation. Thus, stable computation is immediately useful for negative (impossibility) results: Showing a function cannot be stably computed implies it is not rate-independent in the desired intuitive sense. In the next section (Section 4), we address the other (sufficiency) direction and connect stable computation to another notion of computation based on convergence in the limit as \(t \rightarrow \infty\) that provides very strong guarantees for the desired rate-independent behavior of our constructions. Based on this connection, stable computation is taken as the primary definition of rate-independent computation in this work.
First, to formally define what it means for such a CRN to compute a function in any sense, we single out some aspects of the CRN as semantically meaningful. Formally, a chemical reaction computer (CRC) is a tuple \(\mathcal {C}= (\Lambda ,R,\Sigma ,\Gamma)\), where \((\Lambda ,R)\) is a CRN, \(\Sigma \subsetneq \Lambda\), written as \(\Sigma = \lbrace X_1,\ldots ,X_k\rbrace\),18 is the set of input species, and \(\Gamma \subsetneq \Lambda \setminus \Sigma\) is the set of output species. Input and output values can also be encoded indirectly via combinations of species. An important encoding for the purposes of this article will be the dual-rail representation (discussed in Section 5), which can handle both positive and negative quantities and allows for easier composition of CRN “modules.” Since we focus on single-output functions, we will have either a single output species \(\Gamma = \lbrace Y\rbrace\), or in the case of dual-rail computation, two output species \(\Gamma = \lbrace Y^+, Y^-\rbrace\).
We now define output-stable states and stable computation. Intuitively, output-stable states are “ideal” output states for rate-independent computation: The output is correct and no rate law can change it. Stable computation is then defined with respect to output-stable states by requiring that the correct output-stable state remains reachable no matter what “devious rate laws” may do. Although it is not obvious that the notion of output-stable states remains pertinent when transferred from the discrete setting to the continuous one (see the discussion at the beginning of Section 4), theorem 3.3 below and the subsequent results of Section 4 show that output stability remains crucial.
Definition 3.1.
A state \(\mathbf {o}\in \mathbb {R}_{\ge 0}^\Lambda\) is output-stable if, for all \(\mathbf {o}^{\prime }\) such that \(\mathbf {o}\rightsquigarrow \mathbf {o}^{\prime }\), \(\mathbf {o}\upharpoonright \Gamma = \mathbf {o}^{\prime } \upharpoonright \Gamma\), i.e., once \(\mathbf {o}\) is reached, no reactions can change the concentration of any output species.
Note that for a single output species Y, Definition 3.1 says that \(\mathbf {o}(Y) = \mathbf {o}^{\prime }(Y)\) for all \(\mathbf {o}^{\prime }\) such that \(\mathbf {o}\rightsquigarrow \mathbf {o}^{\prime }\). For the sake of brevity and readability, subsequently, we will state many definitions and formal theorem/lemma statements assuming there is only a single output species Y. In each case, there is a straightforward modification of the definition or result, so it applies to CRCs with multiple output species as well.
Definition 3.2.
Let \(f:\mathbb {R}_{\ge 0}^k \rightarrow \mathbb {R}_{\ge 0}\) be a function and let \(\mathcal {C}\) be a CRC. We say that \(\mathcal {C}\)stably computesf if, for all \(\mathbf {x}\in \mathbb {R}_{\ge 0}^k\), 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) = f(\mathbf {x})\).
To extend our results to functions with l outputs, we can compute l separate functions \(f_j:\mathbb {R}_{\ge 0}^k\rightarrow \mathbb {R}_{\ge 0}\) for \(1 \le j \le l\) with l independent CRCs and then combine them into a single CRC with l output species. In particular, we can use reactions like \(X_i \rightarrow X_i^1 + \cdots + X_i^l\) to copy input \(X_i\) to each of the l CRCs.
We now capture in a theorem the intuition that for a CRC to compute a function rate-independently in any reasonable sense, it must stably compute the function. The theorem says that if a CRC does not stably compute, then, no matter what you do, an adversary can “fight back” and make the output substantially (\(\epsilon\)) wrong. The proof uses the definitions of partial states and partial reachability, as well as Theorem E.1, which are in Appendix E.
Theorem 3.3.
Suppose a CRC does not stably compute \(f:\mathbb {R}_{\ge 0}^k \rightarrow \mathbb {R}_{\ge 0}\). Then there is \(\epsilon \gt 0\), input state \(\mathbf {x}\) and state \(\mathbf {z}\) reachable from \(\mathbf {x}\) such that for all \(\mathbf {o}\) reachable from \(\mathbf {z}\) there is \(\mathbf {o}^{\prime }\) reachable from \(\mathbf {o}\) such that \(| {\mathbf {o}^{\prime }(Y) - f(\mathbf {x})}| \gt \epsilon .\)
Proof.
We prove the contrapositive. Suppose that for all \(\epsilon \gt 0\), for any given input state \(\mathbf {x}\) and state \(\mathbf {z}\) such that \(\mathbf {x}\rightsquigarrow \mathbf {z}\), there exists a state \(\mathbf {o}\) such that \(\mathbf {z}\rightsquigarrow \mathbf {o}\) and for all \(\mathbf {o}^{\prime }\) such that \(\mathbf {o}\rightsquigarrow \mathbf {o}^{\prime }\), \(| {\mathbf {o}^{\prime }(Y) - f(\mathbf {x})} | \le \epsilon\).
For any input state \(\mathbf {x}\) and any \(\mathbf {z}\) reachable from \(\mathbf {x}\), first, we argue that there is an infinite sequence of states \(\mathbf {x}_{1/2}, \mathbf {x}_{1/3}, \mathbf {x}_{1/4}, \ldots\) such that \(\mathbf {x}\rightsquigarrow \mathbf {z}\rightsquigarrow \mathbf {x}_{1/2} \rightsquigarrow \mathbf {x}_{1/3} \rightsquigarrow \mathbf {x}_{1/4} \rightsquigarrow \dots\), and for all \(n \ge 2\), for all \(\mathbf {o}^{\prime }\) such that \(\mathbf {x}_{1/n} \rightsquigarrow \mathbf {o}^{\prime }\), \(| {\mathbf {o}^{\prime }(Y) - f(\mathbf {x})} | \le 1/n\). In other words, there is a sequence of states we can visit, where the adversary has less and less freedom to push the output away from the target value \(f(\mathbf {x})\). This is true by induction on n, choosing \(\epsilon = 1/(n+1)\), \(\mathbf {z}= \mathbf {x}_{1/n}\) and \(\mathbf {o}= \mathbf {x}_{1/(n+1)}\) in the above assumption.
By the definition of \(\mathbf {x}_{1/n}\), we see that \(|\mathbf {x}_{1/n}(Y) - f(\mathbf {x})| \le 1/n\), so \(\mathbf {x}_{1/n}(Y)\) converges to \(f(\mathbf {x})\) as \(n \rightarrow \infty\). Let \(\Delta = \left\lbrace Y\right\rbrace\) and let \(\mathbf {p}\in \mathbb {R}_{\ge 0}^\Delta\) be the partial state with \(\mathbf {p}(Y) = f(\mathbf {x})\). Then, we see that \(\mathbf {z}\rightsquigarrow ^\infty _{\mathrm{ss}}\mathbf {p}\) via \((\mathbf {x}_{1/i})_{i=2}^\infty\), so by theorem E.1, we can find a partition of \(\Lambda\) into \(\Lambda _\mathrm{b}\) and \(\Lambda _\mathrm{u}\) with \(Y \in \Lambda _\mathrm{b}\), a state \(\mathbf {y}\) (\(\mathbf {d}\) in Theorem E.1), and a subsequence \((\mathbf {x}^{\prime }_{i})_i\) of \(\left(\mathbf {x}_{1/n}\right)_n\) so \(\mathbf {z}\rightsquigarrow \mathbf {y}\) and the subsequence \((\mathbf {x}^{\prime }_{i})\) has the property that \(\mathbf {x}^{\prime }_i(S) \rightarrow \infty\) for all \(S \in \Lambda _\mathrm{u}\) and \(\mathbf {x}^{\prime }_i(S) \rightarrow \mathbf {y}(S)\) for all \(S \in \Lambda _\mathrm{b}\). Note that because \(Y \in \Lambda _\mathrm{b}\), we have \(\mathbf {y}(Y) = \lim _{i \rightarrow \infty } \mathbf {x}^{\prime }_i(Y) = f(\mathbf {x})\).
We now claim that \(\mathbf {y}\) is also output-stable, which is sufficient to prove the lemma as follows: Since \(\mathbf {y}\) is correct (\(\mathbf {y}(Y)=f(\mathbf {x})\)) and reachable from \(\mathbf {z}\), an arbitrary state reachable from input state \(\mathbf {x}\), this establishes that the CRC in fact stably computes f.
Suppose for the sake of contradiction that \(\mathbf {y}\) is not output-stable. Then there is some \(\mathbf {y}^{\prime }\) with \(\mathbf {y}\rightsquigarrow \mathbf {y}^{\prime }\) and \(|\mathbf {y}^{\prime }(Y) - \mathbf {y}(Y)| \gt \epsilon\). Let \(\mathbf {z}_i = \mathbf {x}^{\prime }_i - \mathbf {y}/2 + \mathbf {y}^{\prime }/2\). Because \(\mathbf {x}^{\prime }_i(S) \rightarrow \mathbf {y}(S)\) for all \(S \in \Lambda _\mathrm{b}\) and \(\mathbf {x}^{\prime }_i(S) \rightarrow \infty\) for all \(S \in \Lambda _\mathrm{u}\), there is some N so for all \(i \gt N\), we have \(\mathbf {x}^{\prime }_i \ge \mathbf {y}/2\). Then, by additivity of \(\rightsquigarrow\), for all \(i \gt N\),
where \(|\mathbf {z}_i(Y) - \mathbf {x}^{\prime }_i(Y)| = |\mathbf {y}^{\prime }(Y) - \mathbf {y}(Y)|/2 \gt \epsilon /2\). Since \((\mathbf {x}^{\prime }_{i})\) is a subsequence of \(\left(\mathbf {x}_{1/n}\right)\), if we take i large enough so \(\mathbf {x}^{\prime }_i\) is \(\mathbf {x}_{1/n_i}\) with \(1/n_i \lt \epsilon /4\), then we see that
but \(\mathbf {z}_i\) is reachable from \(\mathbf {x}_{1/n_i}\), which by definition can only reach states \(\mathbf {z}\) with \(|\mathbf {z}(Y) - f(\mathbf {x})| \le 1/n_i\), giving a contradiction.□
In some places, we will talk about CRCs with multiple output species representing multi-valued functions:
Corollary 3.4.
Suppose a CRC does not stably compute \(f: \mathbb {R}_{\ge 0}^k \rightarrow \mathbb {R}_{\ge 0}^l\). Then there is an \(\epsilon \gt 0\), an input state \(\mathbf {x}\) and a state \(\mathbf {z}\) reachable from \(\mathbf {x}\) so for all \(\mathbf {o}\) reachable from \(\mathbf {z}\) there is \(\mathbf {o}^{\prime }\) reachable from \(\mathbf {o}\) such that \(|\mathbf {o}^{\prime }(Y_i) - f(\mathbf {x})_i| \gt \epsilon\) for some \(1 \le i \le l\).
Proof.
This is similar to the proof of Theorem 3.3, but make \(\Delta = \left\lbrace Y_1, \ldots , Y_l\right\rbrace\) the set of all output species.□
4 Fair Computation
In the discrete model of CRN kinetics, if the set of states reachable from any input state is finite (i.e., the molecular counts are bounded as a function of the input state), then stable computation as in Definition 3.2 (the correct output-stable state is always reachable) is equivalent to the condition that the CRC is correct under standard stochastic kinetics with probability 1 (the correct output state state is actually reached) [7]. In the continuous CRN model, however, it might seem that the idea of stable computation is not strong enough to achieve intuitively “rate-independent” computation. There are at least two reasons for the concern.
First, it is possible that the output-stable state is always reachable but the mass-action trajectory does not converge to it. For example, consider the following CRC stably computing the identity function \(f(x) = x\):
\begin{align*} X + X &\rightarrow Y + Y \\ Y + X &\rightarrow X + X, \end{align*}
with X as the input species and Y as the output species. From any reachable state, we can reach the output-stable state with all X converted to Y. However, the mass-action trajectory converges to a dynamic equilibrium with \(k_2/(2 k_1 + k_2)\) fraction of X, where \(k_1\), \(k_2\) are rate constants of the two reactions.19 This shows that in general, showing that a CRC stably computes is not sufficient to claim that it computes rate-independently in any intuitive sense.
The second difficulty lies with the notion of output-stable states. While the notion of output-stable states is natural for discrete CRNs where we want the system to actually reach that state to end the computation, convergence to the output-stable state in continuous CRNs will typically be only in the limit \(t \rightarrow \infty\). For example, consider the following CRC:
\begin{align*} X &\rightarrow M + Z \\ M &\rightarrow Y \\ Z + Y &\rightarrow Z + M \\ Z &\rightarrow \emptyset . \end{align*}
The CRC stably computes \(f(x) = x\), since from any reachable state, we can reach the output-stable state \(\mathbf {o}\) with \(f(x)\) amount of Y by converting any remaining X to M, converting any remaining M to Y and completely draining Z. Note that \(\mathbf {o}\) is output-stable, since without Z, Y cannot be converted back to M. Further, under mass-action kinetics (for any choice of rate constants), the CRC converges to \(\mathbf {o}\), since as Z drains, the rate of the third reaction converges to 0. However, at every finite time \(0 \lt t\lt \infty\) in the mass-action trajectory, since Z is present, the state with zero amount of Y is reachable, so an adversary could substantially perturb the output. Thus, one would not call this CRC rate-independent to adversarial perturbations.
While in the previous section, we argued that stable computation is necessary for an intuitive notion of rate-independent computation, the examples above seem to suggest that it is not sufficient and that basing rate-independent computation entirely on stable computation could be ill-founded.
In this section, we develop an alternative approach to defining a very strong notion of “rate-independent” computation not based on stable computation, an approach we term “fair computation.” The approach is based on delineating a very broad class of rate laws, possibly adversarial, that still lead to convergence to the correct output. Based on the previous section (theorem 3.3), it is not surprising that CRCs that fail to stably compute also fail to fairly compute. What is more surprising, however, is that there is a strong connection in the other direction for a class of CRCs (feedforward)—for these CRCs stable computation implies fair computation. All our constructions will be in this class; thus, we obtain very strong rate-independence guarantees in the positive results part of this work. Combined with the results of the previous section, stable computation can thus be used as an easy-to-analyze proxy for proving both positive and negative results on rate-independent computation.
Intuitively, a CRC is said to fairly compute if it converges to the correct output for a broad class of rate laws, with the class being broad enough to capture adversarial behavior. To define the broad class of rate laws for fair computation, we start with the previously defined notion of valid rate schedules, which captures a very general class of chemical kinetics. Nonetheless, we must add an additional condition. In our original definition (Definition 2.22), the reaction rate \(\mathbf {f}_\alpha\) can vary arbitrarily over time as long as \(\alpha\) is applicable whenever \(\mathbf {f}_\alpha\) is positive. There is no requirement the other way—that a reaction must occur with positive rate if it is applicable—allowing for a greater variety of paths (e.g., segment paths with zero flux through some reactions). But, since there is nothing to prevent an adversary from “starving” reactions when they are applicable, preventing convergence, we now need to impose an additional requirement that we call fairness. We formalize this as a strictly positive lower bound \(\mathbf {H}\) on the reaction rate at states where the reaction is applicable. In particular, while the reaction rate vector \(\mathbf {f}(t)\) is a function of the time t, \(\mathbf {H}(\mathbf {c})\) is a function of the state\(\mathbf {c}\). We allow this lower bound to be violated occasionally, so long as it holds for an infinite measure of time. (For example, a fair rate schedule could starve applicable reactions on the unit time intervals \([0,1], [2,3], [4,5], \ldots\))
Definition 4.1.
Suppose \(\mathbf {f}: \mathbb {R}_{\ge 0}\rightarrow \mathbb {R}_{\ge 0}^R\) is a valid rate schedule starting at \(\mathbf {c}\). We say that \(\mathbf {f}\) is fair if there is a continuous function \(\mathbf {H}: \mathbb {R}_{\ge 0}^\Lambda \rightarrow \mathbb {R}_{\ge 0}^R\) such that, for all reactions \(\alpha \in R\), \(\mathbf {H}_\alpha (\mathbf {c}) \gt 0\) if and only if \(\alpha\) is applicable at \(\mathbf {c}\), and for some subset of times \(T_\alpha \subseteq \mathbb {R}_{\ge 0}\) of infinite measure, \(\mathbf {f}_\alpha (t) \ge \mathbf {H}_\alpha ({\boldsymbol \rho }(t))\) for all \(t \in T_\alpha\).
Requiring that the lower bound \(\mathbf {H}\) be continuous as a function of the state helps to ensure that if \(\mathbf {H}\) converges to zero then the point of convergence is a static state (no reaction is applicable): By continuity, the point of convergence must have \(\mathbf {H}= \mathbf {0}\), which implies that no reaction is applicable by definition 4.1. Note also that \(T_\alpha\) can depend on \(\alpha\); for example, Definition 4.1 allows for \(T_\alpha\) and \(T_\beta\) to be disjoint for different reactions \(\alpha ,\beta\) (i.e., whenever we run one applicable reaction, we starve another applicable reaction).
All typically considered rate laws such as mass-action, Michaelis-Menten, Hill function kinetics, and so on, are fair. We explicitly note this for mass-action CRNs with non-divergent trajectories:
Lemma 4.2.
The valid rate schedule for mass-action CRNs, as defined in Lemma 2.24, is fair if well-defined for all times.
Proof.
In the notation of lemma 2.24, take \(\mathbf {H}= \mathbf {A}\). Because we always assume that the rate constants of a mass action system are positive, \(\mathbf {A}_\alpha (\mathbf {c})\) is positive if and only if \(\alpha\) is applicable at \(\mathbf {c}\). Also, Lemma 2.24 shows that \(\mathbf {f}_\alpha (t) = \mathbf {A}_\alpha ({\boldsymbol \rho }(t))\) for all \(t \ge 0\), so certainly \(\mathbf {f}_\alpha (t) \ge \mathbf {A}_\alpha ({\boldsymbol \rho }(t))\).□
CRCs satisfying the definition below converge to the correct output despite actions of a very powerful adversary. Intuitively, the adversary is allowed to control the rates of all the reactions throughout the computation as long as applicable reactions are not entirely prevented from occurring.
Definition 4.3.
We say a CRC fairly computes a function \(f:\mathbb {R}_{\ge 0}^k \rightarrow \mathbb {R}_{\ge 0}\) if, for every input state \(\mathbf {x}\), every fair rate schedule starting at \(\mathbf {x}\) with trajectory \({\boldsymbol \rho }\) obeys \(\lim \nolimits _{t\rightarrow \infty } ({\boldsymbol \rho }(t) \upharpoonright \lbrace Y\rbrace) = f(\mathbf {x}).\)
There is a natural generalization of the above definition when \(f: \mathbb {R}_{\ge 0}^k \rightarrow \mathbb {R}_{\ge 0}^l\) has multiple outputs: Require the trajectory \({\boldsymbol \rho }\) to obey \(\lim _{t \rightarrow \infty } ({\boldsymbol \rho }(t) \upharpoonright \left\lbrace Y_1, \ldots , Y_l\right\rbrace) = f(\mathbf {x})\), where \(\left\lbrace Y_1, \ldots , Y_l\right\rbrace\) is the set of all output species.
We now want to establish the connection between fair computation and stable computation. The following lemma shows that in one direction, the connection is immediate: Fair computation is at least as strong as stable computation.
Lemma 4.4.
Any CRC that fairly computes a function f also stably computes f.
Proof.
We will prove the contrapositive: Suppose that a CRC \(\mathcal {C}\) does not stably compute f. We want to show that \(\mathcal {C}\) does not fairly compute f, and to do this, we will find an initial state \(\mathbf {x}\) and a fair rate schedule \(\mathbf {f}\) starting at \(\mathbf {x}\) so the trajectory \({\boldsymbol \rho }\) does not obey \(\lim _{t \rightarrow \infty }({\boldsymbol \rho }(t) \upharpoonright \left\lbrace Y\right\rbrace) = f(\mathbf {x})\).
Since we assumed that \(\mathcal {C}\) does not stably compute f, by Theorem 3.3, we know that there is some \(\epsilon \gt 0\) and some input state \(\mathbf {x}\) and state \(\mathbf {z}\) so \(\mathbf {x}\rightsquigarrow \mathbf {z}\) and every \(\mathbf {o}\) reachable from \(\mathbf {z}\) there is an \(\mathbf {o}^{\prime }\) reachable from \(\mathbf {o}\) with \(|\mathbf {o}^{\prime }(Y) - f(\mathbf {x})| \gt \epsilon\). In other words, it is possible to reach a state \(\mathbf {z}\), after which the output can be always made incorrect by some amount \(\epsilon\). By Lemma 2.26, we know that there is a valid rate schedule \(\mathbf {f}^{\prime }\) that reaches from \(\mathbf {x}\) to \(\mathbf {z}\) in finite time \(t_0\). Set \(\mathbf {f}= \mathbf {f}^{\prime }\) on the interval \([0, t_0)\).
Now choose some positive rate constant \(k_\alpha \gt 0\) for every reaction \(\alpha\), and on the interval \([t_0, t_0 + 1)\), set \(\mathbf {f}\) equal to the mass-action rate schedule corresponding to the rate constants \(k_\alpha\) and initial state \(\mathbf {z}\) (this is a valid rate schedule by Lemma 2.24). By Theorem 2.27, we know that the state \(\mathbf {o}\) that is reached at time \(t_0 + 1\) is segment reachable from \(\mathbf {z}\), so by the definition of \(\mathbf {z}\) there must be some \(\mathbf {o}^{\prime }\) with \(\mathbf {o}\rightsquigarrow \mathbf {o}^{\prime }\), where \(|\mathbf {o}^{\prime }(Y) - f(\mathbf {x})| \gt \epsilon\). Applying Lemma 2.26 again, we see that there is a valid rate schedule that reaches from \(\mathbf {o}\) to \(\mathbf {o}^{\prime }\) by time \(t_1\), so set \(\mathbf {f}\) equal to this rate schedule on \([t_0 + 1, t_1)\).
Alternating a mass-action rate schedule with an adversarial rate schedule in this way, we can find a sequence of times \(t_0, t_1, \ldots \in \mathbb {R}_{\ge 0}\) so \(\lim _{n \rightarrow \infty } t_n = \infty\) and a valid rate schedule \(\mathbf {f}\) so the associated trajectory \({\boldsymbol \rho }\) satisfies \(|{\boldsymbol \rho }(t_n)(Y) - f(\mathbf {x})| \gt \epsilon\) for every \(n \ge 1\). This proves that \(\lim _{t \rightarrow \infty } ({\boldsymbol \rho }(t) \upharpoonright \left\lbrace Y\right\rbrace) \ne f(\mathbf {x})\). Furthermore, by construction, we know that on every interval \([t_n, t_n + 1)\) the rate schedule \(\mathbf {f}\) is equal to the mass-action rate schedule with fixed rate constants \(k_\alpha\). Since \(\bigcup _{n \in \mathbb {N}} [t_n, t_n + 1)\) has infinite measure, Lemma 4.2 then shows that \(\mathbf {f}\) is fair. Since \(\mathbf {f}\) is a fair rate schedule failing to converge to the correct output, \(\mathcal {C}\) does not fairly compute f.□
It is natural to wonder if the converse of lemma 4.4 holds, i.e., whether every CRC that stably computes f also fairly computes f. This is not true in general, but the following section shows that it is true for the CRNs we will construct.
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 consumesS 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:
Definition 4.5.
A CRN \(\mathcal {C}= (\Lambda , R)\) is feedforward if \(\Lambda =\left\lbrace S_1, \ldots , S_n\right\rbrace\) and its stoichiometry matrix satisfies \(\mathbf {M}(i,j) \gt 0 \Rightarrow \exists (i^{\prime } \lt i) \mathbf {M}(i^{\prime },j) \lt 0\).
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.
Observation 4.6.
Every reaction in a feedforward CRN consumes some 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:
Lemma 4.7.
For a feedforward CRN, for each reaction \(\alpha\), there is a constant \(K_\alpha\) independent of \(\mathbf {u}\) such that \(\Vert \mathbf {M}\mathbf {u}\Vert \lt \epsilon\) implies that \(\mathbf {u}_\alpha \lt K_\alpha \epsilon\).
Proof.
We show by induction that we can find such a \(K_\beta\) for every reaction \(\beta\) that consumes a species. First suppose \(\beta\) is a reaction that consumes the first species \(S_1\) in the feedforward ordering. By the feedforward property, no reaction can produce \(S_1\), so for each reaction \(\alpha\), we have \(\mathbf {M}(S_1,\alpha) \le 0\). Thus,
Since the left-hand side of this equation is less than \(\epsilon\), and the right-hand side is a sum of non-negative terms, each term on the right-hand side must be less than \(\epsilon\). Since \(\beta\) consumes \(S_1\), \(|\mathbf {M}(S_1, \beta)| \ge 1\), so \(\mathbf {u}_\beta \lt \epsilon\) and we can take \(K_\beta = 1\). This establishes the base case for \(S_1\).
Now assume inductively that we have found an appropriate constant \(K_\alpha\) for every reaction \(\alpha\) that consumes a species \(S_i\) for \(i \lt n\), and suppose that \(\beta\) consumes \(S_n\). Then, because \(|(\mathbf {M}\mathbf {u})(S_n)| \lt \epsilon\), we know that
where for the second inequality we have used the feedforward condition to conclude that every reaction producing \(S_n\) must consume \(S_i\) for some \(i \lt n\), and therefore have flux bounded by \(K_\alpha \epsilon\) by inductive assumption. Like before, the leftmost term in the inequality is a sum of non-negative terms, so if we take
then \(\mathbf {u}_\beta \lt K_\beta \epsilon\). This shows that an appropriate \(K_\beta\) exists for every reaction \(\beta\) that consumes a species. But by Observation 4.6, every reaction consumes a species, so we are done.□
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.
Lemma 4.8.
Consider a feedforward CRN, and suppose \(\mathbf {f}\) is a valid rate schedule. Then the corresponding trajectory \({\boldsymbol \rho }\) converges to a state \(\mathbf {y}\) in the limit time \(t \rightarrow \infty\). Further, if \(\mathbf {f}\) is fair, then \(\mathbf {y}\) is static, i.e., no reaction is applicable in \(\mathbf {y}\).
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
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,
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.,
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
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
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
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,
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.
Lemma 4.9.
Any feedforward CRC that stably computes a function f also fairly computes f.
Proof.
By Lemma 4.8, under any fair rate schedule, the CRC converges from the initial state \(\mathbf {x}\) to a static state \(\mathbf {y}\). By Theorem 2.27, since \(\mathbf {y}\) is reachable from \(\mathbf {x}\) under a valid rate schedule, we know \(\mathbf {x}\rightsquigarrow \mathbf {y}\). Since \(\mathbf {y}\) is static, it is output-stable. Then, because the CRC stably computes f and \(\mathbf {x}\rightsquigarrow \mathbf {y}\), we must have \(\mathbf {y}(Y) = f(\mathbf {x})\) or else it would have stabilized to an incorrect output. So, the CRC fairly computes f.□
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.
Lemma 4.10.
If the CRC is feedforward, then the mass-action trajectory \({\boldsymbol \rho }:\mathbb {R}_{\ge 0}\rightarrow \mathbb {R}_{\ge 0}^\Lambda\) as defined in Lemma 2.24 is defined at all times.
Proof.
By the Escape Lemma [36, Lemma 9.19] it suffices to show that \({\boldsymbol \rho }\) remains in a compact subset of \(\mathbb {R}_{\ge 0}^\Lambda\) for all times when it is defined. Let \(V(\mathbf {x})\) be as in the proof of Lemma 4.8. We know by the argument of Lemma 4.8 that \(\frac{d}{dt} V({\boldsymbol \rho }(t)) \le 0\), so for all t, we know that \({\boldsymbol \rho }(t)\) lies in the subset of \(\mathbb {R}_{\ge 0}^\Lambda\) where \(V(\mathbf {x}) \le V(\rho (0))\). Since V is a linear function with strictly positive coefficients, this is a compact subset of \(\mathbb {R}_{\ge 0}^\Lambda\).□
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.
Corollary 4.11.
Consider a feedforward CRC stably computing a function f. Then for any input state \(\mathbf {x}\), for any state \(\mathbf {z}\) reachable from \(\mathbf {x}\), for any choice of reaction rate constants, the mass-action trajectory starting at \(\mathbf {z}\) is defined for all times and converges to an output-stable state \(\mathbf {y}\) in the limit \(t \rightarrow \infty\) such that \(\mathbf {y}(Y) = f(\mathbf {x})\).
5 The Computational Power of Stable Computation
This section presents the main results of our article, delineating the computational power of stable computation. As justified in Sections 3 and 4, we use stable computation as our primary notion of rate-independence. In section 5.1, we discuss the input and output representation of negative quantities and the composition of CRN modules via the “dual-rail representation.” In Section 5.2, we summarize our results on the computational power of stable computation for direct and dual-rail representations. These results are proven in subsequent sections, with positive (section 5.3 for dual-rail, Section 5.4 for direct) and negative directions (Section 5.4 for direct, section 5.6 for dual-rail) separately.
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\).
Definition 5.1.
We say that a CRC stably dual-computes\(f:\mathbb {R}^k \rightarrow \mathbb {R}\) if it stably computes a dual-rail representation \(\hat{f}:\mathbb {R}_{\ge 0}^k \times \mathbb {R}_{\ge 0}^k \rightarrow \mathbb {R}_{\ge 0}\times \mathbb {R}_{\ge 0}\) of f.
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:
Definition 5.2.
A CRC \(\mathcal {C}= (\Lambda ,R,\Sigma ,\Gamma)\) is output-oblivious if none of its output species \(\Gamma\) is a reactant in any reaction. In other words, for every \(\alpha =\langle \mathbf {r},\mathbf {p}\rangle \in R\) and \(Y \in \Gamma\), \(\mathbf {r}(Y) = 0\).
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.
We say a CRC fairly dual-computes a function \(f:\mathbb {R}^k \rightarrow \mathbb {R}\) if it fairly computes a dual-rail representation of f.
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:
Lemma 5.4.
If a CRC fairly dual-computes a function \(f: \mathbb {R}^k \rightarrow \mathbb {R}\), then it stably dual-computes f.
Proof.
The proof is similar to the proof of Lemma 4.4, except one uses Corollary 3.4 in the place of Theorem 3.3.□
Lemma 5.5.
If a feedforward CRC stably dual-computes a function \(f: \mathbb {R}^k \rightarrow \mathbb {R}\), then it fairly dual-computes f.
Proof.
The proof is the same as Lemma 4.9, since a static state is output-stable regardless of how many output species there are.□
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:
Definition 5.6.
A function \(f: \mathbb {R}^k \rightarrow \mathbb {R}\) is rational linear if there exist \(a_1,\ldots ,a_k \in \mathbb {Q}\) such that \(f(\mathbf {x}) = \sum _{i=1}^k a_i \mathbf {x}(i).\) A function \(f: \mathbb {R}^k \rightarrow \mathbb {R}\) is rational affine if there exist \(a_1,\ldots ,a_k, c \in \mathbb {Q}\) such that \(f(\mathbf {x}) = \sum _{i=1}^k a_i \mathbf {x}(i) + c\), i.e., f is a rational constant c plus a rational linear function.
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.7.
A function \(f: \mathbb {R}^k \rightarrow \mathbb {R}\) is piecewise rational linear (affine) if there is a finite set of partial rational linear (affine) functions \(f_1,\ldots ,f_p:\mathbb {R}^k \dashrightarrow \mathbb {R}\), with \(\bigcup _{j=1}^p {\rm dom} \;f_j = \mathbb {R}^k\), such that, for all \(j \in \lbrace 1,\ldots ,p\rbrace\) and all \(\mathbf {x}\in {\rm dom} \;f_j\), \(f(\mathbf {x}) = f_j(\mathbf {x})\). In this case, we say that \(f_1,\ldots ,f_p\) are the components of f.
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.
Theorem 5.9.
For a function \(f: \mathbb {R}_{\ge 0}^k \rightarrow \mathbb {R}_{\ge 0}\), the following are equivalent:
(1)
f is fairly computable by a CRC.
(2)
f is stably computable by a CRC.
(3)
f is positive-continuous piecewise rational linear.
(4)
f is stably computable by a feedforward CRC.
(5)
f is fairly computable by a feedforward CRC.
Proof.
(1) implies (2) is Lemma 4.4. (2) implies (3) is Lemma 5.38. (3) implies (4) is Lemma 5.16. (4) implies (5) is Lemma 4.9. (5) implies (1) is obvious.□
The following is the dual-rail analog of Theorem 5.9.
Theorem 5.10.
For a function \(f: \mathbb {R}^k \rightarrow \mathbb {R}\), the following are equivalent:
(1)
f is fairly dual-computable by a CRC.
(2)
f is stably dual-computable by a CRC.
(3)
f is continuous piecewise rational linear.
(4)
f is stably dual-computable by a feedforward, output-oblivious CRC.
(5)
f is fairly dual-computable by a feedforward, output-oblivious CRC.
Proof.
(1) implies (2) is Lemma 5.4. (2) implies (3) is Lemma 5.40. (3) implies (4) is Lemma 5.15. (4) implies (5) is Lemma 5.5. (5) implies (1) is obvious.□
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.
Theorem 5.11 ([42], Theorem 2.1).
Let \(D \subseteq \mathbb {R}^k\) be convex. For every continuous piecewise affine function \(f:D \rightarrow \mathbb {R}\) with components \(g_1,\ldots ,g_p\), there exists a family \(S_1,\ldots ,S_q \subseteq \lbrace 1,\ldots ,p\rbrace\) such that, for all \(\mathbf {x}\in D\), \(f(\mathbf {x}) = \max \nolimits _{i \in \lbrace 1,\ldots ,q\rbrace } \min \nolimits _{j \in S_i} g_j(\mathbf {x}).\)
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:
Lemma 5.12.
Rational linear functions are stably dual-computable by a feedforward, output-oblivious CRC.
Proof.
Let \(g:\mathbb {R}^k \rightarrow \mathbb {R}\) be a rational linear function \(g(\mathbf {x}) = \sum _{i=1}^k a_i \mathbf {x}(i)\). By clearing denominators, there exist \(n_1,\ldots ,n_k \in \mathbb {Z}\) and \(d \in \mathbb {Z}^+\) such that \(g(\mathbf {x}) = \frac{1}{d} \sum _{i=1}^k n_i \mathbf {x}(i).\) The following reactions compute a dual-rail representation of g with input species \(X_1^+,\ldots ,X_k^+, X_1^-,\ldots ,X_k^-\) and output species \(Y^+,Y^-\). For each i such that \(n_i \gt 0\), add the reactions
It is straightforward to verify that \(\hat{g}\) really is a dual-rail representation of g. To see that the above CRC stably computes \(\hat{g}\), define the functions \(p, q: \mathbb {R}_{\ge 0}^\Lambda \rightarrow \mathbb {R}\) so
It is also straightforward to verify that both p and q are preserved by all of the reactions in the above CRC. This shows that for all \(\mathbf {c}\) and \(\mathbf {d}\) with \(\mathbf {c}\rightsquigarrow \mathbf {d}\), we have \(p(\mathbf {d}) = p(\mathbf {c})\) and \(q(\mathbf {d}) = q(\mathbf {c})\). Now observe that from any state it is always possible to reach a state \(\mathbf {o}\) that only has positive concentrations of \(Y^+\) and \(Y^-\) by executing the reactions above to completion in the order in which they are listed. Such a state is evidently output-stable, so for any input state \(\mathbf {x}\) and any \(\mathbf {c}\) reachable from \(\mathbf {x}\), there is an output-stable state \(\mathbf {o}\) reachable from \(\mathbf {c}\). Furthermore, since \(\mathbf {o}\) only has positive concentrations of \(Y^+\) and \(Y^-\), we know that \(p(\mathbf {o}) = \mathbf {o}(Y^+)\) and \(q(\mathbf {o}) = \mathbf {o}(Y^-)\). Since \(\mathbf {o}\) is reachable from \(\mathbf {x}\), we see that
This shows that a dual-rail representation of any rational linear function can be stably dual-computed by a CRC. Moreover, the CRC above is clearly output-oblivious, and it is feedforward under the ordering \(X_1^+ \lt \cdots \lt X_n^+ \lt X_1^- \lt \cdots X_n^- \lt W^+ \lt W^- \lt Y^+ \lt Y^-\).□
Lemma 5.13.
Min is stably dual-computable by a feedforward, output-oblivious CRC.
Proof.
The following reactions stably compute a dual-rail representation of \(\min\) with input species \(X_1^+\), \(X_2^+\), \(X_1^-\), \(X_2^-\) and output species \(Y^+\), \(Y^-\).
It is straightforward to verify that \(\hat{f}\) is really a dual-rail representation of \(\min\). To see that the above CRC stably computes \(\hat{f}\),22 define the functions \(p, q, \delta : \mathbb {R}_{\ge 0}^\Lambda \rightarrow \mathbb {R}\) so
It is also straightforward to verify that the above three functions are preserved by all of the reactions in the given CRC. Note that by running 5.3, then 5.2, then 5.1 to completion it is always possible to reach a state \(\mathbf {o}\) with \(\mathbf {o}(X_1^-) = \mathbf {o}(X_2^-) = 0\) and also \(\mathbf {o}(X_1^+) = 0\) or \(\mathbf {o}(X_2^+) = 0\). Such a state is evidently output-stable, so for any input state \(\mathbf {x}\) and any \(\mathbf {c}\) reachable from \(\mathbf {x}\) there is an output-stable state \(\mathbf {o}\) reachable from \(\mathbf {c}\). Since \(\mathbf {o}\) is reachable from \(\mathbf {x}\), we know that
Now suppose without loss of generality that \(x_1^+ + x_2^- \le x_2^+ + x_1^-\) (the analysis of the other case is similar). Then \(\delta (\mathbf {x}) \le 0\), so \(\mathbf {o}(X_1^+) - \mathbf {o}(X_2^+) = \delta (\mathbf {o}) \le 0\). If \(\mathbf {o}(X_1^+)\) were positive, then by the definition of \(\mathbf {o}\), we would know that \(\mathbf {o}(X_2^+) = 0\), contradicting the fact that \(\delta (\mathbf {o}) \le 0\). Thus, \(\mathbf {o}(X_1^+) = 0\) and
Solving the above equation for \(\mathbf {o}(Y^+)\) shows that \(\mathbf {o}(Y^+) = x_1^+ + x_2^- = \min (x_1^+ + x_2^-, x_2^+ + x_1^-)\). This shows that the above CRC stably computes \(\hat{f}\). The CRN is clearly output-oblivious and it is feedforward with the ordering \(X_1^- \lt X_2^- \lt X_1^+ \lt X_2^+ \lt Y^- \lt Y^+\).□
Corollary 5.14.
Max is stably dual-computable by a feedforward, output-oblivious CRC.
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
□
Lemma 5.15.
Let \(D \subseteq \mathbb {R}^k\) be convex, and let \(f:D \rightarrow \mathbb {R}\) be a continuous piecewise rational linear function. Then f is stably dual-computed by a feedforward, output-oblivious CRC.
Proof.
By Theorem 5.11, we know that any such continuous piecewise rational linear function can be represented as a composition of \(\max\), \(\min\), and rational linear functions. Moreover, the \(\min\) and \(\max\) functions with two arguments can be composed in a tree of depth \(\log l\) to compute the minimum or maximum functions with input arity l. Since multiple rational linear functions may use the same inputs, we also need a fan-out module (likewise feedforward, output-oblivious) that copies a single input to multiple outputs:
where p is the number of linear components of the piecewise linear function f. Because output-oblivious CRCs are composable, and because the composition of feedforward, output-oblivious CRCs is again feedforward and output-oblivious, we can compose the CRCs from Lemma 5.12, Lemma 5.13, and Corollary 5.14 to produce a feedforward, output-oblivious CRC that computes f.□
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.)
Lemma 5.16.
Every positive-continuous piecewise rational linear function \(f:\mathbb {R}_{\ge 0}^k \rightarrow \mathbb {R}_{\ge 0}\) is stably computable by a feedforward CRC.
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.1 Siphons and Output Stability.
To characterize stable function computation for CRCs, we will crucially rely on the notion of siphons, which we recall from Section 2.4, Definition 2.18. Lemma 5.18 shows the underlying relationship between output stability and siphons.
Let \(\mathcal {C}=(\Lambda ,R,\Sigma ,\lbrace Y\rbrace)\) be a CRC. We call a siphon \(\Omega\)stabilizing if, for any state \(\mathbf {d}\), \(\mathbf {d}\upharpoonright \Omega = \mathbf {0}\) implies that \(\mathbf {d}\) is output-stable. In other words, “draining” \(\Omega\) (removing all of its species) causes the output to stabilize.
Lemma 5.17.
If \(\mathbf {c}\) is an output-stable state, then \(\Omega _\mathbf {c}= \Lambda \setminus \mathsf {P}(\mathbf {c})\) is a stabilizing siphon.
Proof.
By Lemma 2.21, \(\Omega _\mathbf {c}\) is a siphon for any state \(\mathbf {c}\). We show the contrapositive that if \(\Omega _\mathbf {c}\) is not stabilizing, then \(\mathbf {c}\) is not output-stable. Suppose there is some particular state \(\mathbf {d}\) with \(\mathbf {d}\upharpoonright \Omega _\mathbf {c}= \mathbf {0}\) that is not output-stable, i.e., for some \(\mathbf {e}\) such that \(\mathbf {d}\rightsquigarrow \mathbf {e}\), we have \(\mathbf {d}(Y) \ne \mathbf {e}(Y)\). Because of this, there must be some reaction \(\alpha\) applicable at a state \(\mathbf {d}^{\prime }\) reachable from \(\mathbf {d}\) such that \(\alpha\) changes the amount of Y (in other words, \(\mathbf {M}_{Y, \alpha } \ne 0\)). Since \(\Omega _\mathbf {c}\) is a siphon absent in \(\mathbf {d}\), by Lemma 2.20, \(\mathbf {d}^{\prime } \upharpoonright \Omega _{\mathbf {c}} = \mathbf {0}\), so all reactants of \(\alpha\) must be contained in \(\mathsf {P}(\mathbf {c})\). By Lemma 2.10, we can find a state \(\mathbf {c}^{\prime }\) such that \(\mathbf {c}\rightsquigarrow \mathbf {c}^{\prime }\) and \([\mathbf {c}^{\prime }] = \mathsf {P}(\mathbf {c})\). Since \(\alpha\) is applicable at \(\mathbf {c}^{\prime }\), we see that \(\mathbf {c}\) is not output-stable.□
The next lemma shows the key property of siphons that we will use to reason about stably computing CRC’s: that they characterize the output-stable states, i.e., the only way for the output to stabilize is to drain some stabilizing siphon.
Lemma 5.18.
There is a set of stabilizing siphons \(\mathcal {S}\) such that a state \(\mathbf {c}\) is output-stable if and only if \(\exists \Omega \in \mathcal {S}\) such that \(\mathbf {c}\upharpoonright \Omega = \mathbf {0}\).
Proof.
Take \(\mathcal {S}\) to be \(\left\lbrace \Omega _\mathbf {c}\subseteq \Lambda \ |\ \mathbf {c}\text{ is output-stable}\right\rbrace\), where \(\Omega _\mathbf {c}= \Lambda \setminus \mathsf {P}(\mathbf {c})\). By Lemma 5.17, each \(\Omega \in \mathcal {S}\) is a stabilizing siphon, so if \(\exists \Omega \in \mathcal {S}\) with \(\mathbf {d}\upharpoonright \Omega = \mathbf {0,}\) then \(\mathbf {d}\) is output-stable. However, if \(\mathbf {d}\) is some output-stable state, then by the definition of \(\Omega _\mathbf {d}\) (the set of species that cannot be produced from \(\mathbf {d}\)), we have \(\mathbf {d}\upharpoonright \Omega _{\mathbf {d}} = \mathbf {0}\). By our construction of \(\mathcal {S}\), we know that \(\Omega _{\mathbf {d}} \in \mathcal {S}\).□
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
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,
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
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:
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.5.4 Rationality.
Recall Definition 5.6 defining rational linear functions and Definition 5.7 defining piecewise rational linear functions.
The main ideas of this section are as follows: To show that a function is piecewise rational linear, we need to show that it is rational linear on some finite set of domains that cover the input space.
A linear function \(f: \mathbb {R}^n \rightarrow \mathbb {R}\) that sends \(\mathbb {Q}^n\) to \(\mathbb {Q}\) is necessarily rational linear. Since a linear function on \(\mathbb {R}^n\) is completely determined by its behavior on any open ball, we can check this condition “locally” on any domain that contains an open ball. (Since f is continuous and all of the points of domains that do not contain an open ball are limit points of the other domains, we can ignore domains that do not contain open balls.) The fact that the function sends \(\mathbb {Q}^n\) to \(\mathbb {Q}\) on such a domain is ultimately a consequence of the fact that the stoichiometry matrix of a CRN has only integer coefficients, so it preserves rationality.
Recall that \(\Psi\), defined in Definition 2.4, is the space of all prepaths, and \(\Gamma _\infty\), defined in Definition 2.5, is the space of all paths.
Definition 5.33.
A path \({\boldsymbol \gamma }\in \Gamma _\infty\) is a rational path if it has rational initial concentrations and all of its segments have rational fluxes. In other words, \(\mathbf {x}_0({\boldsymbol \gamma }) \in \mathbb {Q}^\Lambda\) and \(\mathbf {u}_i({\boldsymbol \gamma }) \in \mathbb {Q}^R\) for all i.
Note that, since the stoichiometry matrix is an integer-valued matrix, it is automatically the case that \(\mathbf {o}({\boldsymbol \gamma })\) and every \(\mathbf {x}_i({\boldsymbol \gamma })\) is in \(\mathbb {Q}^\Lambda\) for any rational path \({\boldsymbol \gamma }\).
Definition 5.34.
We say that two prepaths \({\boldsymbol \gamma }, {\boldsymbol \gamma }^{\prime } \in \Psi\)have the same sign if for all species S, reactions \(\alpha\), and \(i \in \mathbb {N}\), it is the case that \(\operatorname{sgn}\mathbf {x}_i({\boldsymbol \gamma })_S = \operatorname{sgn}\mathbf {x}_i({\boldsymbol \gamma }^{\prime })_S\) and \(\operatorname{sgn}\mathbf {u}_i({\boldsymbol \gamma })_\alpha = \operatorname{sgn}\mathbf {u}_i({\boldsymbol \gamma }^{\prime })_\alpha\).
Lemma 5.35.
Let \({\boldsymbol \gamma }\in \Gamma _\infty\) be a finite piecewise linear path. Then for any \(\varepsilon \gt 0\), there is a rational path \({\boldsymbol \gamma }^{\prime } \in \Gamma _\infty\) such that \({\boldsymbol \gamma }^{\prime }\) has the same sign as \({\boldsymbol \gamma }\) and \(||{\boldsymbol \gamma }^{\prime } - {\boldsymbol \gamma }|| \lt \varepsilon\). If \({\boldsymbol \gamma }\) already has rational initial concentrations, then \({\boldsymbol \gamma }^{\prime }\) can be chosen with the same initial concentrations.
Proof.
First, let N be the largest natural number such that \(\mathbf {u}_N({\boldsymbol \gamma }) \ne 0\). (Such an N exists, since \({\boldsymbol \gamma }\) is finite.) For any reaction \(\alpha \in R\) and any \(n \in \mathbb {N}_{\gt 0}\) such that \(\mathbf {u}_n({\boldsymbol \gamma })_\alpha = 0\), set \(\mathbf {u}_n({\boldsymbol \gamma })_\alpha = 0\). Now for each species \(S \in \Lambda\) such that \(\mathbf {x}_k({\boldsymbol \gamma })_S = 0\) and \(k \le N\), consider the following linear equation:
Aggregating these equations for all \(0 \le k \le N\) and all \(S \in \Lambda\) such that \(\mathbf {x}_k({\boldsymbol \gamma })_S = 0\) gives a system of equations, linear in \(\mathbf {x}_0({\boldsymbol \gamma }^{\prime })_S\) and \(\mathbf {u}_n({\boldsymbol \gamma }^{\prime })_\alpha\), with rational coefficients. This equation has a real-valued solution, namely, \(\mathbf {x}_0({\boldsymbol \gamma }^{\prime }) = \mathbf {x}_0({\boldsymbol \gamma })\) and \(\mathbf {u}_n({\boldsymbol \gamma }^{\prime }) = \mathbf {u}({\boldsymbol \gamma })\), so by Lemma C.1 (proven in Appendix C), it must have a solution with rational coefficients that is \(\delta\)-close for any \(\delta \gt 0\). By taking \(\delta\) small enough, we can of course make \(\delta \lt \varepsilon\), but we can also guarantee that \(\mathbf {x}_0({\boldsymbol \gamma }^{\prime })_S\) is positive whenever \(\mathbf {x}_0({\boldsymbol \gamma })_S\) is positive and similarly for \(\mathbf {u}_n({\boldsymbol \gamma }^{\prime })_\alpha\). We have therefore specified a \({\boldsymbol \gamma }^{\prime } \in \Psi\) with the same sign as \({\boldsymbol \gamma }\). Since \({\boldsymbol \gamma }^{\prime }\) has the same sign as \({\boldsymbol \gamma }\), and, since \({\boldsymbol \gamma }\) is a valid path, we conclude that \({\boldsymbol \gamma }^{\prime }\) is also a valid path, so \({\boldsymbol \gamma }^{\prime } \in \Gamma _\infty\).
If \(\mathbf {x}_0({\boldsymbol \gamma })\) is already in \(\mathbb {Q}^\Lambda\), then the same argument applies, with the modification that you fix \(\mathbf {x}_0({\boldsymbol \gamma }^{\prime }) = \mathbf {x}_0({\boldsymbol \gamma })\), and instead solve the inhomogeneous system of equations
when \(\mathbf {x}_k({\boldsymbol \gamma })_S = 0\).□
Lemma 5.36.
Let \(\Omega\) be a stabilizing siphon. If \(\Sigma (\Omega)\) contains an open ball, then f is rational linear when restricted to inputs in \(\Sigma (\Omega)\).
Proof.
Let B be the open ball contained in \(\Sigma (\Omega)\) and let \(\mathbf {x}\) be in \(\mathbb {Q}^\Lambda \cap B\). We know that there is a piecewise linear path \({\boldsymbol \gamma }\) starting at \(\mathbf {x}\) such that \(f(\mathbf {x}) = \mathbf {o}({\boldsymbol \gamma })_Y\). By Theorem 2.15, we may assume without loss of generality that \({\boldsymbol \gamma }\) is finite. By Lemma 5.35, there is a rational path \({\boldsymbol \gamma }^{\prime }\) with the same sign as \({\boldsymbol \gamma }\) such that \(\mathbf {x}_0({\boldsymbol \gamma }^{\prime }) = \mathbf {x}_0({\boldsymbol \gamma }) = \mathbf {x}\). Because \(\mathbf {o}({\boldsymbol \gamma })\) is an output-stable state, some siphon \(\Omega\) is drained at \(\mathbf {o}({\boldsymbol \gamma })\). Since \({\boldsymbol \gamma }^{\prime }\) has the same sign as \({\boldsymbol \gamma }\), we know that \(\Omega\) is also drained at \(\mathbf {o}({\boldsymbol \gamma }^{\prime })\), so \(\mathbf {o}({\boldsymbol \gamma }^{\prime })\) is also output-stable. We must then have that \(f(\mathbf {x}) = \mathbf {o}({\boldsymbol \gamma }^{\prime })_Y\), but by the construction of \({\boldsymbol \gamma }^{\prime }\), we know that \(\mathbf {o}({\boldsymbol \gamma }^{\prime })_Y \in \mathbb {Q}\). Since \(B \subseteq \Sigma (\Omega)\), we know that \(f|_B\) is linear by Lemma 5.23. Since B is an open ball, we know that \(\mathbb {Q}^\Lambda \cap B\) contains a basis for \(\mathbb {Q}^\Lambda\), so \(f|_B\) is a linear function that maps \(\mathbb {Q}^\Lambda\) to \(\mathbb {Q}\). Since every linear function \(\mathbb {R}^\Lambda \rightarrow \mathbb {R}\) that sends \(\mathbb {Q}^\Lambda\) to \(\mathbb {Q}\) is rational linear, we are done.□
Recall that a closed domain is the closure of an open set.
Lemma 5.37.
Let \(f: X \rightarrow Y\) be a continuous function defined piecewise on closed sets, so X is covered by finitely many closed sets \(D_1 \ldots D_k\) and there are continuous functions \(g_i: X\rightarrow Y\) such that \(f|_{D_i} = g_i|_{D_i}\). Then there are (possibly empty) closed domains \(E_1 \ldots E_k\) that cover X such that \(f|_{E_i} = g_i|_{E_i}\).
Proof.
We show how to convert each \(D_i\) that is not a closed domain to a corresponding \(E_i\) that is. Let \(D_i\) be some set that is not a closed domain. Let
Note that \(D_1, \ldots D_{i - 1}, X \setminus D_+, D_{i + 1}, \ldots D_k\) cover X. Let \(E_i\) be the closure of \(X \setminus D_+\). Clearly, \(E_i\) is a closed domain. Since our original sets \(D_1 \ldots D_k\) cover X, we know that \(X \setminus D_+ \subseteq D_i\), and, since \(D_i\) is closed, this implies that \(E_i \subseteq D_i\). Because of this, we also have that \(f|_{E_i} = g_i|_{E_i}\). Finally, because \(X \setminus D_+ \subseteq E_i\), we know that the sets \(D_1, \ldots D_{i - 1}, E_i, D_{i + 1}, \ldots D_k\) cover X.□
The following is the main result of Section 5.5, showing a limitation on the computational power of CRCs stably computing functions in the direct sense.
Lemma 5.38.
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)\). Then f is positive-continuous and piecewise rational linear.
Proof.
By Lemma 5.32, we know that f is positive-continuous and piecewise linear. By a general topological argument one could show that any function with these properties has domains of definition that are closed relative to \(\mathbb {R}_{\gt 0}^\Sigma\), but, since by Lemma 5.29 the domains we constructed earlier already have this property, we will not give the general proof here. By Lemma 5.37, we can replace the closed sets that give the domains of definition of f by closed domains. If some of the domains produced by Lemma 5.37 are empty, then we can simply ignore them in what follows. Since all of the nonempty domains are the closures of nonempty open sets, they must each contain some open ball. By Lemma 5.36, f is a rational linear function when restricted to each of these closed domains, so f is piecewise rational linear on \(\mathbb {R}_{\gt 0}^\Sigma\).
For any proper subset U of the input species, one can apply the above argument to the reduced CRN that discards all species not producible from the given inputs to show that f is continuous and piecewise rational linear on \(D_U\). This shows that f is a positive-continuous piecewise rational linear function on all of \(\mathbb {R}_{\ge 0}^\Sigma\).□
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.
Proposition 5.39.
Let \(f:\mathbb {R}^k \rightarrow \mathbb {R}\) be stably dual computed by a CRC. Let \(\Omega\) be a stabilizing siphon. Then f restricted to inputs that have a dual rail representation in \(\Sigma (\Omega)\) is linear.
Proof.
A dual-rail computing CRC can be thought to directly compute two separate functions \(\hat{f}^+, \hat{f}^-: \mathbb {R}_{\ge 0}^{2k} \rightarrow \mathbb {R}_{\ge 0}\) such that \(\hat{f} = \hat{f}^+ - \hat{f}^-\) where \(\hat{f}\) is a dual rail representation of f. By Lemma 5.23, we know that \(\hat{f}^+\) and \(\hat{f}^-\) are rational linear when restricted to \(\Sigma (\Omega)\). The proposition follows, because linearity is closed under subtraction.□
The following is our main negative result for dual-rail CRC’s, a dual-rail analog of Lemma 5.38:
Lemma 5.40.
Let \(f:\mathbb {R}^k \rightarrow \mathbb {R}\) be stably dual-computable by a CRC. Then f is continuous and piecewise rational linear.
Proof.
Let \(\mathcal {C}\) be the CRC stably computing a dual-rail representation \(\hat{f}\) of f, with input species \(X_1^+\), \(\ldots\), \(X_k^+\), \(X_1^-\), \(\ldots\), \(X_k^-\) and output species \(Y^+,Y^-\).
Similarly to the proof of Proposition 5.39, a dual-rail computing CRC can be thought to directly compute two separate functions \(\hat{f}^+, \hat{f}^-: \mathbb {R}_{\ge 0}^{2k} \rightarrow \mathbb {R}_{\ge 0}\) such that \(\hat{f} = \hat{f}^+ - \hat{f}^-\) where \(\hat{f}\) is a dual rail representation of f. Since \(\hat{f}^+\) and \(\hat{f}^-\) are stably computed by a CRC, Lemma 5.38 implies that they are both piecewise rational linear, and, since piecewise rational linear functions are closed under subtraction, this implies that f is also piecewise rational linear. It remains to show that f is continuous.
For any input \(\mathbf {x}^{\prime } \in \mathbb {R}^k\) to f, there is an initial state \(\mathbf {x}\in \mathbb {R}^\Sigma _{\gt 0}\) representing \(\mathbf {x}^{\prime }\), with strictly positive concentrations of all input species. (For example, if \(\mathbf {x}^{\prime }(1) = 5\), then we can choose \(\mathbf {x}(X^+_1) = 6\) and \(\mathbf {x}(X^-_1) = 1\).) Let \(\mathbf {x}^{\prime }_1,\mathbf {x}^{\prime }_2,\ldots \in \mathbb {R}^k\) be any sequence of inputs to f such that \(\lim _{i\rightarrow \infty } \mathbf {x}^{\prime }_i = \mathbf {x}^{\prime }\). Let us represent each input \(\mathbf {x}^{\prime }_i\) in the sequence by an initial state \(\mathbf {x}_i \in \mathbb {R}^\Sigma _{\gt 0}\) such that \(\lim _{i\rightarrow \infty } \mathbf {x}_i = \mathbf {x}\). Then the sequence of initial states has the property that all of the \(\mathbf {x}_i\)’s obey \([\mathbf {x}_i]=\Sigma\). By Lemma 5.32, f is continuous on the domain in which all input species are positive, which includes the input represented by \(\mathbf {x}\) and the inputs represented by all of the \(\mathbf {x}_i\)’s. Therefore, \(f(\mathbf {x}^{\prime }) = \lim _{i\rightarrow \infty } f(\mathbf {x}^{\prime }_i)\), so f is continuous.□
6 Extensions
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].
Definition 6.1.
An infinite game is a pair of a set X of possible moves and a set \(A \subseteq X^\mathbb {N}\) called the payoff set.
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.
Definition 6.2.
For an infinite game \((X, A)\) a winning strategy for player II is a collection of functions \(f_{2k}: X^k \rightarrow X\) so the sequence
is never in A for every choice of \(x_1, x_3, \ldots\). A winning strategy for player I is defined similarly.
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.
Definition 6.3.
For a CRC \(\mathcal {C}\) and function \(f: \mathbb {R}_{\ge 0}^k \rightarrow \mathbb {R}_{\ge 0}\), an initial state \(\mathbf {x}\in \mathbb {R}_{\ge 0}^\Sigma\), the Strong Demon Strong Chemist game is the game where the moveset \(M(\mathbf {f}_1, \ldots , \mathbf {f}_n)\) includes all the moves \(\mathbf {f}_{n+1}\) such that \(\mathbf {f}_1 \circ \cdots \circ \mathbf {f}_{n+1}\) is a valid rate schedule starting at \(\mathbf {x}\), and a sequence of moves \(\left(\mathbf {f}_1, \mathbf {f}_2, \ldots \right)\) consistent with the moveset is in the payoff set iff \(\lim _{t \rightarrow \infty } {\boldsymbol \rho }(t) \upharpoonright \left\lbrace Y\right\rbrace \ne f(\mathbf {x})\) for trajectory \({\boldsymbol \rho }\) corresponding to \(\mathbf {f}_1 \circ \mathbf {f}_2 \circ \ldots\) starting at \(\mathbf {x}\).
Lemma 6.4.
The following are equivalent for a CRC \(\mathcal {C}\) and function \(f: \mathbb {R}_{\ge 0}^k \rightarrow \mathbb {R}_{\ge 0}\):
(1)
The Chemist has a winning strategy for the Strong Demon Strong Chemist game associated with \(\mathcal {C}\) and f for every initial state \(\mathbf {x}\in \mathbb {R}_{\ge 0}^\Sigma\).
(2)
\(\mathcal {C}\) stably computes f.
Proof.
First suppose that \(\mathcal {C}\) stably computes. Then, by Theorem 2.27, we know that for any first move \(\mathbf {f}_1\) the Demon plays, the state \({\boldsymbol \rho }(1)\) will be reachable from the initial state \(\mathbf {x}\). By the definition of stable computation there must be a state \(\mathbf {z}\) reachable from \(\mathbf {x}\) so \(\mathbf {z}\) is output-stable with the correct amount of output. By Lemma 2.26, there is a valid rate schedule \(\mathbf {f}_2\) that goes from \(\mathbf {x}\) to \(\mathbf {z}\), and by “rescaling” the rate schedule by increasing the rate of every reaction, we can guarantee that \(\mathbf {f}_2\) is completed within a single unit of time. Since \(\mathbf {z}\) is output-stable, every further move \(\mathbf {f}_4, \mathbf {f}_6, \ldots\) that the Chemist makes can just be taken to be the zero rate schedule (where no reaction occurs). This specifies a winning strategy for the Chemist for the Strong Demon Strong Chemist game associated with \(\mathcal {C}\), f, and \(\mathbf {x}\).
However, suppose that \(\mathcal {C}\) does not stably compute. Then let \(\epsilon\), \(\mathbf {x}\), and \(\mathbf {z}\) be the ones given by Theorem 3.3. The Demon can play a rate schedule that goes from \(\mathbf {x}\) to \(\mathbf {z}\). After this, every time it is the Demon’s turn the CRC will be in a state \(\mathbf {o}\) that is reachable from \(\mathbf {z}\), so the Demon will always be able to play a valid rate schedule that takes the CRN to \(\mathbf {o}^{\prime }\) as in Theorem 3.3 with \(|\mathbf {o}^{\prime }(Y) - f(\mathbf {x})| \gt \epsilon\). This shows that the Demon has a winning strategy for the Strong Demon Strong Chemist game associated with \(\mathcal {C}\), f and \(\mathbf {x}\), so certainly the Chemist does not have a winning strategy.□
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.
Definition 6.5.
A moveset \(M(\mathbf {f}_1, \ldots , \mathbf {f}_n)\) is valid if it includes only moves \(\mathbf {f}_{n+1}\) such that \(\mathbf {f}_1 \circ \cdots \circ \mathbf {f}_{n+1}\) is a valid rate schedule starting at \(\mathbf {x}\). Further, we say that a valid moveset \(M(\mathbf {f}_1, \ldots , \mathbf {f}_n)\) is Chemist-fair if there is some \(\mathbf {H}: \mathbb {R}_{\ge 0}^\Lambda \rightarrow \mathbb {R}_{\ge 0}^R\) (with the same properties as in Definition 4.1) such that for any odd n, the moveset contains at least one \(\mathbf {f}_{n+1}\) such that the trajectory \({\boldsymbol \rho }\) corresponding to \(\mathbf {f}= \mathbf {f}_1 \circ \cdots \circ \mathbf {f}_{n+1}\) starting at \(\mathbf {x}\) satisfies \(\mathbf {f}_\alpha (t) \ge \mathbf {H}_\alpha ({\boldsymbol \rho }(t))\) for every reaction \(\alpha\) and for all \(t \in [n,n+1)\).
Definition 6.6.
Consider a CRC \(\mathcal {C}\), function \(f: \mathbb {R}_{\ge 0}^k \rightarrow \mathbb {R}_{\ge 0}\), and an initial state \(\mathbf {x}\in \mathbb {R}_{\ge 0}^\Sigma\). Let \(M(\mathbf {f}_1, \ldots , \mathbf {f}_n)\) be any valid moveset. We say the Strong Demon M-Chemist is the game where a sequence of moves \(\left(\mathbf {f}_1, \mathbf {f}_2, \ldots \right)\) consistent with the moveset is in the payoff set iff \(\lim _{t \rightarrow \infty } {\boldsymbol \rho }(t) \upharpoonright \left\lbrace Y\right\rbrace \ne f(\mathbf {x})\) for trajectory \({\boldsymbol \rho }\) corresponding to \(\mathbf {f}_1 \circ \mathbf {f}_2 \circ \cdots\) starting at \(\mathbf {x}\).
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:
Definition 6.7.
A Strong Demon complexity class\(\mathsf {SD}\) is a set of functions \(f: \mathbb {R}_{\ge 0}^k \rightarrow \mathbb {R}_{\ge 0}\) for all \(k \in \mathbb {N}\), that is characterized in the following way: There is a (non-empty) set of Chemist-fair movesets \(\mathcal {M}\) so \(f: \mathbb {R}_{\ge 0}^k \rightarrow \mathbb {R}_{\ge 0}\) is in \(\mathsf {SD}\) if and only if there is a CRC \(\mathcal {C}\) such that for each \(M \in \mathcal {M}\) and initial state \(\mathbf {x}\) the Chemist has a winning strategy for the associated Strong Demon M-Chemist game.
Theorem 6.8.
Every Strong Demon complexity class is the same as the Strong Demon Strong Chemist class.
Proof.
If a Chemist has a winning strategy in a Strong Demon M-Chemist game, then clearly the Chemist has a winning strategy in the Strong Demon Strong Chemist game. Thus, every Strong Demon complexity class is a subset of the Strong Demon Strong Chemist class. For the other direction, let \(\mathsf {SD}\) be some Strong Demon complexity class, and let \(\mathcal {M}\) be the nonempty set of Chemist-fair movesets associated with \(\mathsf {SD}\). Suppose that \(f: \mathbb {R}_{\ge 0}^k \rightarrow \mathbb {R}_{\ge 0}\) is in the Strong Demon Strong Chemist class. Then, by Lemma 6.4 and Theorem 5.9, we know that f is fairly computable by a CRC \(\mathcal {C}\). For any moveset \(M \in \mathcal {M}\), let us show that the Chemist has a winning strategy for the Strong Demon M-Chemist game associated with \(\mathcal {C}\), and f. Since M is Chemist-fair, we know that there is some function \(\mathbf {H}: \mathbb {R}_{\ge 0}^\Lambda \rightarrow \mathbb {R}_{\ge 0}^R\) as in Definition 6.5. On their turn, the Chemist simply plays any move in the moveset M that satisfies the condition on \(f_{n+1}\) in definition 6.5. This strategy is winning for the Chemist, because the total rate schedule \(\mathbf {f}\) obtained by concatenating the \(\mathbf {f}_i\)’s is fair, since on the set \(T = \bigcup _{n = 1}^\infty [2n-1, 2n)\), we know that \(\mathbf {f}_\alpha (t) \ge \mathbf {H}_\alpha ({\boldsymbol \rho }(t))\). Since \(\mathcal {C}\) was chosen to fairly compute f this proves that \(\lim _{t \rightarrow \infty } ({\boldsymbol \rho }(t)\upharpoonright \left\lbrace Y\right\rbrace) = f(\mathbf {x})\). This shows that \(f \in \mathsf {SD}\), so any Strong Demon M-Chemist complexity class is the same as the Strong Demon Strong Chemist complexity class.□
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.23 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
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.
7 Conclusion and Open Questions
We characterized the class of functions computable in a manner that is absolutely robust to reaction rates in the continuous model of chemical kinetics. Such rate-independent computation must rely solely on reaction stoichiometry—which reactants, and how many of each, become which products, and how many of each? We considered two methods of encoding inputs and outputs: direct and dual-rail. The dual-rail encoding permits easier composition of modules and can represent negative values; we characterized its computational power as continuous, piecewise rational linear. The direct encoding, however, allows computing functions that are discontinuous at the faces of the nonnegative orthant. For both encodings, we showed matching negative results (showing that nothing more can be computed) and positive results (describing CRNs computing any function in the class).
Since rate-independent computation does not require difficult-to-achieve tuning of parameters or reaction conditions, it may be significantly more “engineerable” than rate-dependent computation. More generally, our work also helps uncover the multifaceted sources of chemical computational power by disentangling the control of stoichiometry from reaction rates.
We now describe some natural open questions.
Reaction complexity of stably computable functions. An interesting question regards the description complexity of functions stably computable by CRNs. Some piecewise linear functions have a number of pieces exponential in the number of inputs; for example, \(f(x_1,\ldots ,x_{2k}) = \min (x_1,x_2) + \min (x_3,x_4) + \cdots + \min (x_{2k-1}, x_{2k})\) has \(2^k\) linear pieces. If we express this function in \(\max \min g_i\) form of Theorem 5.11, then we need \(2^k\) different linear \(g_i\), and thus the construction in the proof of Lemma 5.15 would require exponentially many species and reactions. However, this particular f has a more succinct CRN that stably computes it, namely, the reactions
\[\begin{eqnarray*} X_1 + X_2 &\rightarrow & Y \\ X_3 + X_4 &\rightarrow & Y \\ \vdots \\ X_{2k-1} + X_{2k} &\rightarrow & Y. \end{eqnarray*}\]
Given a positive-continuous, piecewise rational linear function f, how can we tell whether it has a more compact CRN stably computing it than our construction? If it does, then how can we arrive at it?
Requiring “always” fair rate schedules. Lemma 4.9 shows that for a feedforward CRC that stably computes a function f, any fair rate schedule converges to the correct output of f, where fair essentially means that applicable reactions must have positive rate for an infinite subset of times. In other words, the adversary is allowed for some times outside this subset for the rate schedule to be unfair: to “starve” some applicable reactions by keeping their rates at 0 despite all reactants being present. Since such rate schedules seem physically implausible, it is natural to consider a modified definition of computation, one that requires every rate schedule that is always-fair (i.e., for all time, applicable reactions must have positive rate) to converge to the correct output. The following CRC shows that these two requirements can result in different behaviors for a given CRC:
\[\begin{eqnarray*} X &\rightarrow & X+C \\ X &\rightarrow & X^{\prime } \\ C+X^{\prime } &\rightarrow & C+Y. \end{eqnarray*}\]
With initial concentration x of X, every always-fair rate schedule converges to concentration x of Y, so it computes the identity function \(f(x)=x\) under this modified definition. (Note that always-fairness requires the first two reactions to have positive rate at time 0, so C immediately becomes present, at which point it is inevitable to convert all X to \(X^{\prime }\) and all \(X^{\prime }\) to Y.) However, the CRC does not stably compute f: An initially unfair adversary can execute the second reaction to completion, starving the first reaction until it becomes inapplicable. This schedule never produces any C, so Y stays at 0, yet the schedule is fair by the original definition, since, once all reactions become inapplicable, subsequent rate 0 of all reactions for all time vacuously satisfies the definition of fair. It is an interesting question is whether, under the modified definition of always-fair, some CRC can compute a function that is not stably computable.
Arbitrary but fixed rate constants. A related notion of rate-independence is one where the form of rate-law cannot vary, but the constant parameters can, e.g., mass-action rates, where an adversary picks the rate constants (possibly depending on the initial input). For example, consider the following CRN with input species \(A,B,C,X\) and output species Y.
\[\begin{eqnarray*} A + X &\rightarrow & A + Y \\ B + Y &\rightarrow & B + X \\ A + B &\rightarrow & C \\ C + Y &\rightarrow & C + X \\ 3C &\rightarrow & \emptyset . \end{eqnarray*}\]
Let \(a, b, c, x\) denote the initial concentrations of species A, B, C, X. This system does not stably compute any function in the model defined in this article because on input \(a=b\), it can stabilize to any value of output y between 0 and x.
In contrast, consider the above system under mass-action kinetics, where an adversary picks the rate constants, but they remain constant over time. Because A and B are required to produce C, and at least one of them goes to 0 by the third reaction, the concentration of C approaches 0 as time goes to infinity by the final reaction, no matter the rate constants. If \(a\gt b\), then also the concentration of B approaches 0 but the concentration of A remains bounded away from 0. Therefore, the output y converges to x, regardless of what the rate constants are. Similarly, if \(b\gt a\), then the output y approaches 0. When \(a=b\), the concentrations of \(A, B, C\) approach 0 at different rates. The concentrations of \(A, B\) are \(\Theta (\frac{1}{t})\) at time t (rate of bimolecular decay), and the concentration of C is \(\Theta (\frac{1}{\sqrt {t}})\) (rate of trimolecular decay). As a result, the effective rate of conversion of Y to X via the channel \(C + Y \rightarrow C + X\) is \(\Theta (\frac{1}{\sqrt {t}})\). Since this is \(\omega (\frac{1}{t})\) the output y always converges to 0 regardless of the rate constants (in our particular case, the concentration of Y is \(e^{-\Theta (\sqrt {t})}\)). From the above, this CRN computes \(f(a, b, c, x) = x\) when \(a \gt b\) and \(f(a, b, c, x) = 0\) when \(a \le b\), no matter what the rate constants are. This function is discontinuous at points where \(a=b\), so it is not positive-continuous, thus not stably computable by any CRN under our model of rate-independence. It remains open to classify what functions can be computed by mass-action CRNs in which rate constants are chosen adversarially.
Axiomatic derivation of reachability. An important contribution of this article is to develop segment-reachability as the “correct” notion of rate-independent reachability. While we justify our definition of segment-reachability by its equivalence to valid rate schedules (Theorem 2.27), one could imagine taking an axiomatic approach: Any notion of rate-independent reachability ought to satisfy certain constraints. For instance, it should be transitive. Also, if \(\mathbf {d}\) is straight-line reachable from \(\mathbf {c}\), then it is evidently allowed by stoichiometry for the CRN to evolve from \(\mathbf {c}\) to \(\mathbf {d}\), so \(\mathbf {d}\) should be reachable from \(\mathbf {c}\). Of course, there are many relations that satisfy these two constraints: For instance, the relation where every state is reachable from every other state satisfies these two constraints—but this is evidently too permissive. Indeed, our notion of segment-reachability is the transitive closure of the relation of straight-line reachability (Corollary 2.16), so it is minimal among all of the relations with these two properties. However, since our goal is to develop a notion of reachability that is as unrestricted as possible, a more compelling axiomatic derivation would follow from simple set of natural conditions for which segment-reachability is maximal.
Absolute inhibition. In biology it is not uncommon to have rate laws with explicit inhibitors; the higher the concentration of the inhibitor, the slower the reaction. Formally, general rate laws have been described in which reactants are partitioned into consumed species, species that increase the reaction rate (catalysts), and species that decrease the rate (inhibitors) [26]. In our model, while inhibitors can be modeled mechanistically as sequestering reacting species in a non-reactive form (e.g., \(A + B \rightarrow C\) is inhibited by I via the reaction \(A + I \mathop {\rightleftharpoons }\limits AI\)), such inhibition just slows down the reaction rate but does not prevent the reaction entirely. In contrast, one can imagine a definition of reachability in which a reaction is applicable only if all of its reactants are present and all of its inhibitors are absent. It remains an open question whether the computational power of rate-independent computation changes if such absolute inhibition is allowed. We note that this change drastically changes the notion of reachability because it is no longer additive: It might be that \(\mathbf {x}\rightsquigarrow \mathbf {y}\), but \(\mathbf {x}+ \mathbf {c}\not\rightsquigarrow \mathbf {y}+ \mathbf {c}\) if \(\mathbf {c}\) contains some inhibitors of reactions occurring in the first path. Further, for discrete CRNs, absolute inhibition dramatically expands computational power of stable computation to Turing universality [21, 38] (i.e., the ability to compute any function computable by any algorithm). thus it seems reasonable to conjecture that it also expands the computational power in the continuous setting.
Decision problems. In the discrete CRN model, particularly the subset of it known as population protocols [6], a major focus of research is on decision problems with a yes/no output, a.k.a. predicates. The typical output convention partitions the set \(\Lambda\) of species into two disjoint subsets \(\Lambda = \Lambda _\mathsf {Y}\cup \Lambda _\mathsf {N}\), the “yes voters” and “no voters.” The goal of stable computation in this setting is to reach a configuration with a unanimous, correct vote (e.g., if the correct answer is yes, then only species in \(\Lambda _\mathsf {Y}\) are present), which is also stable in the sense that no incorrect voter is producible. In this setting, it has been shown that exactly the semilinear predicates can be stably decided [6, 7] by discrete CRNs.
The concept can be similarly defined with continuous CRNs using our notion of segment-reachability. Say that a CRN with voting species defined as above stably computes a predicate \(\phi : \mathbb {R}_{\ge 0}^k \rightarrow \lbrace \mathsf {Y},\mathsf {N}\rbrace\) if, for any initial configuration \(\mathbf {x}\in \mathbb {R}_{\ge 0}^k\), for any configuration \(\mathbf {c}\) such that \(\mathbf {x}\rightsquigarrow \mathbf {c}\), there is a configuration \(\mathbf {y}\) such that \(\mathbf {c}\rightsquigarrow \mathbf {y}\), and \(\mathbf {y}\) is stably correct, meaning that for all \(\mathbf {y}^{\prime }\) such that \(\mathbf {y}\rightsquigarrow \mathbf {y}^{\prime }\), we have \(\emptyset \ne [\mathbf {y}^{\prime }] \subseteq \Lambda _{\phi (\mathbf {x})}\). In other words, stable computation leads to a nonempty configuration with only “correct” votes (according to \(\phi\)), and this is also true of every configuration reachable from there. We also say in this case that the CRN stably decides the set \(\phi ^{-1}(\mathsf {Y}) \subseteq \mathbb {R}_{\ge 0}^k\) of inputs that map to output \(\mathsf {Y}\). What is the class of sets \(S \subseteq \mathbb {R}_{\ge 0}^k\) that are stably decidable by this definition?
For discrete CRNs, the question of “how long” a system takes to stably compute a function has received much interest [1, 3, 25, 28]. In general, asking questions of time-complexity of continuous computation seems more difficult than in the discrete setting. For general polynomial ODEs, the breakthrough work of Bournez, Graça, and Pouly [12] established a surprisingly tight connection between the length of the trajectory and the Turing machine computation time. Closer to our domain of interest, prior work has studied asymptotic convergence speed for the composition of simple CRN motifs as a function of the number of (feedforward) layers [46]. Although not explicitly stated in terms of rate-independent computation, these modules compute in the rate-independent manner as studied here. It is interesting to ask whether these techniques could be adopted to our constructions to articulate and help resolve questions of computation time.
Acknowledgments
We thank Manoj Gopalkrishnan, Elisa Franco, Damien Woods, and the organizers and participants of the American Mathematical Institute workshop on Mathematical Problems Arising from Biochemical Reaction Networks for insightful discussions. We are grateful to anonymous reviewers for insightful comments and suggestions that have greatly improved this article.
Footnotes
1
Although the finite density of matter physically restricts what the largest concentration of any species could realistically be, standard models of chemical kinetics focus on systems that are far from this bound, mathematically allowing concentrations to be arbitrarily large.
2
A species acts catalytically in a reaction if it is both a reactant and product: e.g., C in reaction \(A + C \rightarrow B + C\). Note that executing this reaction without C does not by itself violate condition (1).
3
See Section 2.4 for examples showing that in the continuous setting conditions (2) and (3) are not mutually redundant.
4
The exact statement of Kurtz’s convergence result [34] is beyond the scope of this article. It considers taking a discrete CRN with initial integer molecular counts given by vector \(\mathbf {c}\in \mathbb {N}^{k}\) in volume \(V \gt 0\), then “scaling up” by factor \(n \in \mathbb {N}\), i.e., considering the discrete CRN with initial state \(n \cdot \mathbf {c}\) in volume \(n \cdot V\). The result, stated very roughly, is that with high probability the n-scaled CRN has a trajectory (dividing discrete counts by \(n \cdot V\) to convert to units of concentration) that stays close to the real-valued mass-action concentration trajectory, but only for time \(O(\log n)\). An example CRN where this time bound is tight is \(A+B \rightarrow A+B+Y,\quad A \rightarrow 2A,\quad B \rightarrow \emptyset\), starting with \(1A, 1B\). In mass-action, the concentrations of A and B at time t are, respectively, \(e^t\) and \(e^{-t}\), whose product is the constant 1, so the first reaction produces Y at a unit rate forever. However, scaling up to \(n A, n B\), the discrete CRN consumes all B in \(O(\log n)\) time, at which point production of Y halts. See References [24, 35] for other example CRNs for which the two models diverge after sufficient time.
5
It is generally supposed that chemical reactions would follow mass-action if properly decomposed into truly elementary reactions and the solution is well-mixed. For example, Michaelis-Menten and Hill-function kinetics can be derived as a limiting case of mass-action when the reaction is initiated and completed at vastly different time scales.
6
For example, consider the reaction \(A \rightarrow 2B\), with ODEs \(\dot{a} = -a\) and \(\dot{b} = 2a\). One can imagine a “chemical” adversary adjusting the rate of the reaction to speed it up or slow it down, but what the adversary cannot control is that to consume x amount of A requires producing exactly 2x amount of B and vice versa. This connection between the rates of consumption of A and production of B does not have an obvious counterpart in more general polynomial ODEs and analog computational models.
7
Our notion of valid rate schedules in definition 2.22 is even more general than a reaction system in that a valid rate schedule does not require a reaction’s rate to be a function of species concentrations, for instance, allowing an adversary to visit the same state twice but apply different reaction rates each time.
8
It is customary to define, for each reaction, a rate constant\(k \in \mathbb {R}_{\gt 0}\) specifying a constant multiplier on the mass-action rate (i.e., the product of the reactant concentrations), but as we are studying CRNs whose output is independent of the reaction rates, we leave the rate constants out of the definition.
9
We allow high-order reactions; i.e., those that have more than two reactants. Such higher order reactions could be eliminated from our constructions using the transformation that replaces \(S_1 + S_2 + \cdots + S_n \rightarrow P_1 + \cdots + P_m\) with bimolecular reactions \(S_1 + S_2 \mathop {\rightleftharpoons }\limits S_{12}, S_{12} + S_3 \mathop {\rightleftharpoons }\limits S_{123}, S_{123} + S_4 \mathop {\rightleftharpoons }\limits S_{1234}, \ldots , S_n + S_{12 \ldots n-1} \rightarrow P_1 + \cdots + P_m\).
10
Note that \(\mathbf {M}\) does not fully specify \(\mathcal {C}\), since catalysts are not modeled: reactions \(Z + X \rightarrow Z + Y\) and \(X \rightarrow Y\) both correspond to the column vector \((-1,1,0)^\top\).
11
Another property of \(\rightsquigarrow\) that we use extensively is scale-invariance: If \(\mathbf {x}\rightsquigarrow \mathbf {y}\), then \(\lambda \mathbf {x}\rightsquigarrow \lambda \mathbf {y}\) for any \(\lambda \ge 0\), which is essentially responsible for the convexity of Lemma 2.7. This does not hold for discrete CRN reachability when \(\lambda \lt 1\), even when the scaled discrete states are well-defined, e.g., the reaction \(X+X \rightarrow Y\) is applicable in state \(\mathbf {x}= \lbrace 2X\rbrace\) but not in state \(0.5 \mathbf {x}= \lbrace 1 X\rbrace\) in the discrete model.
12
Indeed, the mass-action ODE corresponding to the CRN \(2X \rightarrow 3X\) is \(\frac{dx}{dt} = x^2\), which is solved by \(x(t) = \frac{1}{C - t}\), where \(C = 1/x(0)\). This goes to infinity as t approaches C.
13
Note that a more general definition would say \(\mathbf {d}\) is mass-action reachable from \(\mathbf {c}\) if there exist positive rate constants such that the trajectory starting at \(\mathbf {c}\) passes through or approaches \(\mathbf {d}\). Note, however, that this relation is not transitive: For some CRNs, \(\mathbf {c}\) reaches to \(\mathbf {d}\) under one set of rate constants, and \(\mathbf {d}\) reaches to \(\mathbf {x}\) under another set of rate constants, yet no single assignment of rate constants takes the CRN from \(\mathbf {c}\) to \(\mathbf {x}\).
14
It is obvious in the discrete CRN model, and in an intuitive physical sense, that if producing a species initially absent causally requires another species also initially absent and vice versa, then neither species can ever be produced. However, it requires care to prove this for mass-action ODEs. Consider the CRN \(2X \rightarrow 3X\). The corresponding mass-action ODE is \(dx/dt = x^{2}\), and has the property that starting with \(x(0) = 0\), it cannot become positive, i.e., the only solution with \(x(0)=0\) is \(x(t) = 0\) for all \(t \ge 0\). However, the very similar non-mass-action ODE \(dx/dt = x^{1/2}\) has a perfectly valid solution \(x(t) = t^2/4\), which starts at 0 but becomes positive, despite the fact that at \(t=0\), \(dx/dt = 0\). (Though \(x(t) = 0\) for all \(t \ge 0\) is another valid solution.) The difference is that mass-action polynomial rates are locally Lipschitz (have bounded rates of change, unlike \(x^{1/2}\), whose derivative goes to \(\infty\) as \(x \rightarrow 0\)) and so are guaranteed to have a unique solution by the Picard-Lindelöf theorem.
15
An alternative to Definition 2.22 would start with a differentiable trajectory \({\boldsymbol \rho }\) and total flux \(\mathbf {F}\) (related via \({\boldsymbol \rho }(t) = \mathbf {M}\cdot \mathbf {F}(t) + \mathbf {c}\)) and define \(\mathbf {f}= d\mathbf {F}/dt.\) However, requiring differentiable \({\boldsymbol \rho }\) and \(\mathbf {F}\) rules out many natural cases, such as the rate schedules implicit in segment-reachability (Definition 2.2), whose trajectories are not differentiable at cusp points \(\mathbf {b}_i\) in between straight lines and whose rate schedules are not even continuous.
16
Consider a rate schedule that takes concentrations negative: for instance, starting with \(1 X\) and applying reaction \(\alpha : X \rightarrow Y\) with \(\mathbf {F}_\alpha (t) = 2\) for some \(t \gt 0\). To see that this contradicts condition (2) when we naturally generalize notation \([\mathbf {c}]\) to possibly negative \(\mathbf {c}\) (for any \(\mathbf {c}\in \mathbb {R}^\Lambda\), \([\mathbf {c}] = \lbrace S \in \Lambda \ |\ \mathbf {c}(S) \gt 0 \rbrace\)), suppose that \({\boldsymbol \rho }_X(t^{\prime }) \lt 0\) for some species X at some time \(t^{\prime }\). Let \(t_0\) be the supremum of all the times less than \(t^{\prime }\) where \({\boldsymbol \rho }_X(t) \ge 0\). Recall \(\mathbf {F}_\alpha (t)\) is a locally absolutely continuous (and therefore continuous) function. Thus, \({\boldsymbol \rho }(t)\) is also a continuous function so \({\boldsymbol \rho }_X(t_0) \ge 0\). Moreover, for all \(t_0 \lt t \lt t^{\prime }\), we know that \({\boldsymbol \rho }_X(t) \lt 0\) by our choice of \(t_0\). So, by condition (2) \(\mathbf {f}_\alpha (t) = 0\) for all \(\alpha\) where X is a reactant, and therefore (recall \(\mathbf {M}(X,\alpha)\) means the net consumption of X in reaction \(\alpha\))
a contradiction, since \({\boldsymbol \rho }(t^{\prime })\) was assumed to be negative. See also Reference [26, Proposition 2.8], where the term strict is equivalent to condition (2).
17
Although Lemma 2.24 has the precondition that the mass-action trajectory be defined for all time, states reached in finite time by diverging mass-action CRNs can also be segment-reached. For example, for the CRN \(2X \rightarrow 3X\) (with rate constant 1) starting in \(\lbrace 1 X\rbrace\), which diverges as \(t \rightarrow 1\), all states on the trajectory prior to time \(t = 1\) are segment reachable: For each such state, we can construct a valid rate schedule that obeys mass-action until reaching that state and then is constant for all later time.
18
We assume a canonical ordering of \(\Sigma =\lbrace X_1,\ldots ,X_k\rbrace\) so a vector \(\mathbf {x}\in \mathbb {R}_{\ge 0}^k\) (i.e., an input to f) can be viewed equivalently as a state \(\mathbf {x}\in \mathbb {R}_{\ge 0}^\Sigma\) of \(\mathcal {C}\) (i.e., an input to \(\mathcal {C}\)). Note that we have defined valid initial states to contain only the input species \(\Sigma\); other species must have initial concentration 0. Our results would change slightly if we relaxed this assumption—see Section 6.2.
19
The discrepancy between stable computation and correctness under mass-action kinetics shows a major difference between the discrete and continuous CRN models. In the example above, with n total molecules, the discrete CRN model does a random walk on the number of X that is biased toward the dynamic equilibrium point \(n k_2 / (2k_1 + k_2)\). Despite the bias upward when X is below this value, there is always a positive probability to decrease X, so with probability 1, X will eventually reach 0.
20
By Lebesgue’s fundamental theorem of calculus [43, Theorem 6.11, Theorem 6.14], applied to \(\mathbf {f}_\alpha (t)\), we know that \({\boldsymbol \rho }(t)\) is locally absolutely continuous, and almost everywhere differentiable with derivative \(\mathbf {M}\cdot \mathbf {f}(t)\).
21
Note that \(\int _0^\infty \frac{d}{dt} V({\boldsymbol \rho }(t)) dt = \lim _{t \rightarrow \infty } V({\boldsymbol \rho }(t)) - V(\mathbf {c})\) by the fundamental theorem of calculus [43, Theorem 6.10]. To use the fundamental theorem of calculus, we must ensure that \(V({\boldsymbol \rho }(t))\) is locally absolutely continuous—this follows from the fact that the trajectory \({\boldsymbol \rho }(t)\) is defined in terms of a (Lebesgue) integral.
22
This analysis of the CRC for the min function here is directly based on the definition of stable computation. Recently a powerful framework has been developed [52], based on a wide class of so-called noncompetitive CRCs in which no species consumed in a reaction is a reactant in another reaction (not even as a non-consumed catalyst). For such CRCs, the task of proving correctness of stable computation is greatly simplified. Since the CRC computing min is noncompetitive, that framework could be applied here to yield a simpler proof of correctness. We use our direct proof here for the sake of making the current article self-contained.
23
For example, to simulate initial context \(\lbrace 3/2\ Z_1, 1/2\ Z_2 \rbrace\) from \(\lbrace 1\ Z\rbrace\), add the reactions \(2Z \rightarrow Z^{\prime }\) and \(Z^{\prime } \rightarrow 3 Z_1 + Z_2\).
Appendices
A Forward-invariance of Absent Siphons in Mass-action Systems
This section gives an alternate proof of the following result used in Section 2.4, originally due to Angeli, De Leenheer, and Sontag [5].
Lemma 2.19 ([5], Proposition 5.5).
Fix any assignment of positive mass-action rate constants. Let \(\Omega \subseteq \Lambda\) be a set of species. Then \(\Omega\) is a siphon if and only if, for any state \(\mathbf {c}\) such that \(\Omega \cap [\mathbf {c}] = \emptyset\) and any state \(\mathbf {d}\) that is mass-action reachable from \(\mathbf {c}\), \(\Omega \cap [\mathbf {d}] = \emptyset\).
Proof.
Let \(\mathcal {C}= (\Lambda ,R)\) be the CRN, with positive mass-action rate constants assigned.
To see the forward direction, let \(\Omega\) be a siphon, and let \(\mathbf {c}\) be a state such that \(\Omega \cap [\mathbf {c}] = \emptyset\). Consider the reduced CRN \(\mathcal {C}_\mathrm{red}= (\Lambda _\mathrm{red}, R_\mathrm{red})\) where we remove all species in \(\Omega\) (i.e., \(\Lambda _\mathrm{red}= \Lambda \setminus \Omega\)) and all reactions referencing them (i.e., \(R_\mathrm{red}= \lbrace \langle \mathbf {r},\mathbf {p}\rangle \in R \mid [\mathbf {r}] \cap \Omega = \emptyset \text{ and } [\mathbf {p}] \cap \Omega = \emptyset \rbrace\)). Let \(\mathbf {c}_\mathrm{red}= \mathbf {c}\upharpoonright \Lambda _\mathrm{red}\). Let \({\boldsymbol \rho }^{\prime }:\mathbb {R}_{\ge 0}\rightarrow \mathbb {R}_{\ge 0}^{\Lambda _\mathrm{red}}\) be the mass-action trajectory of \(\mathcal {C}_\mathrm{red}\) starting at \(\mathbf {c}_\mathrm{red}\). Define the trajectory \({\boldsymbol \rho }: \mathbb {R}_{\ge 0}\rightarrow \mathbb {R}_{\ge 0}^\Lambda\) of \(\mathcal {C}\) by \({\boldsymbol \rho }_S(t) = {\boldsymbol \rho }^{\prime }_S(t)\) if \(S \in \Lambda _\mathrm{red}\) and \({\boldsymbol \rho }_S(t) = 0\) otherwise, i.e., \({\boldsymbol \rho }\) keeps all of \(\Omega\) at 0 and otherwise follows \({\boldsymbol \rho }^{\prime }\).
We claim that \({\boldsymbol \rho }\) is a solution to the mass-action ODEs of \(\mathcal {C}\). Since all polynomials, such as those defining mass-action ODEs, are locally Lipschitz (have bounded derivatives in some open set around every point), the Picard-Lindelöf Theorem implies that the mass-action ODEs have a unique solution. Hence, \({\boldsymbol \rho }\) is the only solution to the mass-action ODEs of \(\mathcal {C}\) starting at \(\mathbf {c}\). Since \({\boldsymbol \rho }_S(t) = 0\) for all \(S \in \Omega\) and \(t \ge 0\), this implies that any \(\mathbf {d}\) mass-action reachable from \(\mathbf {c}\) obeys \(\Omega \cap [\mathbf {d}] = \emptyset\).
To show that \({\boldsymbol \rho }\) is a valid solution to the mass-action ODEs of \(\mathcal {C}\), we need to check that
for all times \(t \ge 0\), where \(\mathbf {A}({\boldsymbol \rho }(t))\) is the vector of reaction rates at time \(t \ge 0\) (Section 2.4) and \(\mathbf {M}\) is the stoichiometry matrix (Section 2.2) converting reaction rates to species derivatives. First, let us show that \(\mathbf {A}_\alpha ({\boldsymbol \rho }(t)) = 0\) for all \(t \ge 0\) and \(\alpha \not\in R_\mathrm{red}\). If \(\alpha \notin R_\mathrm{red}\), then there is some species in \(\Omega\) that is either a reactant or product of \(\alpha\), and by the fact that \(\Omega\) is a siphon, we know that there must necessarily be a reactant S of \(\alpha\) that is in \(\Omega\). Because \({\boldsymbol \rho }_S(t) = 0\) for all \(t \ge 0\) by definition and because for mass action ODEs \(\mathbf {A}_\alpha ({\boldsymbol \rho }(t)) = 0\) if \({\boldsymbol \rho }_S(t) = 0\) for any reactant S of \(\alpha\), this shows that \(\mathbf {A}_\alpha ({\boldsymbol \rho }(t)) = 0\) for all \(t \ge 0\) and \(\alpha \not\in R_\mathrm{red}\). Thus, the right-hand side of Equation (A.1) consists entirely of the contributions of reactions in \(\mathcal {C}_\mathrm{red}\),
From this and the fact that \({\boldsymbol \rho }^{\prime }\) is a solution of the mass-action ODEs of \(\mathcal {C}_\mathrm{red}\) it follows that \({\boldsymbol \rho }\) is a solution of the mass-action ODEs for \(\mathcal {C}\).
To see the reverse direction, let \(\alpha = \langle \mathbf {r},\mathbf {p}\rangle\) be a reaction with a product \(S \in \Omega\); it suffices to show that \(\alpha\) has a reactant in \(\Omega\). If S itself is a reactant in \(\alpha\), then we are done, so assume otherwise; then \(\alpha\) produces S.
Let \(\mathbf {c}\) be a state with \([\mathbf {c}] = \Lambda _\mathrm{red}\), i.e., exactly species not in \(\Omega\) are present. In particular, \(\mathbf {c}(S) = 0\). We claim that \(\mathbf {A}_\alpha (\mathbf {c}) = 0\), i.e., \(\alpha\) has rate 0 in \(\mathbf {c}\). To see why, for the sake of contradiction, suppose \(\mathbf {A}_\alpha (\mathbf {c}) \gt 0\). To have \(\mathbf {d}(S) = 0\) for all \(\mathbf {d}\) mass-action reachable from \(\mathbf {c}\), S as a function of time is the constant 0, so \(dS/dt = 0\) in \(\mathbf {c}\). Then, to balance \(\alpha\)’s production of S in \(\mathbf {c}\) to maintain \(dS/dt = 0\), there must be some other reaction \(\beta\) with \(\mathbf {M}(S, \beta) \lt 0\) (so \(\beta\) consumes S) and \(\mathbf {A}_\beta (\mathbf {c}) \gt 0\). Since \(\beta\) consumes S, we know that S is a reactant in \(\beta\). But, since \(\mathbf {c}(S) = 0\), the rate of any reaction consuming S is 0 in \(\mathbf {c}\), a contradiction. Thus, \(\mathbf {A}_\alpha (\mathbf {c}) = 0\).
Since \([\mathbf {c}] = \Lambda _\mathrm{red}\) (all species outside of \(\Omega\) are present), to have mass-action rate 0 in \(\mathbf {c}\), \(\alpha\) must have at least one reactant in \(\Omega\) (recall reaction rate constants in a mass-action system are strictly positive). So, \(\Omega\) is a siphon.□
The definition of a valid rate schedule (Definition 2.22, part (2)) requires that if some reactant is 0, then the reaction rate is 0. However, the converse implication (if a reaction rate is 0, then some reactant must be 0) holds for mass-action but not more general rate schedules such as segment-reachability, which are allowed to “starve” applicable reactions by holding their rates at 0. The reverse direction of the proof of Lemma 2.19 uses this converse implication, but the forward direction uses only the more general implication of Definition 2.22, part (2). For the forward direction, the key property used from mass-action is its determinism: It has unique solutions, so to show that the siphon remains absent in all possible reachable states it suffices to show that there is just one solution where the siphon remains absent.
B Max-min Representation of Continuous Piecewise Linear Functions
Here, we prove a slight generalization of Ovchinnikov’s theorem [42]. In Ovchinnikov’s original paper, he only considers piecewise affine functions (in Ovchinnikov’s terminology, piecewise “linear” functions) that are defined on closed domains (that is, closures of open subsets of \(\mathbb {R}^n\)). However, the key proof techniques of Reference [42] did not crucially use this fact. In fact, we apply theorem 5.11 on non-closed domains such as the sets \(D_U\) in the proof of Lemma 5.16. For completeness, we prove the variant of the theorem not requiring D to be closed.
Theorem 5.11 ([42], Theorem 2.1).
Let \(D \subseteq \mathbb {R}^k\) be convex. For every continuous piecewise affine function \(f:D \rightarrow \mathbb {R}\) with components \(g_1,\ldots ,g_p\), there exists a family \(S_1,\ldots ,S_q \subseteq \lbrace 1,\ldots ,p\rbrace\) such that, for all \(\mathbf {x}\in D\), \(f(\mathbf {x}) = \max \nolimits _{i \in \lbrace 1,\ldots ,q\rbrace } \min \nolimits _{j \in S_i} g_j(\mathbf {x}).\)
To prove the theorem, we first prove three lemmas. The first technical lemma is implicit in Reference [42]. The second and third lemmas correspond to Lemmas 2.1 and 2.2 of Reference [42]. The proofs we give of the second and third lemmas are almost identical in content to the proofs of the corresponding lemmas in Reference [42], except for the fact that we consider piecewise affine functions defined over more general subsets of \(\mathbb {R}^n\). The same is true for our proof of Theorem 5.11, which is again almost identical to the proof of Theorem 2.1 in Reference [42].
Lemma B.1.
If \(f: [a,b] \rightarrow \mathbb {R}\) is a continuous piecewise affine function with components \(\left\lbrace g_1, \ldots , g_n\right\rbrace\), then there are finitely many numbers \(a = x_0 \lt x_1 \lt \ldots \lt x_m = b\) such that f is affine on \([x_k, x_{k + 1}]\) for all k.
Proof.
Without loss of generality, we can assume that all of the component functions \(g_i\) are distinct affine functions. For each i between 1 and n, let \(D_i\) be the set of \(x \in [a,b]\) such that \(g_i(x) = f(x)\). For each i, both \(g_i\) and f are continuous, so \(D_i\) is closed. Let S be the subset of \([a,b]\) consisting of points \(x \in [a,b]\) where x is a member of more than one \(D_i\). For each pair \(i \ne j\), we know that \(g_i\) and \(g_j\) are distinct affine functions, so there can be at most one \(x \in [a,b]\) such that \(g_i(x) = g_j(x)\). This implies that \(D_i\) and \(D_j\) can intersect in at most one point, so S must be a finite set. Let \(x_0\cdots x_m\) be the elements of \(S \cup \left\lbrace a,b\right\rbrace\).
Now for a given k write I for the interval \((x_k, x_{k + 1})\) and consider the restriction of f to I. Pick a random point \(c \in I\) and suppose \(f(c) = g_l(c)\). Then clearly \(D_l \cap I\) is nonempty. Because \(D_l\) is closed in \([a,b]\), by definition \(D_l \cap I\) is closed relative to I. Because \(S \cap I = \emptyset\), we also know that
\begin{equation*} D_l \cap I = I \setminus \left(\bigcup _{i \ne l} D_i \right), \end{equation*}
so \(D_l \cap I\) is open relative to I. But the only subset of an interval that is both open and closed is the whole interval, so \(I \subseteq D_l\). Therefore, \(f = g_l\) when restricted to I, and by continuity, we see that f is affine on the closure of I as well.□
Note that we define piecewise affine functions to have only finitely many components—without this assumption, the above lemma is false.
Lemma B.2.
Let \(f: [a,b] \rightarrow \mathbb {R}\) be a continuous piecewise affine function. Let \(\left\lbrace g_1, \ldots , g_n\right\rbrace\) be its set of components. Then there is some k such that
We will first prove the result for \(f(a) = f(b) = 0\) and then show how this implies the general case. Given this assumption, if one of the \(g_i\) is the zero function, then we are done. If not, since all of the component functions \(g_i\) are affine, each \(g_i\) can have at most one zero. Since f has finitely many components, this implies that f has finitely many zeros. Let c be the smallest number such that \(c \gt a\) and \(f(c) = 0\).
By Lemma B.1, we know that there are some x and y with \(a \lt x \lt y \lt c\) and component functions \(g_k\) and \(g_l\) such that \(f = g_k\) on \([a, x]\) and \(f = g_l\) on \([y, c]\). If the slope of either \(g_k\) or \(g_l\) is non-negative, then we are done:
and similarly for \(g_l\). But \(g_k\) and \(g_l\) cannot both have negative slope, for then \(f(x) = g_k(x) \lt 0\) and \(f(y) = g_l(y) \gt 0\), so by the intermediate value theorem there would be some z between x and y such that \(f(z) = 0\). This contradicts our assumption that c was the smallest number with \(c \gt a\) and \(f(c) = 0\). This concludes the proof assuming that \(f(a)=f(b)=0.\)
To deduce the result for a general continuous piecewise affine function from this special case, subtract the affine function
Let D be a convex subset of \(\mathbb {R}^n\) and let \(f: D \rightarrow \mathbb {R}\) be a continuous piecewise affine function. If the components of f are \(\left\lbrace g_1 \ldots g_n\right\rbrace\), then for every pair of vectors \(\mathbf {a}\) and \(\mathbf {b}\) in D, there is some k such that
Because D is convex, the straight-line interval between \(\mathbf {a}\) and \(\mathbf {b}\) is contained in D. Apply Lemma B.2 to the restriction of f to this interval.□
Because there is always some component function \(g_j\) with \(g_j(\mathbf {b}) = f(\mathbf {b})\), we see that \(F_\mathbf {b}(\mathbf {b}) = f(\mathbf {b})\) for every \(\mathbf {b}\in D\). Also, by Lemma B.3, we know that for every \(\mathbf {a}\in D\), there is some component function \(g_k \in S_\mathbf {b}\) with \(g_k(\mathbf {a}) \le f(\mathbf {a})\), so \(F_\mathbf {b}(\mathbf {a}) \le f(\mathbf {a})\) for every pair \(\mathbf {b}, \mathbf {a}\in D\). This implies that
Since \(\left\lbrace 1 \ldots p\right\rbrace\) is a finite set, it has only finitely many subsets, so each \(S_\mathbf {b}\) is equal to one of finitely many sets \(S_j\). We can therefore replace the maximum over all \(\mathbf {b}\in D\) in Equation (B.1) with a maximum over finitely many functions.□
C Finding Rational Solutions to Systems of Linear Equations
It is well-known that a system of linear equations with rational coefficients has a rational solution if and only if it has a real solution. The following result shows the slightly generalized claim that rational solutions exist arbitrarily close to all real solutions (i.e., the rational solutions are dense in the real solutions).
Lemma C.1.
Let \(A\mathbf {x} = \mathbf {b}\) be a system of linear equations, where A is a matrix with rational coefficients and \(\mathbf {b}\) is a vector with rational coefficients. If the equation has a solution \(\mathbf {x}\) with real coefficients, then for any \(\varepsilon \gt 0\), it has a solution with \(\mathbf {x}^{\prime }\) rational coefficients such that \(||\mathbf {x}^{\prime } - \mathbf {x}|| \lt \varepsilon\).
Proof.
Let n be the number of rows of A and the length of \(\mathbf {b}\). Let m be the number of columns of A and the length of \(\mathbf {x}\). Because A has rational entries, using elementary row and column operations it can be decomposed as PNQ where P is an \(m \times m\) invertible rational matrix, Q is an \(n \times n\) invertible rational matrix, and
\begin{equation*} N_{ij} = {\left\lbrace \begin{array}{ll}1 & i = j \text{ and } i \le r \\ 0 & \text{otherwise,} \end{array}\right.} \end{equation*}
where r is the rank of M. Let \(\mathbf {y} = Q\mathbf {x}\) and let \(\mathbf {c} = P^{-1} \mathbf {b}\), so \(N\mathbf {y} = \mathbf {c}\). If \(r = n\), then all of the entries of \(\mathbf {y}\) must be rational, since \(\mathbf {y}_i = \mathbf {c}_i\) for all i and all of the entries of \(\mathbf {c}_i\) are rational. Then the entries of \(\mathbf {x}\) must also all be rational, since \(\mathbf {x} = Q^{-1}\mathbf {y}\) and Q has all rational entries. As a result, if \(r = n\), then we can just take \(\mathbf {x}^{\prime } = \mathbf {x}\).
is the operator norm of \(||Q^{-1}||\). Now let \(\mathbf {y}^{\prime }\) be a vector such that \(\mathbf {y}^{\prime }_i = \mathbf {y}_i\) for \(i \le r\) and \(\mathbf {y}^{\prime }_i\) is a rational number such that \(|\mathbf {y}^{\prime }_i - \mathbf {y}_i| \lt \delta\) for \(i \gt r\). All of the components of \(\mathbf {y}^{\prime }\) are rational: \(\mathbf {y}^{\prime }_i\) is rational by construction for \(i \gt r\), and \(\mathbf {y}^{\prime }_i = \mathbf {y}_i = \mathbf {c}_i\) is rational for \(i \le r\). Moreover, the fact that \(\mathbf {y}^{\prime }_i = \mathbf {c}_i\) for \(i \le r\) shows that \(N\mathbf {y}^{\prime } = \mathbf {c}\).
If we take \(\mathbf {x}^{\prime } = Q^{-1}\mathbf {y}^{\prime }\), then all of the components of \(\mathbf {x}^{\prime }\) are rational, and \(M\mathbf {x}^{\prime } = \mathbf {b}\), since
This shows that \(\mathbf {x}^{\prime }\) is our desired solution.□
D Bounding Reaction Fluxes in Straight-line Reachability
This section is devoted to proving Lemma D.4. Intuitively, it shows that if a CRN can reach from state \(\mathbf {c}\) to state \(\mathbf {d}\) by a straight line (of length \(\Vert \mathbf {d}- \mathbf {c}\Vert\)), then the reaction fluxes required can be bounded by \(O(\Vert \mathbf {d}- \mathbf {c}\Vert)\). This is nontrivial, since one can have reactions that cancel, e.g., \(X \rightarrow Y\) and \(Y \rightarrow X\). The same straight line from \(\mathbf {c}\) to \(\mathbf {d}\) could result from arbitrarily large but equal fluxes of each reaction (plus some other reactions). Lemma D.4 states that we never need arbitrarily large fluxes to get from \(\mathbf {c}\) to \(\mathbf {d}\).
Definition D.1.
A convex polyhedral coneC is a subset of a vector space V such that there exist vectors \(v_1, \ldots , v_k \in V\) so C is exactly the set of \(x \in V\) that can be written as \(x = \sum _i \lambda _iv_i\) with all \(\lambda _i \ge 0\). Such a set \(S = \left\lbrace v_1 \ldots v_k\right\rbrace\) is called a spanning set for \(C,\) and we say that C is spanned byS. Given a set \(S = \left\lbrace v_1 \ldots v_k\right\rbrace\) of vectors in V, we write \(C_S\) for the convex polyhedral cone spanned by S.
Lemma D.2.
Let C be a convex cone with spanning set \(S = \left\lbrace v_1 \ldots v_k\right\rbrace\). Let \(\mathcal {I}\) be the collection of all linearly independent subsets of S. Then
Suppose to the contrary that there was some \(x \in C\) not contained in \(C_{S^{\prime }}\) for any \(S^{\prime } \in \mathcal {I}\). Let T be a minimal subset of S such that \(x \in C_T\). Then, since \(x \in C_T\), we can write x as
with some \(\lambda ^{\prime }_i = 0\). This will imply that \(x \in C_{T \setminus \left\lbrace v_i\right\rbrace }\), contradicting the minimality of T. By assumption, since \(x \in T\), we know that \(T \notin \mathcal {I}\), so there is some nontrivial linear relationship
Let us show that \(\lambda _j + C\mu _j \ge 0\) for all j. Because \(\lambda _i \gt 0\) for all i, we know that \(C \gt 0\). As a result, if \(\mu _j \ge 0,\) then necessarily \(\lambda _j + C\mu _j \ge 0\). However, if \(\mu _j \lt 0\), then \(C \le -\lambda _j/\mu _j\), so
where \(\lambda ^{\prime }_{i_0} = 0\). Thus, \(x \in C_{T\setminus v_{i_0}}\), contradicting the minimality of T and therefore the existence of x.□
Lemma D.3.
Let C be a convex cone with spanning set \(S = \left\lbrace v_1\ldots v_k\right\rbrace\). Then there is a constant K depending only on S so for any \(x \in C\), there is some representation of x as
\begin{equation*} x = \sum _i \lambda _i v_i \end{equation*}
with \(0 \le \lambda _i \le K||x||\). Additionally, if there is some subset \(T \subseteq S\) such that \(x \in C_T\), then the above representation can be chosen with \(\lambda _i = 0\) for any \(v_i \notin T\).
Proof.
Let \(\mathcal {I}\) be the collection of linearly independent subsets of S. If \(S^{\prime } \in \mathcal {I}\) and \(x \in C_{S^{\prime }}\) then because the vectors \(v_1 \ldots v_k\) in \(S^{\prime }\) are linearly independent, there is a unique way to write
and by Reference [33, Lemma 2.4-1], we know there is some constant \(K_{S^{\prime }}\) so \(0 \le \lambda _i \le K_{S^{\prime }} \Vert x\Vert\). Take \(K = \max _{S^{\prime } \in \mathcal {I}} K_{S^{\prime }}\). Applying Lemma D.2 to \(C_T\), we can find some linearly independent collection of vectors \(T^{\prime }\) such that \(T^{\prime } \subseteq T\) and \(x \in C_{T^{\prime }}\). Then
\begin{equation*} x = \sum _{v_i \in T^{\prime }}\lambda _iv_i \end{equation*}
with \(0 \le \lambda _i \le K_{T^{\prime }}\Vert x\Vert\). Because \(T^{\prime } \subseteq T \subseteq S\), we know that \(T^{\prime } \in \mathcal {I}\), so \(K_{T^{\prime }} \le K\), and, since \(T^{\prime } \subseteq T\), we know that above sum only ranges over vectors \(v_i \in T\).□
Finally, we have the main result of this section.
Lemma D.4.
Fix a CRN \(\mathcal {C}\) and suppose a flux vector \(\mathbf {u}\in \mathbb {R}_{\ge 0}^R\) is applicable at a state \(\mathbf {c}\in \mathbb {R}_{\ge 0}^\Lambda\). Let \(\ell = ||M\mathbf {u}||\) be the length of the straight-line segment in \(\mathbb {R}^\Lambda\) given by \(M\mathbf {u}\). Then there exists a constant K, depending only on \(\mathcal {C}\) and independent of \(\mathbf {u}\) and \(\mathbf {c}\), such that there exists a flux vector \(\mathbf {u}^{\prime }\) where \(\mathbf {u}^{\prime }\) is also applicable at \(\mathbf {c}\), the length of \(\mathbf {u}^{\prime }\) is bounded as \(||\mathbf {u}^{\prime }|| \le K\ell\), and \(M\mathbf {u}= M\mathbf {u}^{\prime }\).
Proof.
Let \(\left\lbrace v_1\ldots v_m\right\rbrace\) be the standard basis vectors of \(\mathbb {R}^R\). Because \(\mathbf {u}= \sum _i \lambda _i v_i\) with all \(\lambda _i \ge 0\), we know that \(M\mathbf {u}\) is contained in the convex polyhedral cone spanned by \(\left\lbrace Mv_1 \ldots Mv_m\right\rbrace\). By Lemma D.3, we know that there are \(\lambda ^{\prime }_i\) such that \(M\mathbf {u}= \sum _i \lambda ^{\prime }_i Mv_i\) and \(0 \le \lambda ^{\prime }_i \le K||M\mathbf {u}||\), and moreover \(\lambda _i^{\prime } = 0\) whenever \(\lambda _i = 0\). Let \(\mathbf {u}^{\prime } = \sum _i \lambda ^{\prime }_iv_i\). Then \(\mathbf {u}^{\prime }\) is still applicable at \(\mathbf {c}\) because \(\lambda ^{\prime }_i \gt 0\) implies \(\lambda _i \gt 0\) and we assumed that \(\mathbf {u}\) was applicable at \(\mathbf {c}\). Also, note that
so \(\mathbf {u}^{\prime }\) is our desired vector in \(\mathbb {R}_{\ge 0}^R\).□
E Partial States and Reachability
In this section, we define a notion of “partial” states and reachability, which are used in the proof of theorem 3.3. If \(\Delta \subsetneq \Lambda\), then we say \(\mathbf {p}\in \mathbb {R}_{\ge 0}^\Delta\) is a partial state. Given a state \(\mathbf {c}\in \mathbb {R}_{\ge 0}^\Lambda\), recall that \(\mathbf {c}\upharpoonright \Delta\) is \(\mathbf {c}\) restricted to \(\Delta\), i.e., the partial state \(\mathbf {p}= \mathbf {c}\upharpoonright \Delta \in \mathbb {R}_{\ge 0}^\Delta\) such that \(\mathbf {p}(S) = \mathbf {c}(S)\) for all \(S \in \Delta\).
Let \(k \in \mathbb {N}\cup \lbrace \infty \rbrace .\) Given a state \(\mathbf {c}\in \mathbb {R}_{\ge 0}^\Lambda\) and a partial state \(\mathbf {p}\in \mathbb {R}_{\ge 0}^\Delta\), we write \(\mathbf {c}\rightsquigarrow ^k \mathbf {p}\) if there is a sequence of states \(\mathbf {b}_0, \dots , \mathbf {b}_{k} \in \mathbb {R}_{\ge 0}^\Lambda\) such that \(\mathbf {c}= \mathbf {b}_0 \rightarrow ^1\mathbf {b}_1 \rightarrow ^1\mathbf {b}_2 \rightarrow ^1\cdots \rightarrow ^1\mathbf {b}_{k}\), with \(\mathbf {p}= \mathbf {b}_k \upharpoonright \Delta\) if \(k \in \mathbb {N}\), or \(\mathbf {p}= \lim \nolimits _{i \rightarrow \infty } (\mathbf {b}_i \upharpoonright \Delta)\) if \(k = \infty\). We write \(\mathbf {c}\rightsquigarrow \mathbf {p}\) via \(\left(\mathbf {b}_i\right)_{i=1}^k\) if \(\mathbf {c}\rightsquigarrow ^k \mathbf {p}\) for some \(k \in \mathbb {N}\cup \lbrace \infty \rbrace\) with intermediate states \(\mathbf {b}_0,\mathbf {b}_1,\dots\) as above, or simply \(\mathbf {c}\rightsquigarrow \mathbf {p}\) when the intermediate states \(\mathbf {b}_0,\mathbf {b}_1,\dots\) are implicit. We write \(\mathbf {c}\rightsquigarrow ^\infty _{\mathrm{ss}}\mathbf {p}\) via \(\left(\mathbf {b}_i\right)_{i=1}^k\) if \(\mathbf {c}= \mathbf {b}^{\prime }_0 \rightarrow ^1 \mathbf {b}^{\prime }_1 \rightarrow ^1 \dots\), where for some subsequence \(\mathbf {b}_1,\mathbf {b}_1,\dots\) of \(\mathbf {b}^{\prime }_1,\mathbf {b}^{\prime }_1,\dots\), we have \(\mathbf {p}= \lim \nolimits _{i\rightarrow \infty } (\mathbf {b}_i \upharpoonright \Delta)\), i.e., an infinite subsequence of states converges on concentrations in \(\Delta\).
Note that if there is a state \(\mathbf {d}\) such that \(\mathbf {c}\rightsquigarrow \mathbf {d}\) and \(\mathbf {d}\upharpoonright \Delta = \mathbf {p}\), then \(\mathbf {c}\rightsquigarrow \mathbf {p}\), but it is not apparent from the definition that this is the only way for a partial state to be reachable if the number of line segments is infinite. In particular, it could be that the sequence \(\mathbf {b}_0, \mathbf {b}_1, \ldots\) does not converge to any state (even though the partial states \(\mathbf {b}_0 \upharpoonright \Delta , \mathbf {b}_1 \upharpoonright \Delta , \ldots\) converge to \(\mathbf {p}\)), if some concentration values outside of \(\Delta\) do not converge (for instance, they may oscillate or go to infinity).
Our goal now is to show that in fact, if a partial state \(\mathbf {p}\) is reachable (or even merely \(\rightsquigarrow ^\infty _{\mathrm{ss}}\) reachable), then there is a particular state reachable whose restriction to \(\Delta\) is \(\mathbf {p}\).
Theorem E.1.
Let \(\Delta \subsetneq \Lambda\), let \(\mathbf {c}\in \mathbb {R}_{\ge 0}^\Lambda\) be a state, and \(\mathbf {p}\in \mathbb {R}_{\ge 0}^\Delta\) be a partial state. If \(\mathbf {c}\rightsquigarrow ^\infty _{\mathrm{ss}}\mathbf {p}\), then there is a state \(\mathbf {d}\) such that \(\mathbf {c}\rightsquigarrow \mathbf {d}\) and \(\mathbf {p}=\mathbf {d}\upharpoonright \Delta\).
Furthermore, if \(\mathbf {c}\rightsquigarrow ^\infty _{\mathrm{ss}}\mathbf {p}\) via \(\left(\mathbf {b}_i\right)_{i=1}^k\) for \(k \in \mathbb {N}\cup \lbrace \infty \rbrace\), then there is a partition of \(\Lambda\) into \(\Lambda _\mathrm{b}\) and \(\Lambda _\mathrm{u}\), with \(\Delta \subseteq \Lambda _\mathrm{b}\), and a subsequence \(\left(\mathbf {r}_i\right)_i\) of \(\left(\mathbf {b}_i\right)_i\) so \(\lim _{i\rightarrow \infty } \mathbf {r}_i(S) = \infty\) for all \(S \in \Lambda _\mathrm{u}\), and \(\lim _{i\rightarrow \infty }\mathbf {r}_i(S) = \mathbf {d}(S)\) for all \(S \in \Lambda _\mathrm{b}\).
Proof.
The finite case is immediate from the definition of \(\rightsquigarrow\) for partial states (choose \(\mathbf {d}= \mathbf {b}_k\) and \(\Lambda _\mathrm{u}= \emptyset\)), so assume \(\mathbf {c}\rightsquigarrow ^\infty _{\mathrm{ss}}\mathbf {p}\). Then there is an infinite sequence of states (the converging subsequence in the definition of \(\rightsquigarrow ^\infty _{\mathrm{ss}}\)) \(\mathbf {b}_0, \mathbf {b}_1, \dots \in \mathbb {R}_{\ge 0}^\Lambda\) such that \(\mathbf {c}= \mathbf {b}_0\), each \(\mathbf {b}_i \rightsquigarrow ^{k_i} \mathbf {b}_{i+1}\) for some \(k_i \in \mathbb {N}\), and \(\mathbf {p}= \lim _{i \rightarrow \infty } \mathbf {b}_i \upharpoonright \Delta\).
Let \(\Gamma = \Lambda \setminus \Delta\) be the species outside of \(\Delta\), which may not be converging in the subsequence \((\mathbf {b}_i)_i\). Intuitively, our goal will be to make some species in \(\Gamma\) converge and all others simultaneously diverge to infinity. Then, by driving the concentrations of the diverging species sufficiently large, then removing them from the system, we obtain a reduced CRN where all species converge to a single state. We apply Corollary 2.15 to this reduced CRN to find a finite path of \(m+1\) segments reaching this state. Finally, we argue that this path is applicable even in the original CRN, because the species that were removed first had their concentrations driven large enough that the finite path does not have sufficient flux to send any of them to 0. Thus, at the end of this finite path, all species converge to some concentration.
Formally, define the set \(\Gamma _\mathrm{u}\) (the “unbounded” species in \(\Gamma\)) iteratively as follows: If there is any \(S_1 \in \Gamma\) such that \(\limsup _{i\rightarrow \infty } \mathbf {b}_i(S_1) = \infty\), then put \(S_1\) in \(\Gamma _\mathrm{u}\), otherwise \(\Gamma _\mathrm{u}\) is defined to be \(\emptyset\). Then pick a subsequence \(\mathbf {b}^{\prime }_0, \mathbf {b}^{\prime }_1,\ldots\) of \(\mathbf {b}_0,\mathbf {b}_1,\ldots\) such that, for all \(i\in \mathbb {N}\), \(\mathbf {b}^{\prime }_i(S_1) \ge i\); such a subsequence exists, since \(\limsup _{i\rightarrow \infty } \mathbf {b}_i(S_1) = \infty\). Now, from that subsequence, if any species \(S_2\) obeys \(\limsup _{i\rightarrow \infty } \mathbf {b}^{\prime }_i(S_2) = \infty\), then place \(S_2\) in \(\Gamma _\mathrm{u}\) and choose a subsequence \(\mathbf {b}^{\prime \prime }_0,\mathbf {b}^{\prime \prime }_1,\ldots\) of \(\mathbf {b}^{\prime }_0,\mathbf {b}^{\prime }_1,\ldots\) where \(\mathbf {b}^{\prime \prime }_i(S_2) \ge i\) for each \(i \in \mathbb {N}\); note that because \(\mathbf {b}^{\prime \prime }_i = \mathbf {b}^{\prime }_{i^{\prime }}\) for \(i^{\prime } \ge i\), we also have \(\mathbf {b}^{\prime \prime }_i(S_1) \ge i\).
Repeat this alternation of choosing a species \(S_j\) to put in \(\Gamma _\mathrm{u}\) and picking a subsequence, until in the final subsequence \(\mathbf {b}^{\prime \prime \prime }_0,\mathbf {b}^{\prime \prime \prime }_1,\ldots\) every remaining species S obeys \(\limsup _{i\rightarrow \infty } \mathbf {b}^{\prime \prime \prime }_i(S) \lt \infty\). Define the set \(\Gamma _\mathrm{b}= \Gamma \setminus \Gamma _\mathrm{u}\); by construction each species in \(\Gamma _\mathrm{b}\) has bounded concentrations in \(\mathbf {b}^{\prime \prime \prime }_0,\mathbf {b}^{\prime \prime \prime }_1,\dots\) (Though some may have had unbounded concentrations in the original sequence \(\mathbf {b}_0,\mathbf {b}_1,\dots\))
Now, since concentrations of species in \(\Gamma _\mathrm{b}\) are in a closed, bounded (i.e., compact) set, there is a subsequence \(\mathbf {r}_0,\mathbf {r}_1,\dots\) of \(\mathbf {b}^{\prime \prime \prime }_0,\mathbf {b}^{\prime \prime \prime }_1,\dots\) that converges on concentrations for species in \(\Gamma _\mathrm{b}\). Furthermore, all subsequences of \(\mathbf {b}_0,\mathbf {b}_1,\dots\) we have taken so far can be assumed without loss of generality to contain \(\mathbf {b}_0\). Thus, the subsequence \(\mathbf {r}_0,\mathbf {r}_1,\dots\) obeys
for all \(S \in \Gamma _\mathrm{u}\) and \(i\in \mathbb {N}\), \(\mathbf {r}_i(S) \ge i\) (concentrations in \(\Gamma _\mathrm{u}\) increase to infinity simultaneously), and
(3)
for all \(S \in \Gamma _\mathrm{b}\cup \Delta\), \(\lim _{i \rightarrow \infty } \mathbf {r}_i(S)\) exists and is finite (all other concentrations converge).
Consider the “quotient” CRN \(\mathcal {C}_\mathrm{q}= (\Lambda _\mathrm{q}, R_\mathrm{q})\), where \(\Lambda _\mathrm{q}= \Lambda \setminus \Gamma _\mathrm{u}= \Delta \cup \Gamma _\mathrm{b}\), and \(R_\mathrm{q}\) is defined by taking each reaction from R and removing any species from it in \(\Gamma _\mathrm{u}\). For example, if \(A,B \in \Gamma _\mathrm{u}\) and \(C,D,E \not\in \Gamma _\mathrm{u}\), then the reaction \(A + B + 2C \rightarrow 2A + D + E\) becomes \(2C \rightarrow D+E\). Then for each state \(\mathbf {r}\in \mathbb {R}_{\ge 0}^\Lambda\) of the original CRN \(\mathcal {C}\), its partial state \(\mathbf {r}\upharpoonright \Lambda _\mathrm{q}\) is a (normal) state of \(\mathcal {C}_\mathrm{q}.\)
For each \(i\in \mathbb {N}\), let \(\mathbf {q}_i = \mathbf {r}_i \upharpoonright \Lambda _\mathrm{q}\) be the state of \(\mathcal {C}_\mathrm{q}\) corresponding to \(\mathbf {r}_i\). Note that \(\mathbf {q}_0 = \mathbf {c}\upharpoonright \Lambda _\mathrm{q}\). By the definition of \(\Lambda _\mathrm{q}= \Delta \cup \Gamma _\mathrm{b}\) and the convergence of concentrations in both \(\Delta\) (by the definition of \(\mathbf {c}\rightsquigarrow \mathbf {p}\)) and \(\Gamma _\mathrm{b}\) (shown as part (3) above), the sequence \(\mathbf {q}_0, \mathbf {q}_1,\dots\) converges to some state \(\mathbf {q}\in \mathbb {R}_{\ge 0}^{\Lambda _\mathrm{q}}\). This implies that \(\mathbf {q}_i \rightsquigarrow ^\infty _{\mathrm{ss}}\mathbf {q}\) for each \(i \in \mathbb {N}\). Note that \(\mathbf {q}\upharpoonright \Delta = \mathbf {p}\).
Since the \(\mathbf {q}_i\)’s converge to \(\mathbf {q}\), there is some \(i_0\) such that, for all \(i \ge i_0\), \(\Vert \mathbf {q}- \mathbf {q}_i \Vert \le 1\). Since each \(\mathbf {q}_i \rightsquigarrow ^\infty _{\mathrm{ss}}\mathbf {q}\), by Theorem 2.14, \(\mathbf {q}_i \rightsquigarrow ^{m+1} \mathbf {q}\), where \(m = \min \lbrace |\Lambda _\mathrm{q}|, |R_\mathrm{q}| \rbrace\), and for some constant K depending only on \(\mathcal {C}_\mathrm{q}\), the total reaction fluxes do not exceed K along the entire path of \(m+1\) segments.
Let \(C = \max \lbrace \mathbf {r}^{\prime }(S) \mid S \in \Gamma _\mathrm{u}, \langle \mathbf {r}^{\prime },\mathbf {p}^{\prime } \rangle \in R \rbrace\) be the maximum reactant coefficient of any reactant in \(\Gamma _\mathrm{u}\) for any reaction in R. Choose \(i_1 = \max \lbrace C \cdot K, i_0 \rbrace\). Then \(\mathbf {q}_{i_1} \rightsquigarrow \mathbf {q}\) with total reaction flux at most K, since \(i_1 \ge i_0\).
Consider running the same reaction fluxes on the original CRN, from the state \(\mathbf {b}_{i^{\prime }}\), choosing \(i^{\prime }\) such that \(\mathbf {b}_{i^{\prime }} \upharpoonright \Lambda _\mathrm{q}= \mathbf {q}_{i_1}\). By property (2) above, there is sufficient concentration \(C \cdot K\) of each reactant in \(\Gamma _\mathrm{u}\) in state \(\mathbf {b}_{i^{\prime }}\) for these reactions to remain applicable along the entire path. Let \(\mathbf {d}\in \mathbb {R}_{\ge 0}^\Lambda\) be the state reached at the end of this path, then \(\mathbf {b}_{i^{\prime }} \rightsquigarrow \mathbf {d}\). Note that \(\mathbf {d}\upharpoonright \Lambda _\mathrm{q}= \mathbf {q}\), and recall \(\mathbf {q}\upharpoonright \Delta = \mathbf {p}\), so \(\mathbf {d}\upharpoonright \Delta = \mathbf {p}\). Since \(\mathbf {c}\rightsquigarrow \mathbf {b}_{i^{\prime }}\), by transitivity of \(\rightsquigarrow\), we have \(\mathbf {c}\rightsquigarrow \mathbf {d}\), proving the theorem.□
Note that by Corollary 2.15, if \(\mathbf {c}\rightsquigarrow \mathbf {p}\) for a partial state \(\mathbf {p}\), then \(\mathbf {c}\rightsquigarrow ^{m+1} \mathbf {p}\), where \(m = \min \lbrace |\Lambda |,|R|\rbrace\).
References
[1]
Dan Alistarh, James Aspnes, David Eisenstat, Rati Gelashvili, and Ronald L. Rivest. 2017. Time-space trade-offs in population protocols. In 28th Annual ACM-SIAM Symposium on Discrete Algorithms. SIAM, 2560–2579.
Dan Alistarh, James Aspnes, and Rati Gelashvili. 2018. Space-optimal majority in population protocols. In 29th Annual ACM-SIAM Symposium on Discrete Algorithms. 2221–2239.
David Angeli, Patrick De Leenheer, and Eduardo D. Sontag. 2006. On the structural monotonicity of chemical reaction networks. In 45th IEEE Conference on Decision and Control. IEEE, 7–12.
David Angeli, Patrick De Leenheer, and Eduardo D. Sontag. 2007. A Petri net approach to the study of persistence in chemical reaction networks. Mathem. Biosci. 210, 2 (2007), 598–618.
Dana Angluin, James Aspnes, Zoë Diamadi, Michael Fischer, and René Peralta. 2006. Computation in networks of passively mobile finite-state sensors. Distrib. Comput. 18, 4 (2006), 235–253. DOI:
Dana Angluin, James Aspnes, and David Eisenstat. 2006. Stably computable predicates are semilinear. In PODC 2006: 25th Annual ACM Symposium on Principles of Distributed Computing. ACM Press, New York, NY, 292–299. DOI:
Dana Angluin, James Aspnes, and David Eisenstat. 2008. Fast computation by population protocols with a leader. Distrib. Comput. 21, 3 (Sept.2008), 183–199.
Amanda Belleville, David Doty, and David Soloveichik. 2017. Hardness of computing and approximating predicates and functions with leaderless population protocols. In 44th International Colloquium on Automata, Languages, and Programming (Leibniz International Proceedings in Informatics (LIPIcs)), Vol. 80. 141:1–141:14.
Olivier Bournez, Daniel S. Graça, and Amaury Pouly. 2017. Polynomial time corresponds to solutions of polynomial ordinary differential equations of polynomial length. J. ACM 64, 6 (Oct.2017). DOI:
Adam Case, Jack H. Lutz, and Donald M. Stull. 2018. Reachability problems for continuous chemical reaction networks. Natur. Comput. 17, 2 (2018), 223–230.
Ho-Lin Chen, David Doty, and David Soloveichik. 2014. Deterministic function computation with chemical reaction networks. Natur. Comput. 13, 4 (2014), 517–534. DOI:
Ho-Lin Chen, David Doty, and David Soloveichik. 2014. Rate-independent computation in continuous chemical reaction networks. In 5th Conference on Innovations in Theoretical Computer Science.
Yuan-Jyue Chen, Neil Dalchau, Niranjan Srinivas, Andrew Phillips, Luca Cardelli, David Soloveichik, and Georg Seelig. 2013. Programmable chemical controllers made from DNA. Nat. Nanotechnol. 8, 10 (2013), 755–762.
Matthew Cook, David Soloveichik, Erik Winfree, and Jehoshua Bruck. 2009. Programmability of chemical reaction networks. In Algorithmic Bioprocesses, Anne Condon, David Harel, Joost N. Kok, Arto Salomaa, and Erik Winfree (Eds.). Springer, Berlin, 543–584.
Wojciech Czerwiński and Łukasz Orlikowski. 2022. Reachability in vector addition systems is Ackermann-complete. In IEEE 62nd Annual Symposium on Foundations of Computer Science (FOCS). 1229–1240.
Élisabeth Degrand, François Fages, and Sylvain Soliman. 2020. Graphical conditions for rate independence in chemical reaction networks. In Computational Methods in Systems Biology, Alessandro Abate, Tatjana Petrov, and Verena Wolf (Eds.). Springer International Publishing, Cham, 61–78.
David Doty and Eric Severson. 2021. Ppsim: A software package for efficiently simulating and visualizing population protocols. In 19th International Conference on Computational Methods in Systems Biology. 245–253. Retrieved from https://arxiv.org/abs/2105.04702.
François Fages, Steven Gay, and Sylvain Soliman. 2015. Inferring reaction systems from ordinary differential equations. Theoret. Comput. Sci. 599 (2015), 64–78. DOI:
François Fages, Guillaume Le Guludec, Olivier Bournez, and Amaury Pouly. 2017. Strong Turing completeness of continuous chemical reaction networks and compilation of mixed analog-digital programs. In International Conference on Computational Methods in Systems Biology. Springer, 108–127.
Leszek Gąsieniec and Grzegorz Staehowiak. 2018. Fast space optimal leader election in population protocols. In 29th Annual ACM-SIAM Symposium on Discrete Algorithms. SIAM, 2653–2667.
Manoj Gopalkrishnan, Ezra Miller, and Anne Shiu. 2013. A projection argument for differential inclusions, with applications to persistence of mass-action kinetics. Symm., Integrabil. Geom.: Meth. Applic. 9, 0 (2013), 25–25.
Hooman Hashemi, Ben Chugg, and Anne Condon. 2020. Composable computation in leaderless, discrete chemical reaction networks. In 26th International Conference on DNA Computing and Molecular Programming (DNA’20). Schloss Dagstuhl-Leibniz-Zentrum für Informatik.
Shigeru Kondo and Takashi Miura. 2010. Reaction-diffusion model as a framework for understanding biological pattern formation. Science 329, 5999 (2010), 1616–1620.
James I. Lathrop, Jack H. Lutz, Robyn R. Lutz, Hugh D. Potter, and Matthew R. Riley. 2020. Population-induced phase transitions and the verification of chemical reaction networks. In 26th International Conference on DNA Computing and Molecular Programming (Leibniz International Proceedings in Informatics (LIPIcs)), Cody Geary and Matthew J. Patitz (Eds.), Vol. 174. Schloss Dagstuhl–Leibniz-Zentrum für Informatik, Dagstuhl, Germany, 5:1–5:17. DOI:
Jérôme Leroux. 2022. The reachability problem for Petri nets is not primitive recursive. In IEEE 62nd Annual Symposium on Foundations of Computer Science (FOCS). 1241–1252.
Sayed Ahmad Salehi, Keshab K. Parhi, and Marc D. Riedel. 2017. Chemical reaction networks for computing polynomials. ACS Synthet. Biol. 6, 1 (2017), 76–83.
Georg Seelig and David Soloveichik. 2009. Time-complexity of multilayered DNA strand displacement circuits. In International Workshop on DNA-based Computers. Springer, 144–153.
Eric E. Severson, David Haley, and David Doty. 2021. Composable computation in discrete chemical reaction networks. Distrib. Comput. 34, 6 (2021), 437–461.
David Soloveichik, Matthew Cook, Erik Winfree, and Jehoshua Bruck. 2008. Computation with finite stochastic chemical reaction networks. Natur. Comput. 7, 4 (2008), 615–633. DOI:
David Soloveichik, Georg Seelig, and Erik Winfree. 2010. DNA as a universal substrate for chemical kinetics. Proc. Nat. Acad. Sci. 107, 12 (2010), 5393.
Niranjan Srinivas, James Parkin, Georg Seelig, Erik Winfree, and David Soloveichik. 2017. Enzyme-free nucleic acid dynamical systems. Science 358, 6369 (2017).
Marko Vasić, Cameron Chalk, Austin Luchsinger, Sarfraz Khurshid, and David Soloveichik. 2022. Programming and training rate-independent chemical reaction networks. Proc. Nat. Acad. Sci. 119, 24 (2022), e2111552119.
Anderson DJoshi B(2024)Chemical mass-action systems as analog computers: Implementing arithmetic computations at specified speedTheoretical Computer Science10.1016/j.tcs.2024.114983(114983)Online publication date: Nov-2024
ITCS '14: Proceedings of the 5th conference on Innovations in theoretical computer science
Understanding the algorithmic behaviors that are in principle realizable in a chemical system is necessary for a rigorous understanding of the design principles of biological regulatory networks. Further, advances in synthetic biology herald the time ...
Unconventional Computation and Natural Computation
Abstract
We study the model of continuous chemical reaction networks (CRNs), consisting of reactions such as that can transform some continuous, nonnegative real-valued quantity (called a concentration) of chemical species A and B into equal ...
We explore the class of real numbers that are computed in real time by deterministic chemical reaction networks that are integral in the sense that all their reaction rate constants are positive integers. We say that such a reaction network computes a ...
Permission to make digital or hard copies of part or all of this work for personal or classroom use is granted without fee provided that copies are not made or distributed for profit or commercial advantage and that copies bear this notice and the full citation on the first page. Copyrights for third-party components of this work must be honored. For all other uses, contact the owner/author(s).
Anderson DJoshi B(2024)Chemical mass-action systems as analog computers: Implementing arithmetic computations at specified speedTheoretical Computer Science10.1016/j.tcs.2024.114983(114983)Online publication date: Nov-2024
where \(|\mathbf {z}_i(Y) - \mathbf {x}^{\prime }_i(Y)| = |\mathbf {y}^{\prime }(Y) - \mathbf {y}(Y)|/2 \gt \epsilon /2\). Since \((\mathbf {x}^{\prime }_{i})\) is a subsequence of \(\left(\mathbf {x}_{1/n}\right)\), if we take i large enough so \(\mathbf {x}^{\prime }_i\) is \(\mathbf {x}_{1/n_i}\) with \(1/n_i \lt \epsilon /4\), then we see thatbut \(\mathbf {z}_i\) is reachable from \(\mathbf {x}_{1/n_i}\), which by definition can only reach states \(\mathbf {z}\) with \(|\mathbf {z}(Y) - f(\mathbf {x})| \le 1/n_i\), giving a contradiction.□
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 thatfor any \(t \gt M\). This shows that, taking (*) for granted, \({\boldsymbol \rho }\) must converge to \(\mathbf {y}\).
□
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\).
□
denote those input states from which the siphon \(\Omega\) is drainable.
denote those states from which siphon \(\Omega\) is drainable via a single straight line segment.
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}\).
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.□
