Composable Rate-Independent Computation in Continuous Chemical Reaction Networks
Cameron Chalk Niels Kornerup Wyatt Reeves David Soloveichik










Biological regulatory networks depend upon chemical interactions to process information.  
Engineering such molecular computing systems is a major challenge for synthetic biology and related fields.
The chemical reaction network (CRN) model idealizes chemical interactions, allowing rigorous reasoning about the computational power of chemical kinetics.
Here we focus on function computation with CRNs, where we think of the initial concentrations of some species as the input and the equilibrium concentration of another species as the output. Specifically, we are concerned with CRNs that are rate-independent (the computation must be correct independent of the reaction rate law) and composable (f∘ g can be computed by concatenating the CRNs computing f and g).
Rate independence and composability are important engineering desiderata,
permitting implementations that violate mass-action kinetics, or even “well-mixedness", 
and allowing the systematic construction of complex computation via modular design. 
We show that to construct composable rate-independent CRNs, it is necessary and sufficient to ensure that the output species of a module is not a reactant in any reaction within the module.
We then exactly characterize the functions computable by such CRNs as superadditive, positive-continuous, and piecewise rational linear.
Thus composability severely limits rate-independent computation unless more sophisticated input/output encodings are used.






















Composable Rate-Independent Computation in Continuous Chemical Reaction Networks
    Cameron Chalk, Niels Kornerup, Wyatt Reeves, David Soloveichik
The authors were with the University of Texas at Austin.C. Chalk and D. Soloveichik were supported in part by National Science Foundation grants CCF-1618895 and CCF-1652824.

    December 30, 2023
=================================================================================================================================================================================================================================================







§ INTRODUCTION

A ubiquitous form of biological information processing occurs in complex chemical regulatory networks in cells.
The formalism of chemical reaction networks (CRNs) has been widely used for modelling the interactions underlying such natural chemical computation. 
More recently CRNs have also become a useful model for designing synthetic molecular computation.
In particular, DNA strand displacement cascades can in principle realize arbitrary CRNs, thus motivating the study of CRNs as a programming language <cit.>.
The applications of synthetic chemical computation include reprogramming biological regulatory networks, as well as embedding control modules in environments that are inherently incompatible with traditional
electronic controllers for biochemical, nanotechnological, or medical applications.


The study of information processing within biological CRNs, as well the engineering of CRN functionality in artificial systems, motivates the exploration of the computational power of CRNs.
In general, CRNs are capable of Turing universal computation <cit.>; however, we are often interested in restricted classes of CRNs which may have certain desired properties.
Previous work distinguished two programmable features of CRNs: the stoichiometry of the reactions and the rate laws governing the reaction speeds <cit.>.
As an example of computation by stoichiometry alone, 
consider the reaction 2X → Y.
We can think of the concentrations of species X and Y to be the input and output, respectively.
Then this reaction effectively computes f(x) = x/2, as in the limit of time going to infinity, the system converges to producing one unit of Y for every two units of X initially present.
The reason we are interested in computation via stoichiometry is that it is fundamentally rate-independent,
requiring no assumptions on the rate law (e.g., that the reaction occurs at a rate proportional to the product of the concentrations of the reactants).
This allows the computation to be correct independent of experimental conditions such as temperature, chemical background, or whether or not the solution is well-mixed.


Computation does not happen in isolation.
In an embedded chemical controller, inputs would be produced by other chemical systems, and outputs would affect downstream chemical processes.
Composition is easy in some systems (e.g. digital electronic circuits can be composed by wiring the outputs of one to the inputs of the other).
However, in other contexts composition presents a host of problems. 
For example, the effect termed retroactivity, which results in insufficient isolation of modules, has been the subject of much research in synthetic biology <cit.>. 
In this paper, we attempt to capture a natural notion of composable rate-independent computation, and study whether composability restricts computational power.


.375

    (a)
    
    X_1 + X_2    → Y


.48

    (b)
    
    X_1    → Z_1 + Y 
    
    X_2    → Z_2 + Y 
    
    Z_1 + Z_2    → K 
    
    Y + K    →∅.











Above, we see two examples of rate-independent computation.
Example (a) shows y = min(x_1,x_2). The amount of Y eventually produced will be the minimum of the initial amounts of X_1 and X_2, since the reaction will stop as soon as the first reactant runs out. 
Example (b) shows y = max(x_1,x_2). 
(The last reaction generates waste species that are not used in the computation; thus, we describe the products as an empty set.)
The amount of Y eventually produced in reactions 1 and 2 is the sum of the initial amounts of X_1 and X_2. The amount of K eventually produced in reaction 3 is the minimum of the initial amounts of X_1 and X_2. Reaction 4 subtracts the minimum from the sum, yielding the maximum.

Now consider how rate-independent computation can be naturally composed.

Suppose we want to compute min(min(x_1,x_2),x_3).
It is easy to see that simple concatenation of two min modules (with proper renaming of the species) correctly computes this function:

    X_1 + X_2     Y 
    
        Y + X_3     Y',

where Y' represents the output of the composed computation.
In contrast, suppose we want to compute min(max(x_1,x_2),x_3).
Concatenating the modules yields:

    X_1    → Z_1 + Y 
    
        X_2    → Z_2 + Y 
    
        Z_1 + Z_2    → K 
    
        Y + K    →∅
    
        Y + X_3    → Y',

where Y' represents the output of the composed computation.
Observe that depending on the relative rates of reactions 4 and 5, the eventual value of Y' will vary between min(max(x_1,x_2),x_3) and min(x_1+x_2,x_3),
and the composition does not compute in a rate-independent manner.

Why is min composable, but max not? 
The problem arose because the output of the max module (Y) is consumed in both the max module and in the downstream min module. 
This creates a competition between the consumption of the output within its own module and the downstream module.


Towards modularity, we assume the two CRNs to be composed do not share any species apart from the interface between them (i.e., a species Y representing the output of the first network is used as the species representing the input to the second network, and otherwise the two sets of species are disjoint).
We prove that to construct composable rate-independent modules in this manner, it is necessary and sufficient to ensure that the output species of a module is not a reactant in any reaction of that module.
We then exactly characterize the computational power of composable rate-independent computation.

Previously it was shown that without the composability restriction, rate-independent CRNs can compute arbitrary positive-continuous, piecewise rational linear functions <cit.>.
Positive-continuity means that the only discontinuities occur when some input goes from 0 to positive, 
and piecewise rational linear means that the function can by defined by a finite number of linear pieces (with rational coefficients).
Note that non-linear continuous functions can be approximated to arbitrary accuracy.[To approximate arbitrary continuous non-linear functions, piecewise linear functions are not sufficient, but rather we need piecewise affine functions (linear functions with offset). However, affine functions can be computed if we use an additional input fixed at 1.]


We show that requiring the CRN to be composable 
restricts the class of computable functions to be superadditive functions; i.e., functions that satisfy: for all input vectors , ,  f() + f() ≤ f( + ).
This strongly restricts computational power: for example, subtraction or max cannot be computed or approximated in any reasonable sense.
In the positive direction, we show that any
superadditive, positive-continuous, piecewise rational linear function 
can be computed by composable CRNs in a rate-independent manner. 
Our proof is constructive, and we further show that unimolecular and bimolecular reactions are sufficient.


We note that different input and output encodings can change the computational power of rate-independent, composable CRNs.
For example, in the so-called dual-rail convention, input and output values are represented by differences in concentrations of two species (e.g., the output is equal to the concentration of species Y^+ minus the concentration of Y^-).
Dual-rail simplifies composition—instead of consuming the output species to decrease the output value, a dual-rail CRN can produce Y^-—at the cost of greater system complexity.
Dual-rail CRNs can compute the full class of continuous, piecewise rational linear functions while satisfying rate-independence and composability <cit.>.
Note, however, that the dual-rail convention moves the non-superadditive subtraction operation to “outside” the system, and converting from a dual-rail output to a direct output must break composability.
































§ PRELIMINARIES


Let  and  denote the set of nonnegative integers and the set of real numbers, respectively. The set of the first n positive integers is denoted by [n].
Let _≥ 0 be the set of nonnegative real numbers, and similarly _> 0 be the set of positive real numbers.


If Λ is a finite set (in this paper, of chemical species), we write ^Λ to denote the set of functions f:Λ→, and similarly for ^Λ, ^Λ, etc.
Equivalently, we view an element ∈ A^Λ as a vector of |Λ| elements of A, each coordinate “labeled” by an element of Λ.
Given a function f : A → B, we use f|_C to denote the restriction of f to the domain C. We also use the notation ↾Δ to represent  projected onto ^Δ. Thus, ↾Δ = 0⃗ iff (∀ S∈Δ) (S)=0.

If Δ⊆Λ, we view a vector ∈^Δ equivalently as a vector ∈^Λ by assuming (S)=0 for all S ∈Λ∖Δ.



 §.§ Chemical reaction networks

We will start by defining the notation used to describe chemical reactions.













Given a finite set of chemical species Λ, a reaction over Λ is a pair α = ⟨,⟩∈^Λ×^Λ, specifying the stoichiometry of the reactants and products, respectively.[As we are studying CRNs whose output is independent of the reaction rates, we leave the rate constants out of the definition.] 




In this paper, we assume that ≠0⃗, i.e., we have no reactions of the form ∅→….


For instance, given Λ={A,B,C}, the reaction A+2B → A+3C is the pair (1,2,0)(1,0,3).



A (finite) chemical reaction network (CRN) is a pair =(Λ,R), where Λ is a finite set of chemical species,

and R is a finite set of reactions over Λ. 


Next we map language about chemical reaction networks to formal definitions and notation.


A state of a CRN =(Λ,R) is a vector ∈^Λ.



For any ∈^Λ and any S∈Λ, (S) is the concentration of S in .


For any ∈^Λ, the set of species present in  (denoted by []) is {S ∈Λ | (S) > 0 }.



Given a state  and reaction α=, we say that α is applicable in  if [] ⊆ [] (i.e.,  contains positive concentration of all of the reactants).



A reaction produces (consumes) a species S if S appears as a product (reactant).[Note that typically a catalyst species is not considered consumed nor produced by a catalytic reaction. 
For simplicity, as our definition suggests, we say it is both produced and consumed.]





 §.§ Reachability and stable computation

We now follow <cit.> in defining rate-independent computation in terms of reachability between states (this treatment is in turn based on the notion of “stable computation” in distributed computing <cit.>).

Intuitively, we say a state is “reachable” if some rate law can take the system to this state. 
For computation to be rate-independent, 
since unknown rate laws might take the system to any reachable state, the system must be able to reach the correct output from any such reachable state.







To define the notion of reachability,

a key insight of <cit.> allows one to think of reachability via a sequence of straight line segments.
This may be unintuitive, since mass-action[Although the formal definition of mass-action kinetics is outside the scope of this paper, we remind the reader that a CRN with rate constants on each reaction define a system of ODEs under mass-action kinetics.
For example, the two reactions
A + B → A + C and
C + C → B correspond to the following ODEs:
    ȧ = 0,
    ḃ = k_2c^2 - k_1ab, and
    ċ = k_1ab - 2k_2c^2,
where a, b, and c are the concentrations of species A, B, and C over time and k_1, k_2 are the rate constants of the reactions.]  and other rate laws trace out smooth curves.
However, a number of properties are shown which support straight-line reachability as an interpretation which includes mass-action 
reachability as well as reachability  under other rate laws. 






Let  be a CRN defined by (Λ, R). The linear transformation : ^R →^Λ that maps from the unit vector representing a reaction to the net change in species caused by that reaction is the stoichiometry matrix for . 















Note that we can intuitively think of  being a matrix where the columns represent the net change in species caused by each reaction. Under this representation, observe that entries in  will be negative when more of a reactant is consumed than is produced in a reaction. Observe that the image of  represents the possible changes in a state that can occur via the reactions in R. We will formalize this notion with the next few definitions.

For a CRN with the reactions R, we say that any vector ∈^R is a flux vector. We use [] to denote the set {r | (r) > 0}. We say that  is applicable at a state  if every reaction in [] is applicable at .


For a CRN with species Λ and stoichiometry matrix , we say a state ∈_≥ 0^Λ is straight-line reachable from , written , or more precisely as →_, if there is a applicable flux vector  such that +  =.


Intuitively, a single segment means running the reactions applicable at  at a constant (possibly 0) rate specified by  to get from  to . Since applying a flux vector can change the set of species present, a⃗b⃗ does not imply that a⃗ and b⃗ have the same set of applicable reactions. Therefore there can be a state c⃗ that is straight-line reachable from b⃗ but not from a⃗. This leads us to our next definition.

We say state  is 1-segment reachable from  if it is straight line reachable. We say a state  is l-segment reachable if there is a state ' that is (l-1)-segment reachable from  such that '.

Generalizing to an arbitrary number of segments, we obtain our general notion of reachability below.
Note that by the definition of straight-line reachability, only applicable reactions occur in each segment.
The definition of reachability is closely related to exploring the “stoichiometric compatibility class” of the initial state <cit.>.



A state  is reachable from , written →, if ∃ l∈ such that  is l-segment reachable from .
We denote the set of states reachable from c⃗, i.e., {d⃗ |  c⃗→d⃗}, as Post(c⃗).















We think of state  as being reachable from state  if there is a “reasonable” rate law that takes the system from  to .
As desired, previous work showed that if state  is reached from  via a mass-action trajectory, it is also segment-reachable.


  
   If  is mass-action reachable from , then →.



We can now use reachability to  formally define rate-independent computation.








A chemical reaction computer (CRC) is a tuple = (Λ,R,Σ,Y), where (Λ,R) is a CRN, Σ⊂Λ, written as Σ = {X_1,…,X_n}, is the set of input species, and Y ∈Λ∖Σ is the output species.

For simplicity, assume a canonical ordering of Σ={X_1,…,X_n} so that a vector ∈^n (i.e., an input to f) can be viewed equivalently as a state ∈^Σ of  (i.e., an input to ).
















A state ∈^Λ is output stable if, for all ' such that →', (Y) = '(Y),  i.e., once  is reached, no reactions can change the concentration of the output species Y.



 
Let f:^n → be a function and let  be a CRC.
We say that  stably computes f if, for all ∈^n and all  such that →, there exists an output stable state  such that → and (Y) = f().


We can intuitively justify the above definition of reachability and stable computation as capturing the class of computation that is independent of the rate law.
The output stable states are exactly those in which the output cannot be changed by a rate law chosen by an adversary.
If a chemical reaction network does not stably compute a function, 
then some rate law can take the system to a state from which an output stable state is not reachable (including by mass-action by Lemma <ref>).

The results herein extend easily to functions f:^n →^l, i.e., whose output is a vector of l real numbers.
This is because such a function is equivalently l separate functions f_i:^n→.

Also note that initial states contain only the input species Σ; other species must have initial concentration 0.
Section <ref> discusses how allowing some initial concentration of non-input species affects computation.






 §.§ Composability

In this section we define the composition of CRCs and formally relate composability to a CRC not using its output species as a reactant (output-obliviousness).
We show that output-oblivious CRCs are composable, and that any composable CRC can be reduced (simply by removing reactions) to an output-oblivious form.

We define the composition of two CRCs intuitively
as the concatenation of their chemical reactions, such that the output species of the first is the input species of the second:

Given two CRCs _1 = (Λ_1,R_1,Σ_1,Y_1) and _2 = (Λ_2,R_2,Σ_2,Y_2), consider _2' = (Λ_2',R_2',Σ_2',Y_2') constructed by renaming species of _2 such that Λ_1∩Λ_2' = {Y_1} and Y_1 ∈Σ_2'.
The composition of _1 and _2 is the CRC  _2∘1 = (Λ_1∪Λ_2', R_1∪ R_2', Σ_1 ∪Σ_2' ∖{Y_1}, Y_2').
In other words, the composition is constructed by concatenating _1 and _2 such that their only interface is the output species of _1, used as the input for _2.


We say two CRCs are composable if they stably compute the composition of their functions when composed:


A CRC _1 which stably computes f_1 is composable if ∀_2 stably computing f_2, _2∘1 stably computes f_2∘ f_1.



We want to relate composability to the property that a CRC does not use its output species as a reactant:

We call a CRC (Λ,R,Σ,Y) output-oblivious if Y does not appear as a reactant in R.





For simplicity, we focus on single-input, single-output CRCs, but these results easily generalize to multiple input and output settings.

For the proof that the output-oblivious condition is sufficient to guarantee composability, we formalize the idea that the composed CRCs act independently, and do not interfere with each other's execution.
In Lemmas <ref> and <ref> we show how this independence can be used to “reorder” the sequence of reactions of a CRC in way that preserves the concentrations in the final state.
In Lemma <ref>, we take an output-oblivious CRC, compose another CRC downstream, and reorder any sequence of reactions of the composed CRC into a sequence which we can easily argue must have stably computed as expected.














A flux vector u⃗ produces (consumes) a species S if there is an ∈ [u⃗] such that S is a product (reactant) of .


A flux vector u⃗_1 is independent of a flux vector u⃗_2 if u⃗_1 does not consume any species that are produced or consumed by u⃗_2. 


In a CRC = (Λ, R) if flux vector u⃗_1 is independent of flux vector u⃗_2 then:

    
  * If →__2→__1, then →__1 + _2.
    
  * If →__1 + _2, then there is a state  such that →__1→__2.
    
  * If →__1 and →__2, then there is a state  such that →__1.








For 1, since u⃗_1 is independent of u⃗_2, we know that u⃗_2 cannot produce any of the species necessary to make u⃗_1 applicable. Since u⃗_1 was applicable at b⃗, we know that u⃗_1 must be applicable at a⃗, so a⃗→_u⃗_1 + u⃗_2c⃗.

For 2, we need to show that ∈^Λ (has nonnegative concentrations)
and that _2 is applicable at .
Consider any species S ∈Λ such that S is not produced in _2. Since =  - _2 and (S) ≥ 0, we know that S has nonnegative concentration at . Now consider any species S ∈Λ such that S is produced in _2. Since _2 produces S, we know that _1 must not consume it because u⃗_1 is independent of u⃗_2. Since =  + _1, we can conclude that S must have nonnegative concentration at . 
Therefore we can conclude that ∈_≥ 0^Λ.
To see that _2 is applicable at , first observe that 
_2 is applicable at . Then, since _1 is independent of _2, it follows that  _2 must also be applicable at . 

For 3,
we want to show that _1 is still applicable at  and that =  + M_1 is in ^Λ (has nonnegative concentrations).
_1 is still applicable at  since _1 does not consume species involved in _2.
If  had negative concentration on species S, that species must have been consumed by _1 since  has nonnegative concentrations.
Since _2 does not consume any species consumed by _1 and =  + _1, then negative (S) implies negative (S), which contradicts that ∈^Λ, so ∈^Λ.














Given two output-oblivious CRCs _1 and _2,
consider the composition CRC _2∘1.
If c⃗→d⃗ then there is b⃗ such that c⃗→b⃗→d⃗, where c⃗→b⃗ only uses reactions from _1 and b⃗→d⃗ only uses reactions from _2.



Let {v⃗_1, …, v⃗_n} be the flux vectors such that c⃗→_v⃗_1c⃗_1 →_v⃗_2…→_v⃗_nd⃗. We can write v⃗_i = u⃗_1, i + u⃗_2, i, where u⃗_1,i corresponds to the reactions in _1 and u⃗_2, i corresponds to the reactions in _2. Since _1 is output-oblivious, we know that every _1,i is independent of every _2,j and thus we can apply Lemma <ref> item 2 to see that c⃗→_u⃗_1, 1b⃗_1 →_u⃗_2,1c⃗_1 →_u⃗_1, 2b⃗_2 …→_u⃗_2, nd⃗. By repeatedly applying Lemma <ref> items 1 and 2, we can then rearrange the sequence of reactions so that each u⃗_1, i precedes each u⃗_2, j to get c⃗→_u⃗_1, 1b⃗_1 →_u⃗_1,2…→_u⃗_1,nb⃗→_u⃗_2, 1…→_u⃗_2, nd⃗. 



Output-oblivious
CRCs are composable.


Consider the composition _2∘1 =(Λ,R,Σ,Y) of two CRCs _1=(Λ_1,R_1,Σ_1,Y_1) and _2=(Λ_2,R_2,Σ_2,Y_2) that stably compute f_1 and f_2 respectively, and consider an input x ∈_≥0^Σ. Consider some state c⃗ reached from x⃗ in _2 ∘ 1. Let b⃗ be as in Lemma <ref>, so  →_u⃗_1,1…→_u⃗_1,nb⃗→_u⃗_2,1…→_u⃗_2,n , where r⃗_i = (_i, 1, …, _i, n) is a sequence of flux vectors with ⋃_j [_i,j] ⊆ R_i. Since _1 stably computes f_1, we know that there is some _1-output stable state o⃗_1 reachable from b⃗ using a series of flux vectors r⃗ = (_1, _2 …_k) such that ⋃_i[_i] ⊆ R_1. Since r⃗_2 only uses reactions from _2 and _1 is output-oblivious, every flux vector in r⃗ is independent of every flux vector in r⃗_2, so by Lemma <ref> item 3 we know the sequence of flux vectors  is applicable starting at c⃗. Let a⃗ be such that c⃗→_u⃗_1…→_u⃗_ka⃗. Then applying Lemma <ref> items 1 and 2 repeatedly to the flux vectors in r⃗ and r⃗_2, we see that o⃗_1 →_u⃗_2,1…→_u⃗_2,na⃗. Since _2 stably computes f_2, since o⃗_1(Y_1) = f(x), and since a⃗ is reachable from o⃗_1 only using reactions in _2, there must be some o⃗_2 that is _2-output stable such that a⃗→o⃗_2 and o⃗_2(Y_2) = f_2 ∘ f_1(x). We know that o⃗_2 is reachable from c⃗ since c⃗→a⃗→o⃗_2. Finally, since o⃗_1 is _1-output stable, reactions from _1 cannot change the concentrations of species in o⃗_2↾Λ_2, so if o⃗_2 →y⃗, then restricting to _2 we find o⃗_2↾Λ_2→__2y⃗↾Λ_2. Since o⃗_2 is _2-output stable we see that y⃗(Y_2) = o⃗_2(Y_2), so o⃗_2 is _2 ∘ 1-output stable.






Next we show that the output-oblivious condition is effectively necessary for composition.
Technically, there are CRCs which are not output-oblivious but are composable. However, we show that for such CRCs, we can remove reactions until they are output-oblivious, resulting in a CRC which is still composable and computes the same function.
Thus, characterizing what is computable by output-oblivious CRCs does characterize the class of functions computable by composable CRCs.

To prove this, we will want to reach a state that has used up its capacity to produce more of some species S. Intuitively this can be done by producing the maximum amount of S possible by the CRC. 
However, in general this notion is ill-defined.[In the CRN X → Z, X + Z → S + Z, from an initial state  consisting of only species X, note that max_ | →(S) is not well-defined.
To see this, note that some ϵ amount of X must convert to Z to catalyze the second reaction; we can always choose a smaller ϵ' to produce more S.
]
Lemma <ref> proves that there is a state with maximal amount of S. In other words, 
from any state we can always reach a state where afterwards it is impossible to increase the amount of S. The proof of Lemma <ref> is left to the appendix.


For any state c⃗ and any species S, if the amount of S present in any state reachable from c⃗ is bounded above, there is a state d⃗ reachable from c⃗ such that for any state a⃗ reachable from d⃗, we know that a⃗(S) ≤d⃗(S). 





If a CRC  stably computes f and is composable, then we can remove all reactions where the output species appears as a reactant, and the resulting output-oblivious CRC will still stably compute f.



Let _1 = (Λ_1,R_1,Σ_1,Y_1) be a composable CRC stably computing some function f. Let _0 ⊆_1 be the CRN obtained by removing all of the reactions that consume Y_1 from _1. We would like to show that _0 stably computes f. Suppose we compose _1 with _2 consisting of only the reaction Y_1 → Y_2 with input species Y_1 and output species Y_2.

Since _1 is composable and _2 stably computes the identity function, the resulting CRN _2∘1 must stably compute f. For any input vector x⃗ consider a state c⃗ reachable from x⃗.

Assume that _0 could reach a state d⃗ from c⃗ where d⃗(Y_1) > f(x⃗). Then _1 would not be composable because _2∘1 can also reach d⃗ and then applying the reaction in _2 to convert all Y_1 into Y_2 gives us a state d⃗'⃗ with d⃗'⃗(Y_2) > f(x⃗). Since there is no reaction in _2∘1 that consumes Y_2 there is no output stable state reachable from d⃗'⃗ that computes f(x⃗). Therefore there is no state d⃗ such that d⃗(Y_1) ≥ f(x⃗) and 
d⃗ is reachable from c⃗ via reactions of _0. 


Since this implies that the amount of Y_1 in any state reachable from c⃗ is bounded, we can apply Lemma <ref> to say that there is a state d⃗ such that c⃗→d⃗ and for any state a⃗ reachable from d⃗ we know a⃗(Y_1) ≤d⃗(Y_1). Since _0 has no reactions that consume Y_1, this is an output stable state of _0. Now, consider the state b⃗ in _2∘1 obtained by converting all Y_1 in d⃗ into Y_2. Observe that if there were a way to produce Y_1 from b⃗, then there would be a state in _0 reachable from d⃗ that contained more Y_1. Since there are no reactions in _2∘ 1 that consume Y_2 and no reactions that produce Y_1, we know that b⃗ is an output stable state. Since _2 ∘ 1 stably computes f(x⃗), we know b⃗(Y_2) = f(x⃗). Thus we can conclude that d⃗(Y_1) = f(x⃗) and _0 stably computes f(x⃗).





To allow composition of multiple downstream CRCs, we can use the reaction Y → Y_1 + … + Y_n to generate n “copies” of the output species Y, such that each downstream module uses a different copy as input.
Additionally, if the downstream module is output-oblivious, then the composition is also output-oblivious and thus the composition is composable.
These observations allow complex compositions of modules, and will be used in our constructions in Section <ref>.





§ FUNCTIONS COMPUTABLE BY COMPOSABLE CRNS

Here we give a complete characterization of the functions computable by composable CRNs as superadditive, positive-continuous, and piecewise rational linear.



A function f : ^n →^l is superadditive iff ∀a⃗, b⃗∈^n, f(a⃗) + f(b⃗) ≤ f(a⃗ + b⃗).


Note that superadditivity implies monotonicity in our case, since the functions computed must be nonnegative.
As an example, we show that the max function is not superadditive:


The function max(x_1, x_2) is not superadditive.


Pick any x_1, x_2 > 0. Observe that max(x_1,0) + max(0,x_2) = x_1 + x_2. But since x_1 and x_2 are both positive, we know that x_1 + x_2 > max(x_1,x_2). Thus max is not superadditive and by Lemma <ref> there is no composable CRN which stably computes max.



A function f : ^n_≥0→^l is positive-continuous if for all U ⊆ [n], f is continuous on the domain D_U = { x∈^n_≥0 | (∀ i ∈ [n]),  x(i) > 0  i ∈ U }.
I.e., f is continuous on any subset D ⊂^n_≥0 that does not have any coordinate i∈[n] that takes both zero and positive values in D.


Next we give our definition of piecewise rational linear.
One may (and typically does) consider a restriction on the domains selected for the pieces, however this restriction is unneccesary in this work, particularly because the additional constraint of positive-continuity gives enough restriction.


A function f : ^n → is rational linear if there exists a_1,…, a_n ∈ℚ such that f(x) = ∑_i=1^n a_ix(i).
A function f : ^n → is piecewise rational linear if there is a finite set of partial rational linear functions f_1,…,f_p : ^n → with ⋃_j=1^p  f_j = ^n, such that for all j ∈ [p] and all x ∈ f_j, f(x) = f_j(x). We call f_1,…,f_p the components of f.


The following is an example of a superadditive, positive-continuous, piecewise rational linear function:

rCl
f(x⃗) = 
      x_1 + x_2     x_3 > 0 

      min(x_1,x_2)     x_3 = 0.
   


The function is superadditive since for all input vectors a⃗ = (a_1, a_2, a_3), b⃗ = (b_1, b_2, b_3), there are three cases: (1) a_3 = b_3 = 0, in which case both input vectors compute min which is a superadditive function; (2) a_3, b_3 ≠ 0, in which case both input vectors compute x_1 + x_2, which is a superadditive function; (3) without loss of generality, a_3 = 0 and b_3 ≠ 0, in which case f(a⃗) + f(b⃗) = min(a_1, a_2) + b_1 + b_2 ≤ a_1 + a_2 + b_1 + b_2 = f(a⃗ + b⃗).
The function is positive-continuous, since the only points of discontinuity are when x_3 changes from zero to positive.
The function is piecewise rational linear, since min is piecewise rational linear.


A function f : ℝ^n_≥0→ℝ_≥0 is computable by a composable CRC if and only if it is superadditive positive-continuous piecewise rational linear.

We prove each direction of the theorem independently in Sections <ref> and <ref>.



 §.§ Computable Functions are Superadditive Positive-Continuous Piecewise Rational Linear


Here, we prove that a stably computable function must be superadditive positive-continuous piecewise rational linear.
The constraints of positive-continuity and piecewise rational linearity stem from previous work:

[Proven in <cit.>]
If a function f : ℝ^n_≥0→ℝ_≥0 is stably computable by a CRC, then f is positive-continuous piecewise rational linear.


In addition to the constraints in the above lemma, we show in Lemma <ref> that a function must be superadditive if it is stably computed by a CRC.
To prove this, we first note a useful property of reachability in CRNs. 

Given states a⃗,b⃗,c⃗, if a⃗→b⃗ then a⃗ + c⃗→b⃗ + c⃗.

Adding species cannot prevent reactions from occurring.
Thus, we can consider the series of reactions where c⃗ doesn't react to reach the state b⃗ + c⃗ from the state a⃗ + c⃗.

We now utilize this lemma to prove that composably computable functions must be superadditive.


If a function f : ℝ^n_≥0→ℝ_≥0 is stably computable by a composable CRC, then f is superadditive.


Assume 𝒞 stably computes f.

By definition of 𝒞 stably computing f, ∀ initial states x⃗_1, x⃗_⃗2⃗, ∃ o⃗_1, o⃗_2 such that x⃗_1 →o⃗_1 with o⃗_1(Y) = f(x⃗_1) and x⃗_2 →o⃗_2 with o⃗_2(Y) = f(x⃗_2).
Consider 𝒞 on input x⃗_1 + x⃗_2.
By Lemma <ref>, x⃗_1 + x⃗_2 →o⃗_1 + x⃗_2, and again by Lemma <ref>, o⃗_1 + x⃗_2 →o⃗_1 + o⃗_2.
Looking at the concentration of output species Y, we have (o⃗_1 + o⃗_2)(Y) = f(x⃗_1) + f(x⃗_2).
Since  stably computes f, there exists an output stable state o⃗' reachable from initial state x⃗_1 + x⃗_2 and reachable from state o⃗_1 + o⃗_2, with o⃗'⃗(Y) = f(x⃗_1 + x⃗_2).
Since 𝒞 is composable, we can assume that species Y does not appear as a reactant without loss of generality by Lemma <ref>. Thus the concentration of Y in any state reachable from state o⃗_1 + o⃗_2 cannot be reduced from f(x⃗_1) + f(x⃗_2), implying o⃗'⃗(Y) = f(x⃗_1 + x⃗_2) ≥ f(x⃗_1) + f(x⃗_2).
This holds for all input states x⃗_1, x⃗_2, and thus f is superadditive.




No composable CRC computes f(x_1, x_2) = max(x_1,x_2).






 §.§ Superadditive Positive-Continuous Piecewise Rational Linear Functions are Computable


It was shown in <cit.> that every piecewise linear function can be written as a max of mins of linear functions. This fact was exploited in <cit.> to construct a CRN that dual-rail 

computed continuous piecewise rational linear functions. To directly compute a positive-continuous piecewise rational linear function, dual-rail networks were used to compute the function on each domain, take the appropriate max of mins, and then the reaction Y^+ + Y^- →∅ was used to convert the dual-rail output into a direct output where the output species is Y^+. However, this technique is not usable in our case: by Corollary <ref>, we cannot compute the max function, and the technique of converting dual-rail output to a direct output is not output-oblivious. In fact, computing f(Y^+, Y^-) = Y^+ - Y^- is not superadditive, and so by Lemma <ref>, there is no composable CRC which computes this conversion.

Since our functions are positive-continuous, we first consider domains where the function is continuous, and show that it can be computed by composing rational linear functions with min.
Since rational linear functions and min can be computed without using the output species as a reactant, we achieve composability.
We then extend this argument to handle discontinuities between domains.



An open ray ℓ in ℝ^n from the origin through a point x⃗ is the set ℓ = y⃗∈ℝ^n  | y⃗ = t·x⃗, t ∈ℝ_>0. Note that t is strictly positive, so the origin is not contained in ℓ.



We call a subset D ⊆ℝ^n a cone if for all x⃗∈ℝ^n, we know that x⃗∈ D implies the open ray from the origin through x⃗ is contained in D.



Suppose we are given a continuous piecewise rational linear function f: ℝ_>0^n →ℝ_≥ 0. Then we can choose domains for f which are cones which contain an open ball of non-zero radius.


Intuitively, we can consider any open ray from the origin and look at the domains for f along this ray. If the ray traveled through different domains, then there must be boundary points where the function switches domains. But we know that f is continuous, so the domains must agree on their boundaries. Since there is only one line that passes through the origin and any given point, the domains must share the same linear function to be continuous. Thus we can place the ray into one domain corresponding to its linear function.
Applying this argument to all rays gives these domains as cones. This argument is formalized in a proof in the appendix.


















Any superadditive continuous piecewise rational linear function f: ℝ_>0^n →ℝ_≥ 0 can be written as the minimum of a finite number of rational linear functions g_i.



Since f is a continuous piecewise rational linear function, by Lemma <ref>, we can choose domains D_i_i = 1^N for f which are cones and contain an open ball of non-zero radius, such that f|_D_i = g_i|_D_i, where g_i is a rational linear function. Now pick any x⃗∈ℝ_>0^n and any g_j. Then because D_j is a cone containing an open ball of finite radius, it contains open balls with arbitrarily large radii. In particular, it contains a ball with radius greater than |x⃗|, so there exist points y⃗, z⃗∈ D_j such that y⃗ + x⃗ = z⃗. By the superadditivity of f, the linearity of g_j, and the fact that y⃗, z⃗∈ D_j, we see: 

    g_j(y⃗) + f(x⃗) = f(y⃗) + f(x⃗) ≤ f(z⃗) = g_j(x⃗ + y⃗) = g_j(y⃗) + g_j(x⃗)


so that f(x⃗) ≤ g_j(x⃗). Since this is true for all g_j, and since we know that f(x⃗) = g_i(x⃗) for some i, we see that f(x⃗) = min_i g_i(x⃗), as desired.


Lemma <ref> is particularly useful for us since, as seen in the introduction, CRCs computing min are easy to construct, and rational linear functions are relatively straightforward as well.
The next lemma gives details on constructing a CRC to compute f by piecing together CRCs which compute the components (rational linear functions) of f and then computing the min across their outputs.
However, since Lemma <ref> as given applies to continuous functions with domain ^n_>0, so does this lemma; we handle the domain ^n_≥0 later on.


We can construct a composable CRC that stably computes any superadditive continuous piecewise rational linear function f : ℝ^n_>0→ℝ_≥ 0.


By Lemma <ref>, we know that f can be written as the minimum of a finite number of rational linear functions g_i. Observe that a general rational linear function g(x⃗) = a_1 x_1 + a_2 x_2 + … a_n x_n is stably computed by the reactions 
    ∀ i,  k_iX_i → a_ik_i Y
 where k_i is a positive integer such that k_ia_i is also a positive integer. Since f is the minimum of a number of g_i's, we can make a chemical reaction network where we compute each g_i using a copy of the input species, calling the output Y_i (the reaction X_1 → X_1^1 + … + X_1^5 produces five species with concentrations equal to X_1's initial concentration, effectively copying the input species so that the input may be a reactant in several modules without those modules competing). Next, we use the chemical reaction 
    Y_1 + … + Y_n → Y
 to get the minimum of the Y_i's. Since each Y_i obtains the count of the corresponding g_i, this CRN will produce the minimum of the g_i's quantity of Y's. Thus, according to Lemma <ref>, the described CRC stably computes f.
Note that each sub-CRC described in this construction is output-oblivious, and thus composable, so the composition of these modules maintains correctness.


The above construction only handles the domain ^n_>0, where we know our functions are continuous by positive-continuity. However, when extended to the domain ^n_≥0, positive-continuity of our functions allows discontinuity where inputs change from zero to positive.
The challenge, then, is to compute the superadditive continuous piecewise rational linear function corresponding to which inputs are nonzero.




Surprisingly, Lemma <ref> below shows that we can express a superadditive positive-continuous piecewise rational linear function as a min of superadditive continuous piecewise rational linear functions.
The first step towards this expression is to see that, given two subspaces of inputs wherein the species present in one subspace A are a superset of the species present in a subspace B, the function as defined on the subspace A must be greater than the function as defined on the subspace B; otherwise, the function would disobey monotonicity and thus superadditivity, as proven below:




Consider any superadditive positive-continuous piecewise rational linear function f: ℝ_≥0^n →ℝ_≥ 0. Write N = [n], and for each S ⊆ N, let g_S(x⃗) be the superadditive continuous piecewise rational linear function that is equal to f on D_S. If S, T ⊆ N and S ⊆ T, then for all x⃗∈ D_S we know g_S(x⃗) ≤ g_T(x⃗).



Write e⃗_i for the vector of length 1 pointing in the positive direction of the ith coordinate axis. Define the vector v⃗ = ∑_i ∈ T ∖ Se⃗_i. Then for any x⃗∈ D_S and any ϵ∈ℝ_>0, we know that x⃗ + ϵv⃗∈ D_T. Since f is superadditive, it is also monotonic. Suppose that g_T(x⃗) < g_S(x⃗). Because g_T is continuous, taking δ = g_S(x⃗) - g_T(x⃗) > 0, there is some small enough ϵ > 0 such that 
    f(x⃗ + ϵv⃗) = g_T(x⃗ + ϵv⃗) < g_T(x⃗) + δ = g_S(x⃗) = f(x⃗)

contradicting the monotonicity of f. Our assumption must be false, so g_S(x⃗) ≤ g_T(x⃗).












Next we define a predicate for each subset of inputs which is true if all inputs in that subset are nonzero. 
Intuitively, in the CRC construction to follow, this predicate is used by the CRC to determine which inputs are present:


For any set S ⊆ [n], define the S-predicate P_S: ℝ^n_≥ 0→0,1 to be the function given by:


    P_S(x⃗) = 
    1    x⃗(i) > 0 ∀ i ∈ S 
    
    0    otherwise.



A naïve approach might be the following: for each subdomain D_S, the function is continuous, so compute it by CRC according to Lemma <ref>, producing an output Y_S.
Then compute the P_S predicate by CRC, and if the predicate is true (e.g., a species representing P_S has nonzero concentration), use that species to catalyze a reaction which changes the Y_S to Y, the final output of the system.
However, note that if T is a subset of S, P_S and and P_T are both true, so this technique will overproduce Y.

The following technique solves this issue by identifying a min which can be taken over the intermediate outputs Y_S.
In particular, for each S, we compute g_S(x⃗) + ∑_K ⊈S P_K(x⃗) g_K(x⃗), and then take the min of these terms.
When S corresponds to the set of input species with initially nonzero concentrations, then the summation term in this expression is 0, since P_K(x) = 0 for all K⊈S.
When S does not correspond to the set of input species with initially nonzero concentration, then either (1) it is a superset of the correct set I, in which case Lemma <ref> says that g_S(x⃗) ≥ g_I(x⃗) (thus the min of these is g_I(x⃗)) or (2) the summation term added to g_S(x⃗) contains at least g_I(x⃗), and since g_S(x) + g_I(x⃗) ≥ g_I(x⃗), the min of these is g_I(x⃗).
Thus taking the min for all S of g_S(x⃗) + ∑_K ⊈S P_K(x⃗) g_K(x⃗) is exactly g_I(x⃗), where I is the correct set of initially present input species.


Consider any superadditive positive-continuous piecewise rational linear function f: ℝ_≥0^n →ℝ_≥ 0. Write N = [n], and for each S ⊆ N, let g_S(x⃗) be the superadditive continuous piecewise rational linear function that is equal to f on D_S. Then, f(x⃗) = min_S ⊆ N[g_S(x⃗) + ∑_K ⊈S P_K(x⃗) g_K(x⃗)].



For S ⊆ N, let h_S: _≥ 0^n →_≥ 0 be given by 


    h_S(x⃗) = g_S(x⃗) + ∑_K ⊈S P_K(x⃗)g_K(x⃗).


We want to show that f(x⃗) = min_S ⊆ Nh_S(x⃗). To do this, fix x⃗∈ℝ_≥ 0^n and define the set
I = i ∈ N  | x⃗(i) > 0. First, let's show that h_I(x⃗) = f(x⃗). By the definition of I, for all K ⊈I, we know P_K(x⃗) = 0. Thus, ∑_K ⊈I P_K(x⃗) g_K(x⃗) = 0, so h_I(x⃗) = g_I(x⃗) = f(x⃗). Now we must show that h_S(x⃗) ≥ f(x⃗) for all S ⊆ N. There are two cases to consider:



Case 1: S ⊉I


In this case,

    h_S(x⃗)    = g_S(x⃗) + ∑_K ⊈S P_K(x⃗) g_K(x⃗)
       ≥ g_S(x⃗) + P_I(x⃗) g_I(x⃗)
       ≥ P_I(x⃗) g_I(x⃗).

By the definition of I, we know P_I(x⃗) = 1, so P_I(x⃗)g_I(x⃗) = g_I(x⃗) = f(x⃗). Thus we get that h_S(x⃗) ≥ f(x⃗).



Case 2: S ⊇ I


By Lemma <ref>, g_S(x⃗) ≥ g_I(x⃗). As a result,

    h_S(x⃗)    = g_S(x⃗) + ∑_K ⊈S P_K(x⃗) g_K(x⃗)
       ≥ g_S(x⃗) ≥ g_I(x⃗)
       = f(x⃗).


Since for all x⃗∈ℝ_≥ 0^n, we know h_S(x⃗) ≥ f(x⃗) for all S ⊆ N and h_I(x⃗) = f(x⃗) for some I ⊆ N, it follows that f(x⃗) = min_S ⊆ N h_S(x⃗).


Lemma <ref> takes the above Lemma <ref> along with the construction which stably computes on strictly continuous domains from Lemma <ref> to construct a CRC which stably computes on positive-continuous domains.


Given any superadditive positive-continuous piecewise rational linear function f: ℝ_≥0^n →ℝ_≥0, there exists a composable CRC which stably computes f.




The proof follows by identifying that the function can be expressed as a composition of functions (via Lemma <ref>) which are computable by output-oblivious CRCs and are thus composable by Lemma <ref>.
By Lemma <ref>, we know that f(x⃗) = min_S ⊆ N[g_S(x⃗) + ∑_K ⊈S P_K(x⃗) g_K(x⃗)]. 
The first subroutine copies the input species, e.g. X_1 → X_1^1 +…+ X_1^5, in order for each sub-CRC to not compete for input species.
This copying is output-oblivious.
Then for any Q ⊆ [n], P_Q(x) is computed using one set of copies via the reaction: 
    ∑_i ∈ Q X_i → P_Q

noting that although the predicate P_Q(x) is defined to be 0 or 1, it is sufficient in this construction for the concentration of the species representing P_Q(x) to be zero or nonzero. This CRC is output-oblivious.

 We can also compute each g_Q(x⃗) (via Lemma <ref>) using copies of the input molecules. 
 This construction is output-oblivious. To compute P_Q(x)g_Q(x) given the concentration species P_Q as nonzero iff P_Q(x) = 1 as shown above, we simply compute the following (assuming Y_Q is the output of the module computing g_Q(x)):
 
    f(P_Q, Y_Q) =  
          Y_Q    P_Q ≠ 0 
    
          0     P_Q = 0,

 which is computed by this output-oblivious CRC: 
    Y_Q + P_Q → Y + P_Q.

 
 
 
 
The CRC computing min is output-oblivious, as seen in the introduction.
The CRC computing the sum of its inputs is output-oblivious (e.g., X_1 → Y, X_2 → Y computes X_1 + X_2).
Since each CRC shown is output-oblivious and thus composable, we can compose the modules described to construct a CRC stably computing min_S ⊆ N[g_S(x⃗) + ∑_K ⊈S P_K(x⃗) g_K(x⃗)], which is equal to f(x⃗) by Lemma <ref>.



Given any superadditive positive-continuous piecewise rational linear function f: ℝ_≥0^n →ℝ_≥0, there exists a composable CRC with reactions with at most two reactants and at most two products which stably computes f.

To deduce this corollary, note that the reactions with more than two reactants and/or products are used to compute the following functions: computation of a rational linear function, copying inputs, min, and predicate computation.
We can decompose such reactions into a set of bimolecular reactions. 
For example, a reaction X_1 + … + X_n → Y_1 +  … + Y_n can be decomposed into the reactions
X_1 + X_2 → X_12, X_12 + X_3 → X_123, …, X_123… n-1 + X_n → Y_12… n-1 + Y_n, Y_12… n-1→ Y_12… n-2 + Y_n-1, …, Y_12→ Y_1 + Y_2.
We can verify that each affected module stably computes correctly with these expanded systems of reactions, and remains composable.














§ EXAMPLE


In this section, we demonstrate the construction presented in the previous section through an example.
Consider the function shown in Equation <ref> in Section <ref>.
As shown in that section, the function is superadditive, positive-continuous, and piecewise rational linear.
Thus, we can apply our construction to generate a composable CRN stably computing this function. Note that while this CRN is generated from our methodology, we have removed irrelevant species and reactions.
   
2
Making copies of input:


    X_1    → X_1” + X_1”'
    
       X_2    → X_2” + X_2”'
    
       X_3    → X_3'

Using X_3' to make P_3, which catalyzes reactions for the domain X_3 > 0: 

    X_3'    → P_{3}

Computing the sum in Y_{3}:

    X_1”   → Y_{3}
    
       X_2”   → Y_{3}

Computing the min in Y_∅:

    X_1”' + X_2”'    → Y_∅

Making a copy of Y_{3} for use in increasing Y_∅':

    Y_{3}   → Y_{3}' + Y_{3}, ∅

Increase Y_∅' so that it will not be the min when x_3 is present:

    Y_{3}, ∅ + P_{3}   → Y_∅' + P_{3}

Rename Y_∅ to Y_∅' so that it will be summed with the term created by the previous reaction:

    Y_∅   → Y_∅'
    
    
       Y_∅' + Y_{3}'    → Y







§ FUNCTIONS COMPUTABLE BY COMPOSABLE CRNS WITH INITIAL CONTEXT

So far, our CRCs restrict the concentrations of non-input species in the initial state to be zero.
One may consider some (constant) initial concentration of non-input species, called initial context, and how that may affect computation.


A CRC with initial context is denoted ^I,i⃗ = (Λ,R,Σ,Y,I, i⃗)
with Λ, R, Σ, and Y defined as in Definition <ref>, and the initial context species I ⊂Λ∖ ( Σ∪ Y ) and initial context concentrations i⃗∈^I_≥0.
^I,i⃗ stably computes f : ^n → if, for all ∈^n and all  such that + i⃗→, there exists an output stable state  such that → and (Y) = f().

We will show that such CRCs can compute exactly what we could already compute when we set one of our inputs to one. Note that this results in a larger class of computable functions than superadditive positive-continuous piecewise rational linear. For example, we can now compute the non superadditive function f(x_1) = 1. To formalize this idea, we will also need the following definition.






Consider a piecewise affine (linear with an offset) function f: ^n → defined by g_i on domain D_i.
The linear extension f': ^n ×ℝ_>0→ of f is constructed as follows: for all x⃗∈^n, γ∈ℝ_>0 with x⃗/γ∈ D_i, let f'(x⃗,γ) = g_i(x⃗) - g_i(0⃗) + g_i(0⃗) γ.






For example, given the piecewise affine function 
    f(x_1) =  1     x_1 > 1 
     x_1     x_1 ≤ 1
 its linear extension would be 
    f'(x_1,γ) = γ    x_1 > γ
     x_1     x_1 ≤γ.
 Observe that f' is just the minimum of x_1 and γ.
Intuitively the linear extension replaces all constant terms with a linear term in the new variable γ.

We show that initial context for composable CRCs allows only functions whose linear extensions are superadditive, positive-continuous, piecewise rational linear. Intuitively, this gives us the same class of functions that we would have without initial context when we set one of our inputs to the constant value one.




As defined, we allow an arbitrary number of initial context species with differing initial concentrations, but we will focus on the single species case with an initial concentration of one.
This is well motivated: we show one initial context species with concentration one is equivalent in stable computing power to having any number of species with nonnegative rational initial concentrations.


Given a CRC with initial context ^I,i⃗ with i⃗∈ℚ^I_≥0 (rational initial concentrations) which stably computes f, there exists a CRC ^I',i⃗'⃗
with I' = {S'} and i⃗'⃗(S') = 1 (one initial species with concentration one) which stably computes f. 






Let q_i = a_i/b_i for a_i,b_i ∈ℤ_≥ 0 be the initial (rational) concentrations of the initial context in ^I,i⃗.
Observe that the CRN:

    S' → S_1' + S_2' + … + S_k'

can be used to produce k species with concentrations equal to S''s initial concentration.
Then the CRN:

    bS_i' → aS_i

for each i produces a concentration of q_i for species S_i. 

While this schema cannot be used to generate initial context with irrational concentrations, continued fractions can be used to approximate irrational numbers as rational numbers with arbitrarily small error. Thus our restriction to one initial species with a initial concentration one is reasonable to consider for stable computation in this model.
To characterize the functions stably computable with initial context, we first prove some lemmas.
Recall Post(c⃗) is the set of states reachable from c⃗, i.e., {d⃗ |  c⃗→d⃗}.


Given a CRN {Λ, R}, for any γ∈ℝ_≥ 0 and c⃗∈^Λ, Post(γc⃗) = γPost(c⃗).


If c⃗→d⃗, then for any k∈ℝ_≥ 0, kc⃗→ kd⃗. This can be verified by taking the straight line segments to get from c⃗ to d⃗ and scaling them by k.

When γ=0, this lemma is trivial, as the only state reachable from the zero state is the zero state and zero times any state yields the zero state. Thus we only need to consider the case where γ > 0.

Let v⃗∈Post(γc⃗). By the definition of Post this implies that γc⃗→v⃗. This implies that c⃗→1/γv⃗. Therefore 1/γv⃗∈Post(c⃗) and v⃗∈γPost(c⃗). 

Let v⃗∈γPost(c⃗). By the definition of Post this implies that c⃗→1/γv⃗. This implies that γc⃗→v⃗. Which implies that v⃗∈Post(γc⃗).


Let ^I,i⃗ with I = {S_1} and i⃗(S_1) = q_1 stably compute f. Then ^I,γi⃗ for γ∈_>0 stably computes some function g.





Consider running  with an initial concentration of q_1' = γ q_1 for the species S_1. Observe that the initial state x⃗'⃗ = 1/γx⃗, where x⃗ is a state that will stably compute f under the definition of . By Lemma <ref>, we know that the set of reachable states from x⃗'⃗ is equal to the set of reachable states from x⃗ scaled by γ. Thus consider some state c⃗'⃗ reached from x⃗'⃗. Observe that there exists a state c⃗ reachable from x⃗ such that c⃗'⃗ = γc⃗. Consider some output stable state o⃗ reachable from c⃗. Observe that by Lemma <ref> the state o⃗'⃗ = γo⃗ must be reachable from c⃗'⃗. Likewise by Lemma <ref> we know that o⃗'⃗ must be an output stable state. Thus, we know that  must stably compute some function regardless of the initial value for S_1.


We thus know that scaling the value of the initial context retains the fact that  stably computes a function in the region where that species has a positive concentration.
We can then view a single species of initial context as an additional input to  and claim that this CRN stably computes some function, using lemmas from the case of no initial context to prove properties of that function.

Let 𝒞^I,i⃗ be output-oblivious with I = {S_1} and i⃗(S_1) = 1. 𝒞^I,i⃗ stably computes f : ℝ_≥ 0^n →ℝ_≥ 0 only if f is a piecewise affine function whose linear extension is superadditive positive-continuous piecewise rational linear.


Intuitively, we treat the initial context as an input species. 
For 𝒞^I,i⃗ = (Λ, R,Σ, Y, I, i⃗), let 𝒞' = (Λ, R, Σ∪ I, Y).
By assumption, 𝒞^I,i⃗ is output-oblivious, so 𝒞' is output-oblivious.
Further, by Lemma <ref>, 𝒞' stably computes some function for all (positive) initial concentrations of S_1.
So 𝒞' must stably compute a superadditive, positive-continuous, piecewise rational linear function f' on all domains with positive (nonzero) initial concentration of S_1. Since 𝒞^I,i⃗ and 𝒞' share the same CRN, output species, and input species other than S_i, we know that for any x⃗∈ℝ_≥ 0^Σ, f'(x⃗,1) = f(x⃗).
All that remains is to show that f' is the linear extension of f. Pick any γ∈ℝ_>0, x⃗∈^n  such that x⃗ / γ∈ D_i. Let g_i be the affine function used by f on D_i. 
Since f'(x⃗/γ,1) = f(x⃗/γ) and the domains of f' are cones by Lemma <ref>, f'(x⃗,γ) = γ g_i(x⃗ / γ). Expanding out the right hand side gives us that γ f'(x⃗/γ,1) = g_i(x⃗) - g_i(0⃗) + γ g_i(0⃗). By Definition <ref> f' is the linear extension of f.










Let ^I,i⃗ be output-oblivious with i⃗∈ Q^I_≥0. Then ^I,i⃗ stably computes f : ^n → if and only if f is a rational affine function whose linear extension is super-additive, positive-continuous, piecewise rational linear.


 
Lemma <ref> gives us the necessary conditions for composable CRCs with initial context. Now all we need to do is show that there is a composable CRC for any rational affine function f whose linear extension is superadditive, positive-continuous, piecewise rational linear. Since the linear extension of f is superadditive, positive-continuous, piecewise rational linear, Lemma <ref> gives a CRC that computes this function. Observe that changing the input in this CRC representing the linear extension (γ) to the initial context species will give us a new CRC that computes f.




§ DISCUSSION



Instead of continuous concentrations of species, one may consider discrete counts.
This changes which functions are stably computed by CRNs.
Without the composability constraint, <cit.> shows in the discrete model that a function f : ^n → is stably computable by a direct CRN if and only if it is semilinear; i.e., its graph {(x,y) ∈^n × | f(x) = y } is a semilinear subset of ^n ×.
The proof that composably computable functions must be superadditive (Lemma <ref>) holds for the discrete model as well.
Additionally, there exist functions which are superadditive and semilinear but are not computable in the discrete model by a composable CRN.
For example (the proof is omitted):

    f(x_1, x_2) =  
           x_1 - 1     x_1 > x_2 
    
           x_1     x_1 ≤ x_2,

so the class of composably computable functions is slightly more restricted.

Three recent works characterized output-oblivious computation in the discrete model.
With initial context, the characterization for functions on two inputs (f : ^2 →) was given by <cit.>,
and subsequently extended to arbitrarily many inputs <cit.>.
More recently, <cit.> showed that without initial context, computable functions must be superadditive in addition to the constraints presented in <cit.>.

















Our negative and positive results are proven with respect to stable computation, which formalizes our intuitive notion of rate-independent computation.
It is possible to strengthen our positive results to further show that our CRNs converge (as time t →∞) to the correct output from any reachable state under mass-action kinetics (proof omitted). 
It is interesting to characterize the exact class of rate laws that guarantee similar convergence.

Apart from the dual-rail convention discussed in the introduction, other input/output conventions for computation by CRNs have been studied.
For example, <cit.> considers fractional encoding in the context of rate-dependent computation.

As shown by dual-rail, different input and output conventions can affect the class of functions stably computable by CRNs. While using any superadditive positive continuous piecewise rational linear output convention gives us no extra computational power—since the construction in this paper shows how to compute such an output convention directly—it is unclear how these conventions change the power of rate-independent CRNs in general.



Finally we can ask what insights the study of composition of rate-independent modules gives for the more general case of rate-dependent computation.
Is there a similar tradeoff between ease of composition and expressiveness for other classes of CRNs?




























plain






§ APPENDIX


Most definitions and lemmas in this section work towards proving Lemma <ref>.
We also include a proof of Lemma <ref>.

A polyhedron is a subset of ^n of the form {x⃗ | Ax⃗≤b⃗} for some m × n matrix A and some vector b⃗∈^m.



A convex polytope is the convex hull of a finite set of points in ^n.



A polyhedral cone is a set of the form {x⃗ = λ_1x⃗_1 + … + λ_nx⃗_n | λ_1,…,λ_n ≥ 0} for some finite set of points x_1, …, x_n in ^n

The following lemma comes from a previous work. Note that this sum is the Minkowski sum.
 [Proven in <cit.>]
A subset P ⊆^n is a polyhedron if and only if P = Q + C for some convex polytope Q and some polyhedral cone C.



For a given state c⃗ of a CRN 𝒞, the set of states {d⃗ | c⃗→^1 d⃗} that are straight-line reachable from c⃗ is a polyhedron.



If m = |R| is the number of reactions in 𝒞 and n = |Λ| is the number of species in 𝒞, then the stoichiometry matrix can be thought of as a linear transformation from the reaction space ^m to the species space ^n. Let ℛ_c⃗ be the set of basis vectors corresponding to reactions applicable at c⃗. Then since M is a linear transformation, it sends the polyhedral cone defined by ℛ_c⃗ to a polyhedral cone C' in ^n. By Lemma <ref>, the set c⃗ + C' is a polyhedron in ^n, and since the set of states that are straight-line reachable from c⃗ is the intersection of this polyhedron with the set of all vectors with nonnegative components, it is also a polyhedron.



The set of possible species produced from a state c⃗ is 

    𝒫(c⃗) = ⋃_d ∈Post(c⃗) [d⃗]




For a CRN the set of reachable states is closed under convex combination.


Consider some state c⃗ and states S=α⃗_⃗1⃗, …, α⃗_⃗k⃗ reachable from c⃗. Let d⃗ = ∑_i=0^k a_i α_i where ∀ i  a_i>0 and ∑_i=0^k a_i = 1. By lemma <ref> we know that ∀ i  a_iα⃗_⃗i⃗ is reachable from a_i c⃗. Since ∑_i=0^k a_i c⃗ = c⃗, we know that c⃗→∑_i=0^k a_i α⃗_⃗i⃗, which is equal to d⃗.



Given a state c⃗, there is a state d⃗ reachable from c⃗ such that 𝒫(c⃗) = [d⃗]. For such a state, 𝒫(d⃗) = 𝒫(c⃗).



If c⃗ is the zero vector, observe that 𝒫(c⃗) = [c⃗], so setting d⃗=c⃗ makes this hold. Otherwise, for each species S ∈𝒫(c⃗), there is some state d⃗_S reachable from c⃗ with S ∈ [d⃗_S]. Then the state d⃗ = 1/|𝒫(c⃗)|∑_S ∈𝒫(c⃗)d⃗_S is reachable from c⃗ by lemma <ref>. Since d⃗ contains a positive contribution from each d⃗_S, 𝒫(c⃗) = [d⃗]. Since c⃗→d⃗ we know that 𝒫(d⃗) ⊆𝒫(c⃗). Since 𝒫(c⃗) = [d⃗] and [d⃗] ⊆𝒫(d⃗) we know that 𝒫(c⃗) ⊆𝒫(d⃗). Thus we can conclude that 𝒫(d⃗) = 𝒫(c⃗).





If c⃗ is a state such that 𝒫(c⃗) = [c⃗], then any state d⃗ that is reachable from c⃗ is straight-line reachable from c⃗.



Since 𝒫(c⃗) = [c⃗], by the definition of 𝒫(c⃗) we know that the set of applicable reactions from c⃗ is a superset of those applicable from any state reachable from c⃗. Thus we can take the sum of all the straight-line segments used to reach d⃗ from c⃗ and apply them all in a single straight-line segment to get c⃗→^1 d⃗.





Lemma <ref>. 
For any state c⃗ and any species S, if the amount of S present in any state reachable from c⃗ is bounded above, there is a state d⃗ reachable from c⃗ such that for any state a⃗ reachable from d⃗, we know that a⃗(S) ≤d⃗(S).



By Lemma <ref>, there is some c⃗_1 reachable from c⃗ such that [c⃗_1] = 𝒫(c⃗). By Lemma <ref>, we know that the states that are straight-line reachable from c⃗_1 are a polyhedron P. The linear map ^n → sending x⃗↦x⃗(S) maps P to some polyhedral subset of —in particular this is a closed subset. Since we assume that the image of this map is bounded above, we know that this subset attains its maximum M, so there is some d⃗∈ P with d⃗(S) = M. Any state a⃗ that is reachable from d⃗ is also reachable from c⃗_1, so by Lemma <ref> it is contained in P. As a result, a⃗(S) ≤ M = d⃗(S).


Lemma <ref>. 
 Suppose we are given a continuous piecewise rational linear function f: ℝ_>0^n →ℝ_≥ 0. Then we can choose domains for f which are cones which contain an open ball of non-zero radius.





Since f is piecewise rational linear, we can pick a finite set of domains 𝒟 = D_i_i = 1^N for f, such that f|_D_i = g_i|_D_i, where g_i is a rational linear function. Fix a domain D_k, and consider any point x⃗∈ D_k. Since the open ray ℓ_x⃗ from the origin passing through x⃗ is contained in ℝ_>0^n, it is covered by the domains in 𝒟. If we write any point y⃗∈ℓ_x⃗ in the form t ·x⃗, then, for each i, the restriction of g_i to D_i ∩ℓ_x⃗ is of the form g_i(t ·x⃗) = α_i t for some α_i ∈ℝ. Since x⃗∈ D_k ∩ℓ_x⃗, we know that f(1 ·x⃗) = α_k · 1 = α_k

Now suppose that for some s ∈ℝ_> 0 we know that f(s ·x⃗) ≠α_k s. First consider the case where s > 1. Then define the set A = t ∈ [1,s]  | f(t·x⃗) = α_k t and define the set B = t ∈ [1,s]  | f(t·x⃗) ≠α_k t. We know that A is non-empty since 1 ∈ A, so sup A exists - call it t'. From the standard properties of the supremum, we know that there exists a sequence of points t_j_j = 1^∞ such that t_j ∈ A for all j and lim_j →∞ t_j = t'. As a result, from the continuity of f, we see that:


    f(t' ·x⃗) = lim_j →∞ f(t_j ·x⃗) = lim_j →∞α_k t_j = α_k t'


So t' ∈ A. However, by assumption, s ∈ B, so that t' < s. Since t' is an upper bound on A, it must then be the case that (t', s] ⊆ B, so that there exists a sequence of points s_j_j = 1^∞ such that s_j ∈ B for all j and lim_j →∞ s_j = t'. Since there are only finitely many domains in 𝒟, but infinitely many s_j, by the pigeonhole principle infinitely many of the s_j must be contained in a single domain D_k'. Now write the subsequence of points contained in D_k' as s_j'_j' = 1^∞. We still know that lim_j' →∞ s_j' = t', so by the continuity of f and the fact that s_j'∈ D_k', we see that:


    α_kt' = f(t'·x⃗) = lim_j' →∞ f(s_j'·x⃗) = lim_j' →∞α_k's_j' = α_k't'


Since t' > 0, this implies that α_k' = α_k, so that f(s_j'·x⃗) = α_k s_j'. However, this contradicts the fact that we were able to choose s_j'∈ B. As a result, our assumption, that there is some s > 1 such that f(s·x⃗) ≠α_k s, must be false. A similar argument, using the infimum instead of the supremum, shows that there can be no s < 1 such that f(s·x⃗) ≠α_k s. As a result, for every point t ∈ℓ_x⃗, we know f(t ·x⃗) = α_k t. In other words, f|_ℓ_x⃗ = g_k|_ℓ_x⃗, so we can replace D_k with D_k ∪ℓ_x⃗ without issue. Doing this for every x⃗∈ D_k, we can replace D_k with a cone. By enlarging every domain in 𝒟 in this way, we can choose domains for f which are cones.

Since f is continuous, we can replace each D_i ∈𝒟 by its closure, which is again a cone. Suppose that for any D_i ∈𝒟, there is a point x⃗∈ D_i not in the interior of D_i. Then x⃗ is in the closure of the complement of D_i, so there exists a sequence x⃗_k_k = 1^∞ of points in the complement of D_i such that lim_k →∞x⃗_k = x⃗. Since the complement of D_i is covered by the D_j ∈𝒟, where j ≠ i, we know that each x⃗_k lies in one of the D_j. Since there are only finitely many D_j but infinitely many x⃗_k, we know that infinitely many x⃗_k must lie in at least one of the D_j. As a result, x⃗ is in the closure of this D_j, and since D_j is closed, we see that x⃗∈ D_j. Because of this, if D_i has no interior points, then it is completely contained in the other D_j, so we can remove it from the set of domains. After doing this for every D_i which contains no interior points, we can ensure that the domains we have chosen for f all contain an open ball of non-zero radius.































