Programming discrete distributions with chemical reaction networks

Abstract

We explore the range of probabilistic behaviours that can be engineered with Chemical Reaction Networks (CRNs). We give methods to “program” CRNs so that their steady state is chosen from some desired target distribution that has finite support in ${\mathbb {N}}^m$, with $m \ge 1$. Moreover, any distribution with countable infinite support can be approximated with arbitrarily small error under the $L^1$ norm. We also give optimized schemes for special distributions, including the uniform distribution. Finally, we formulate a calculus to compute on distributions that is complete for finite support distributions, and can be compiled to a restricted class of CRNs that at steady state realize those distributions.

1 Introduction

Individual cells and viruses operate in a noisy environment and molecular interactions are inherently stochastic. How cells can tolerate and take advantage of noise (stochastic fluctuations) is a question of primary importance. It has been shown that noise has a functional role in cells (Eldar and Elowitz 2010); indeed, some critical functions depend on the stochastic fluctuations of molecular populations and would be impossible in a deterministic setting. For instance, noise is fundamental for probabilistic differentiation of strategies in organisms, and is a key factor for evolution and adaptation (Arkin et al. 1998). In Escherichia coli, randomly and independently of external inputs, a small sub-population of cells enters a non-growing state in which they can elude the action of antibiotics that can only kill actively growing bacterial cells. Thus, when a population of E. coli cells is treated with antibiotics, the persisted cells survive by virtue of their quiescence before resuming growth (Losick and Desplan 2008). This is an example in which molecular systems compute by producing a distribution. In other cases cells need to shape noise and compute on distributions instead of simply mean values. For example, in Schmiedel et al. (2015) the authors show, both mathematically and experimentally, that microRNA confers precision on the protein expression: it shapes the noise of genes in a way that decreases the intrinsic noise in protein expression, maintaining its expected value almost constant. Thus, although fundamentally important, the mechanisms used by cells to compute in a stochastic environment are not well understood.

Chemical Reaction Networks (CRNs) with mass action kinetics are a well studied formalism for modelling biochemical systems, more recently also used as a formal programming language (Chen et al. 2013). It has been shown that any CRN can be physically implemented by a corresponding DNA strand displacement circuit in a well-mixed solution (Soloveichik et al. 2010). DNA-based circuits thus have the potential to operate inside cells and control their activity. Winfree and Qian have also shown that CRNs can be implemented on the surface of a DNA nanostructure (Qian and Winfree 2014), enabling localized computation and engineering biochemical systems where the molecular interactions occur between few components. When the number of interacting entities is small, the stochastic fluctuations intrinsic in molecular interactions play a predominant role in the time evolution of the system. As a consequence, “programming” a CRN to provide a particular probabilistic response for a subset of species, for example in response to environmental conditions, is important for engineering complex biochemical nano-devices and randomized algorithms. In this paper, we explore the capacity of CRNs to “exactly program” discrete probability distributions. That is, we give methods such that the steady state distribution of a CRN can be chosen from some desired target distribution. We aim to characterize the probabilistic behaviour that can be obtained, exploring both the capabilities of CRNs for producing distributions and for computing on distributions by composing them.

Contributions We show that at steady state CRNs are able to compute any distribution with finite support in ${\mathbb {N}}^m$, with $m\ge 1$. We propose an algorithm to systematically “program” a CRN so that at steady state it produces any given finite support distribution. Moreover, any distribution with countable infinite support can be approximated with arbitrarily small error under the $L^1$ norm. The resulting network has a number of reactions linear in the dimension of the support of the distribution and the output is produced monotonically allowing composition. Since distributions with large support can result in unwieldy networks, we also give optimised networks for special distributions, including a novel scheme for the uniform distribution. We formulate a calculus that is complete for finite support distributions, which can be compiled to a restricted class of CRNs that at steady state compute those distributions. The resulting CRNs are generally more compact with respect to the ones derived from direct approach. The calculus is equivalent to the baricentric algebra presented in Mardare et al. (2016), and allows for modelling of external influences on the species. Our results are of interest for a variety of scenarios in systems and synthetic biology. For example, they can be used to program a biased stochastic coin or a uniform distribution, thus enabling implementation of randomized algorithms and protocols in CRNs.

Preliminary version of this work appeared as Cardelli et al. (2016a). This paper includes an extended description with illustrative examples and proofs of the results.

Related work It has been shown that CRNs with stochastic semantics are Turing complete, up to an arbitrarily small error (Soloveichik et al. 2008). If we assume error-free computation, their computational power decreases: they can decide the class of the semi-linear predicates (Angluin et al. 2007) and compute semi-linear functions (Chen et al. 2014). A first attempt to model distributions with CRNs can be found in Fett et al. (2007), where the problem of producing a single distribution is studied. However, their circuits are approximated and cannot be composed to compute operations on distributions.

2 Chemical reaction networks

A chemical reaction network (CRN) $(\varLambda ,R)$ is a pair of finite sets, where $\varLambda$ is the set of chemical species, $|\varLambda |$ denotes its size, and R is a set of reactions. A reaction $\tau \in R$ is a triple $\tau =(r_{\tau },p_{\tau },k_{\tau })$, where $r_{\tau } \in {\mathbb {N}}^{|\varLambda |}$ is the source complex, $p_{\tau } \in {\mathbb {N}}^{|\varLambda |}$ is the product complex and $k_{\tau } \in {\mathbb {R}}_{>0}$ is the coefficient associated to the rate of the reaction, where we assume $k_{\tau }=1$ if not specified; $r_{\tau }$ and $p_{\tau }$ represent the stoichiometry of reactants and products. Given a reaction $\tau _1=( [1,0,1],[0,2,0],k_1 )$ we often refer to it as $\tau _1 : \lambda _1 + \lambda _3 \, \rightarrow ^{k_1} \, 2\lambda _2$. The net change (or state change) associated to $\tau$ is defined by $\upsilon _{\tau }=p_{\tau } - r_{\tau }$.

We assume that the system is well stirred, that is, the probability of the next reaction occurring between two molecules is independent of the location of those molecules, at fixed volume V and temperature. Under these assumptions a configuration or state of the system $x \in {\mathbb {N}}^{|\varLambda |}$ is given by the number of molecules of each species.

A chemical reaction system (CRS) $C=(\varLambda ,R,x_0)$ is a tuple where $(\varLambda ,R)$ is a CRN and $x_0 \in {\mathbb {N}}^{|\varLambda |}$ represents its initial condition.

2.1 Stochastic semantics

The stochastic semantics of a CRS is given in terms of a continuous time Markov chain (CTMC). Here, we introduce the semantics according to the representation of Markov processes proposed by Ethier and Kurtz (2009, Theorem 4.1 Chapter 6). Such representation is equivalent to the classical model described by the Chemical Master Equation, but much more compact. It allows us to represent the CTMC in terms of stochastic equations, which have a similar structure to the deterministic rate equations. We illustrate the semantics with the help of Example 1. Below we present Poisson processes, as they will be used in the semantics and in the paper. A building block of the mathematical models we use in the paper is a counting process. Intuitively, a counting process Y is a process such that Y(t) counts the number of times that a particular phenomenon has been observed by time t.

Definition 1

(Counting process) Y is a counting process if $Y(0)=0$ and Y is constant except for jumps of $+1.$

Definition 2

(Poisson process) A counting process Y is a Poisson process if:

• Number of observations in disjoint time intervals are independent random variables, that is, $Y(t_k)-$ $Y(t_{k-1}),$ $k \in {\mathbb {N}}$, are independent random variables.

• The distribution of $Y(t+\varDelta t)-Y(t)$ is independent of t.

Theorem 1

(Anderson and Kurtz 2015) If Y is a Poisson process, then there exists a constant $\lambda >0$ such that for $t_2>t_1 \in {\mathbb {R}}_{\ge 0}$ and $k \in {\mathbb {N}}$ it holds that

$$\begin{aligned} Prob(Y(t_2)-Y(t_1)=k)=\frac{(\lambda (t_2-t_1))^k}{k!}e^{-\lambda (t_2 - t_1)} \end{aligned}$$

That is, $Y(t_2)-Y(t_1)$ is Poisson distributed with parameter $\lambda (t_2-t_1).$

If $\lambda =1,$ we call Y a unit Poisson process.

Example 1

Consider the CRN described by the following reactions

$$\begin{aligned} \tau _1: \lambda _1 + \lambda _2 \rightarrow ^{k_1} \lambda _1 + \lambda _1;\quad \tau _2: \lambda _1 + \lambda _2 \rightarrow ^{k_2} \lambda _2 + \lambda _2 \end{aligned}$$

and let $X(0)\in {\mathbb {N}}^2$ be the initial condition. Then, the state of the system at time $t\ge 0$ will be given by X(0) plus the number of times that each reaction have fired between [0, t] multiplied by the respective state change vector. That is,

$$\begin{aligned} X(t)=X(0) + \begin{pmatrix} 1 \\ -1 \end{pmatrix} R_{\tau _1}(t) + \begin{pmatrix} -1 \\ 1 \end{pmatrix} R_{\tau _2}(t) \end{aligned}$$

where $R_{\tau _1}(t),R_{\tau _2}(t)$ are counting processes that count the number of times that the particular reaction has fired until time t. We now assume that $R_{\tau }$ are independent, unit Poisson processes that depend on the propensity rate of $\tau$. More precisely, $R_{\tau }(t)=Y_{\tau }(\int _0^t \alpha (X(s))ds),$ where $Y_{\tau }(\int _0^t \alpha (X(s))ds)$ is a unit Poisson process with intensity $\int _0^t \alpha (X(s))ds$. Intuitively, $\int _0^t \alpha (X(s))ds$ gives the time interval in which counting events for the unit Poisson process. Under this modelling assumptions it holds that (Ethier and Kurtz 2009)

$$\begin{aligned}Prob(Y_{\tau }\left( \int _0^{t+\varDelta t} \alpha _{\tau }(X(s))ds\right) -Y_{\tau }\left( \int _0^{t}\alpha _{\tau }(X(s)ds)>0|\forall s \in [0,t,X(s)\right) \approx \alpha _{\tau }(X(t)) \varDelta t . \end{aligned}$$

That is, the probability that a reaction $\tau$ happens in the next $\varDelta t$, at the first order, is given by the propensity rate of $\tau$ at time t multiplied by $\varDelta t$, exactly as in the classical stochastic representation (Van Kampen 1992) of CRNs. At this point, for our model, we can write its stochastic model as

$$\begin{aligned} X(t)=&X(0) \\&+\begin{pmatrix} 1 \\ -1 \end{pmatrix} Y_{\tau _1}\left( k_{\tau _1}\int _0^t X_{\lambda _1}(s)X_{\lambda _2}(s)ds\right) \\&+ \begin{pmatrix} -1 \\ 1 \end{pmatrix} Y_{\tau _2}\left( k_{\tau _2}\int _0^t X_{\lambda _1}(s)X_{\lambda _2}(s)ds\right) . \end{aligned}$$

Theorem 2 below shows that the forward equation associated with the Markov process described in the previous stochastic equation is exactly the Chemical Master Equation (CME).

Definition 3

Given a CRS $C=(\varLambda ,R,x_0),$ we define its stochastic semantics at time t as

$$\begin{aligned} X^C(t)=x_0 + \sum _{\tau \in R} \upsilon _{\tau } Y_{\tau }\left( \int _0^t \alpha _{\tau }(X^C(s)ds)\right) \end{aligned}$$

(1)

where $Y_{\tau }$ are unit Poisson processes, independent of each other.

Theorem 2

(Ethier and Kurtz 2009) Let $C=(\varLambda ,R,x_0)$ be a CRS and $X^C$ be the stochastic process as defined in Eq. (1). Define $Prob(X^C(t)=x|X^C(0)=x_0)=P^C(t)(x)$. Assume that, for each $\tau \in R$ and $t\in {\mathbb {R}}_{\ge 0},$ $X^C(t)< \infty ,$ then

$$\frac{d P^C(t)(x)}{dt} = \sum _{\tau \in R} P^C(t)(x-\upsilon _{\tau }) \alpha _\tau (X^C(t))-P^C(t)(x)\alpha _\tau (X^C(t)).$$

(2)

$P^C(t)(x)$ represents the transient evolution of $X^C$, and can be calculated exactly by solving directly the Chemical Master Equation or by approximation techniques (Cardelli et al. 2016b, c; Bortolussi et al. 2016).

Definition 4

The steady state distribution (or limit distribution) of $X^C$ is defined as $\pi ^C= \lim _{t \rightarrow \infty } P^C(t).$

When clear from the context, we omit the superscript indicating the CRN and simply write $\pi$ instead of $\pi ^C$. $\pi$ calculates the percentage of time, in the long-run, that X spends in each state $x \in S$. If S is finite, then the above limit distribution always exists and is unique (Kwiatkowska et al. 2007). In this paper we focus on discrete distributions, and will sometimes conflate the term distribution with probability mass function, defined next.

Definition 5

Suppose that $M: S \rightarrow {\mathbb {R}}^m$ with $m >0$ is a discrete random variable defined on a countable sample space S. Then the probability mass function (pmf) $f: {\mathbb {R}}^m \rightarrow [0, 1]$ for M is defined as $f(x) = Prob(s \in S \mid M(s) = x).$

For a pmf $\pi : {\mathbb {N}}^{m} \rightarrow [0,1]$ we call $J=\{y\in {\mathbb {N}}^{m}|\pi (y)\ne 0\}$ the support of $\pi$. A pmf is always associated to a discrete random variable whose distribution is described by the pmf. Sometimes, when we refer to a pmf, we imply the associated random variable. Given two pmfs $f_1$ and $f_2$ with values in ${\mathbb {N}}^m$, $m>0$, we define the $L^1$ norm (or distance) between them as $d_1(f_1,f_2)=\sum _{n\in {\mathbb {N}}^m}(|f_1(n)-f_2(n)|)$. Note that, as $f_1,f_2$ are pmfs, then $d_1(f_1,f_2)\le 2$. It is worth stressing that, given the CTMC X, for each $t \in {\mathbb {R}}_{\ge 0}$, X(t) is a random variable defined on a countable state space. As a consequence, its distribution is given by a pmf. Likewise, the limit distribution of a CTMC, if it exists, is a pmf.

Definition 6

Given $C=(\varLambda ,R)$ and $\lambda \in \varLambda$, we define $\pi _{\lambda } (k)=\sum _{ \{x\in S |x(\lambda )=k\}} \pi (x)$ as the probability that for $t \rightarrow \infty$, in $X^C$, there are k molecules of $\lambda$.

$\pi _{\lambda }$ is a pmf representing the steady state distribution of species $\lambda$.

3 On computing finite support distributions with CRNs

We now show that, for a pmf with finite support in ${\mathbb {N}}$, we can always build a CRS such that, at steady state (i.e. for $t \rightarrow \infty$) the random variable representing the molecular population of a given species in the CRN is equal to that distribution. Such result allows us to approximate any distribution with countable infinite support with arbitrarily small error under the $L^1$ norm. The result is then generalised to distributions with domain in ${\mathbb {N}}^m$, with $m \ge 1$. The approximation is exact in case of finite support.

3.1 Programming pmfs

Definition 7

Given $f: {\mathbb {N}} \rightarrow [0,1]$ with finite support $J=(z_1,\ldots ,z_{|J|})$ such that $\sum _{i=1}^{|J|}f(z_i)=1$, we define the CRS $C_{f}=(\varLambda ,R,x_0)$ as follows. $C_{f}$ is composed of 2|J| reactions and $2|J|+2$ species. For any $z_i \in J$ we have two species $\lambda _i,\lambda _{i,i} \in \varLambda$ such that $x_0(\lambda _i)=z_i$ and $x_0(\lambda _{i,i})=0$. Then, we consider a species $\lambda _z \in \varLambda$ such that $x_0(\lambda _z)=1$, and the species $\lambda _{out} \in \varLambda$, which represents the output of the network and such that $x_0(\lambda _{out})=0$. For every $z_i \in J$, R has the following two reactions: $\tau _{i,1}: \lambda _z \rightarrow ^{f(z_i)} \lambda _{i,i}$ and $\tau _{i,2}: \lambda _i + \lambda _{i,i} \rightarrow \lambda _{out} + \lambda _{i,i}$.

Example 2

Consider the probability mass function $f: {\mathbb {N}}\rightarrow [0,1]$ defined as $f(y)= \left\{ \begin{array}{ll} \frac{1}{6}, &{}\quad {\text {if }}\; y=2\\ \frac{1}{3}, &{}\quad {\text {if }}\; y=5\\ \frac{1}{2}, &{}\quad {\text {if }}\; y=10\\ 0, &{}\quad {\text {otherwise}}\\ \end{array} \right.$. Let $\varLambda =\{\lambda _{1}, \lambda _{2}, \lambda _{3},$ $\lambda _z, \lambda _{1,1}, \lambda _{2,2}, \lambda _{3,3}, \lambda _{out} \}$, then we build the CRS $C=(\varLambda ,R,x_0)$ following Definition 7, where R is given by the following set of reactions:

$$\begin{aligned}&\lambda _z \rightarrow ^{\frac{1}{6}} \lambda _{1,1};\quad \lambda _z \rightarrow ^{\frac{1}{3}} \lambda _{2,2};\quad \lambda _z \rightarrow ^{\frac{1}{2}} \lambda _{3,3};\\&\lambda _{1} + \lambda _{1,1} \rightarrow ^{1} \lambda _{1,1} + \lambda _{out};\quad \lambda _{2} + \lambda _{2,2} \rightarrow ^{1} \lambda _{2,2} + \lambda _{out};\\&\lambda _{3} + \lambda _{3,3} \rightarrow ^{1} \lambda _{3,3} + \lambda _{out}. \end{aligned}$$

The initial condition $x_0$ is $x_0 (\lambda _{out})=x_0(\lambda _{1,1})=x_0(\lambda _{2,2})$ $=x_0(\lambda _{3,3})=0;$ $x_0(\lambda _{1})=2; \, x_0(\lambda _{2})=5;$ $x_0(\lambda _{3})=10;$ $x_0(\lambda _z)=1.$ Theorem 3 ensures $\pi _{\lambda _{out}}=f$.

Theorem 3

Given a pmf $f: {\mathbb {N}} \rightarrow [0, 1]$ with finite support J, the CRS $C_{f}$ as defined in Definition  7 is such that $\pi ^{C_f}_{\lambda _{out}}=f$.

Proof

Let $J=( z_1,..,z_{|J|} )$ be the support of f, and |J| its size. Suppose |J| is finite, then the set of reachable states from $x_0$ is finite by construction and the limit distribution of $X^{C_f}$, the induced CTMC, exists. By construction, in the initial state $x_0$ only reactions of type $\tau _{i,1}$ can fire, and the probability that a specific $\tau _{i,1}$ fires first is exactly:

$$\begin{aligned} \frac{\alpha _{\tau _{i,1}}(x_0)}{\sum _{j=1}^{|J|}\alpha _{\tau _{j,1}}(x_0)}&=\frac{f(z_i) \cdot 1 }{\sum _{j=1}^{|J|}f(z_j) \cdot 1 } \\&=\frac{f(z_i) }{\sum _{j=1}^{|J|}f(z_j) }=\frac{f(z_i)}{1}=f(z_i) \end{aligned}$$

Observe that the firing of the first reaction uniquely defines the limit distribution of $X^{C_{f}}$, because $\lambda _z$ is consumed immediately and only reaction $\tau _{i,2}$ can fire, with no race condition, until $\lambda _i$ are consumed. This implies that at steady state $\lambda _{out}$ will be equal to $x_0(\lambda _{i})$, and this happens with probability $f(x_0(\lambda _{i}))$. Since $x_0(\lambda _{i})=z_i$ for $i \in [1,|J|]$, we have $\pi _{\lambda _{out}}^{C_f}=f$. $\square$

Then, we can state the following corollary of Theorem 3.

Corollary 1

Given a pmf $f: {\mathbb {N}} \rightarrow [0, 1]$ with countable support J, we can always find a finite CRS $C_{f}$ such that $\pi ^{C_f}_{\lambda _{out}}=f$ with arbitrarily small error under the $L^1$ norm.

Proof

Let $J=\{z_1,\ldots ,z_{|J|} \}$. Suppose J is (countably) infinite, that is, $|J| \rightarrow \infty$. Then, we can always consider an arbitrarily large but finite number of points in the support, such that the probability mass lost is arbitrarily small, and applying Definition 7 on this finite subset of the support we have the result.

In order to prove the result consider the function $f'$ with support $J'=\{z_1,\ldots ,z_{k} \}$, $k\in {\mathbb {N}}$, such that $f(z_i)=f'(z_i)$, for all $i \in {\mathbb {N}}_{\le k}$. Consider the series $\sum _{i=1}^{\infty }f(n)$. This is an absolute convergent series by definition of pmf. Then, we have that $\lim _{i\rightarrow \infty }f(i)=0$ and, for any $\epsilon > 0$, we can choose some $\kappa _\varepsilon \in {\mathbb {N}}$, such that:

$$\begin{aligned} \forall k>\kappa _\varepsilon \quad |\sum _{i=1}^k f'(i)-\sum _{i=1}^\infty f(i)| < \frac{\epsilon }{2}. \end{aligned}$$

This implies that for $k>\kappa _\varepsilon$ given $f'_k=\sum _{i=1}^k f'(i)$ we have, $d_1(f'_k,f)<\epsilon$. $\square$

The following remark shows that the need for precisely tuning the value of reaction rates in Theorem 3 can be dropped by introducing some auxiliary species.

Remark 1

In practice, tuning the rates of a reaction can be difficult or impossible. However, it is possible to modify the CRS derived using Definition 7 in such a way the probability value is not encoded in the rates, and we just require that all reactions have the same rates. We can do that by using some auxiliary species $\varLambda _c=\{\lambda _{c_1},\lambda _{c_2},\ldots ,\lambda _{c_{|\varLambda _c|}} \}$. Then, the reactions $\tau _{i,1}$ for $i\in [1,J]$ become $\tau _{i,1}: \lambda _z+\lambda _{c_i} \rightarrow ^{k} \lambda _{i,i}$, for $k\ge 0$, instead of $\tau _{i,1}: \lambda _z \rightarrow ^{f(y_i)} \lambda _{i,i}$, as in the original definition. The initial condition of $\lambda _{c_i}$ is $x_0(\lambda _{c_i})=f(y_i)\cdot L$, where $L \in {\mathbb {N}}$ is such that for $j \in [1,|J|]$ and $J=\{ z_1,\ldots ,z_{|J|} \}$ we have that $f(z_j) \cdot L$ is a natural number, assuming all the $f(z_j)$ are rationals.

Remark 2

In biological circuits the probability distribution of a species may depend on some external conditions. For example, the lambda Bacteriofage decides to lyse or not to lyse with a probabilistic distribution based also on environmental conditions (Arkin et al. 1998). Programming similar behaviour is possible by extension of Theorem 3. For instance, suppose, we want to program a switch that with rate $50+Com$ goes to state $O_1$, and with rate 5000 goes to a different state $O_2$, where Com is an external input. To program this logic we can use the following reactions: $\tau _{1,1} : \lambda _z + \lambda _{c_1} \rightarrow ^{k_1} \lambda _{O_1}$ and $\tau _{1,2} : \lambda _z + \lambda _{c_2} \rightarrow ^{k_1} \lambda _{O_2}$, where $\lambda _{O_1}$ and $\lambda _{O_2}$ model the two logic states, initialized at 0. The initial condition $x_0$ is such that $x_0(\lambda _z)=1$, $x_0(\lambda _{c_1})=50$ and $x_0(\lambda _{c_2})=5000$. Then, we add the following reaction $Com \rightarrow ^{k_2} \lambda _{c_1}$. It is easy to show that if $k_2 \gg k_1$ then we have the desired probabilistic behaviour for any initial value of $Com \in {\mathbb {N}}$. This may be of interest also for practical scenarios in synthetic biology, where for instance the behaviour of synthetic bacteria needs to be externally controlled (Anderson et al. 2006); and, if each bacteria is endowed with a similar logic, then, by tuning Com, at the population level, it is possible to control the fraction of bacteria that perform this task.

In the next theorem we generalize to the multidimensional case.

Theorem 4

Given $f: {\mathbb {N}} ^{m} \rightarrow [0, 1]$ with $m\ge 1$ such that $\sum _{i\in {\mathbb {N}}^m}f(i)=1$, then there exists a CRS $C=(\varLambda ,R,x_0)$ such that the joint limit distribution of $(\lambda _{out_1},$ $\lambda _{out_2},\ldots ,\lambda _{out_m}) \in \varLambda$ approximates f with arbitrarily small error under the $L^1$ distance. The approximation is exact if the support of f is finite.

To prove this theorem we can derive a CRS similar to that in the uni-dimensional case. The firing of the first reaction can be used to probabilistically determine the value at steady state of the m output species, using some auxiliary species.

Example 3

Consider the following probability mass function

$$\begin{aligned} f(y_1,y_2)= \left\{ \begin{array}{l l} \frac{1}{6}, &{}\quad {\text {if }}\; y_1=3 \;{\text { and }}\;y_2=1\\ \frac{1}{3}, &{}\quad {\text {if }}\;y_1=3\;{\text { and }}\;y_2=2\\ \frac{1}{2}, &{}\quad {\text {if }}\; y_1=1\;{\text { and }}\;y_2=5\\ 0, &{}\quad {\text {otherwise}}\\ \end{array} \right. \end{aligned}$$

we present the CRS $C=(\varLambda ,R,x_0)$ that according to its stochastic semantics, for $\lambda _{out_1}, \lambda _{out_2} \in \varLambda$ yields the steady-state distribution $\pi _{\lambda _{out_1},\lambda _{out_2}}$, joint limit distribution of $\lambda _{out_1},\lambda _{out_2}$, exactly equal to f. Let $\varLambda =$ $\{ \lambda _z, \lambda _a, \lambda _b,$ $\lambda _c, \lambda _{1,1}, \lambda _{1,2} \lambda _{2,1}, \lambda _{2,2}, \lambda _{3,1},$ $\lambda _{3,2} \lambda _{out_1},\lambda _{out_2} \}$ and R given by the following set of reactions:

$$\begin{aligned}&\tau _1: \lambda _z \rightarrow ^{\frac{1}{6}} \lambda _{a};\quad \tau _2: \lambda _z \rightarrow ^{\frac{1}{3}} \lambda _{b};\quad \tau _3: \lambda _z \rightarrow ^{\frac{1}{2}} \lambda _{c};\\&\tau _4: \lambda _{1,1} + \lambda _{a} \rightarrow ^{1} \lambda _{a} + \lambda _{out_1};\\&\tau _5: \lambda _{1,2} + \lambda _{a} \rightarrow ^{1} \lambda _{a} + \lambda _{out_2};\\&\tau _6: \lambda _{2,1} + \lambda _{b} \rightarrow ^{1} \lambda _{b} + \lambda _{out_1};\\&\tau _7: \lambda _{2,2} + \lambda _{b} \rightarrow ^{1} \lambda _{b} + \lambda _{out_2};\\&\tau _8: \lambda _{3,1} + \lambda _{c} \rightarrow ^{1} \lambda _{c} + \lambda _{out_1};\\&\tau _9: \lambda _{3,2} + \lambda _{c} \rightarrow ^{1} \lambda _{c} + \lambda _{out_2}; \end{aligned}$$

The initial condition $x_0$ is such that:

$$\begin{aligned}&x_0(\lambda _z)=1;\\&x_0(\lambda _{1,1})=3;\, x_0(\lambda _{1,2})=1;\, x_0(\lambda _{2,1})=3;\\&x_0(\lambda _{2,2})=2; \, x_0(\lambda _{3,1})=1; \, x_0(\lambda _{3,2})=5; \end{aligned}$$

and all other species mapped to zero. The set of reachable states from $x_0$ is finite so the limit distribution exists. The firing of the first reaction uniquely determines the steady state solution. $x_0(\lambda _{i,1})$ and $x_0(\lambda _{i,2})$ for $i \in [1,3]$ are exactly the value of $\lambda _{out_1}$ and $\lambda _{out_2}$ at steady state if the first reaction to fire is $\tau _i$; this happens with probability $f(x_0(\lambda _{i,1}),x_0(\lambda _{i,2}))$. Therefore, we have that, at steady state, the joint distribution of $\lambda _{out_1}$ and $\lambda _{out_2}$ equals f.

3.2 Special distributions

For a given pmf the number of reactions of the CRS derived from Definition 7 is linear in the dimension of its support. As a consequence, if the support is large then the CRSs derived using Theorems 3 and 4 can be unwieldy. In the following we show three optimised CRSs to calculate the Poisson, binomial and uniform distributions. These CRNs are compact and applicable in many practical scenarios. However, using Definition 7 the output is always produced monotonically. In the circuits below this does not happen, but, on the other hand, the gain in compactness is substantial. The first two circuits have been derived from the literature, while the CRN for the uniform distribution is new.

3.2.1 Poisson distribution

The main result of Anderson et al. (2010) guarantees that all the CRNs that respect some conditions (weakly reversible, deficiency zero and irreducible state space, see Anderson et al. 2010) have a distribution given by the product of Poisson distributions. As a particular case, we consider the following CRS composed of only one species $\lambda$ and the following two reactions $\tau _1 : \emptyset \rightarrow ^{k_1} \lambda ; \, \tau _2: \lambda \rightarrow ^{k_2} \emptyset .$ Then, at steady state, $\lambda$ has a Poisson distribution with expected value $\frac{k_1}{k_2}$.

3.2.2 Binomial distribution

We consider the network introduced in Anderson et al. (2010). The CRS is composed of two species, $\lambda _1$ and $\lambda _2$, with initial condition $x_0$ such that $x_0(\lambda _1)+x_0(\lambda _2)=K$ and the following set of reactions: $\tau _1:\lambda _1 \rightarrow ^{k_1} \lambda _2; \tau _2:\lambda _2 \rightarrow ^{k_2} \lambda _1.$ As shown in Anderson et al. (2010), $\lambda _1$ and $\lambda _2$ at steady state have a binomial distribution such that: $\pi _{\lambda _1}(y)=(\frac{K}{y}){c_1}^{y}(1-c_1)^{K-y}$ and $\pi _{\lambda _2}(y)=(\frac{K}{y}){c_2}^{y}(1-c_2)^{K-y}$ .

3.2.3 Uniform distribution

The following CRS computes the uniform distribution over the sum of the initial number of molecules in the system, independently of the initial value of each species. It has species $\lambda _1$ and $\lambda _2$ and reactions:

$$\begin{aligned}&\tau _1:\lambda _{1} \rightarrow ^{k} \lambda _{2}; \quad \tau _2:\lambda _{2} \rightarrow ^{k} \lambda _{1};\\&\tau _3:\lambda _{1} + \lambda _{2} \rightarrow ^{k} \lambda _{1} + \lambda _{1};\quad \tau _4:\lambda _{1} + \lambda _{2} \rightarrow ^{k} \lambda _{2} + \lambda _{2} \end{aligned}$$

For $k >0$, $\tau _1$ and $\tau _2$ implement the binomial distribution. These are combined with $\tau _3$ and $\tau _4$, which implement a Direct Competition (DC) system (Cardelli and Csikász-Nagy 2012). DC has a bimodal limit distribution in 0 and in K, where $x_0(\lambda _1)+x_0(\lambda _2)=K$, with $x_0$ initial condition. This network, surprisingly, according to the next theorem, at steady state produces a distribution which varies uniformly between 0 and K.

Theorem 5

Let $x_0(\lambda _1)+x_0(\lambda _2)=K \in {\mathbb {N}}$. Then, the CRS described above has the following steady state distribution for $\lambda _1$ and $\lambda _2$:

$$\begin{aligned} \pi _{\lambda _1}(y)=\pi _{\lambda _2}(y)=\left\{ \begin{array}{l l} \frac{1}{K+1}, &{}\quad {\text {if }}\;y \in [0,K]\\ 0, &{}{\text {otherwise}}\\ \end{array} \right. . \end{aligned}$$

Proof

We consider a general initial condition $x_0$ such that $x_0(\lambda _1)=K-M$ and $x_0(\lambda _2)=M$ for $0 \le M \le K$ and $K,M \in {\mathbb {N}}$. Because any reaction has exactly 2 reagents and 2 products, we have the invariant that for any configuration x reachable from $x_0$ it holds that $x(\lambda _1)+x(\lambda _2)=K$. Figure 1 plots the CTMC semantics of the system.

Fig. 1

The figure shows the CTMC induced by the CRS implementing the uniform distribution for initial condition $x_0$ such that $x_0(\lambda _1)+x_0(\lambda _2)=K$

For any fixed K the set of reachable states from any initial condition in the induced CTMC is finite (exactly K states are reachable from any initial condition) and irreducible. Therefore, the steady state solution exists, is unique and independent of the initial conditions. To find this limit distribution we can calculate Q, the infinitesimal generator of the CTMC, and then solve the linear equations system $\pi Q=0$, with the constraint that $\sum _{i \in [0,K]}\pi _{i}=1$, where $\pi _{i}$ is the ith component of the vector $\pi$, as shown in Kwiatkowska et al. (2007). Because the CTMC we are considering is irreducible, this is equivalent to solving the balance equations with the same constraint. The resulting $\pi$ is the steady state distribution of the system.

We consider 3 cases, where $(K-j,j)$ for $j \in [0,K]$ represents the state of the system in terms of molecules of $\lambda _1$ and $\lambda _2$.

• Case $j=0$. For the state (K, 0), whose limit distribution is defined as $\pi (K,0),$ we have the following balance equation:

  $$\begin{aligned}&-\,\pi (K,0) Kk+\,\pi (K-1,1)[(K-1)k+k]=0 \implies \\&\pi (K,0)=\pi (K-1,1). \end{aligned}$$

• Case $j \in [1,K-1]$. In Fig. 1 we see that the states and the rates follow a precise pattern: every state is directly connected with only two states and for any transition the rates depend on two reactions, therefore we can consider the balance equations for a general state $(K-j,j)$ for $j\in [1,K-1]$ (for the sake of a lighter notation instead of $\pi (K-j,j)$ we write $\pi ^{j}$):

  $$\begin{aligned}&\pi ^{j-1}[ K +1 -j+(K+1-j)(j-1)] - \pi ^{j}[2(K-j)j+j+K-j] +\,\pi ^{j+1}[j+1+(K-j-1)(j+1)] =0 \\&\quad \quad \quad \quad \quad \quad \quad \implies \\&\pi ^{j-1}[ Kj-j^2+j] -\,\pi ^{j}[2Kj-2j^2 + K] +\,\pi ^{j+1}[Kj+K-j^2-j] = 0 \end{aligned}$$

  It is easy to verify that if $\pi ^{j-1}= \pi ^j = \pi ^{j+1}$ then the equation is proved.

• Case $j=K$. The case for the state (0, K) is similar to the case (K, 0).

We have shown that each reachable state has equal probability at steady state for any possible initial condition. Therefore, because $\sum _{i=0}^K \pi ^i = 1$ and $\pi _{\lambda _i}(y)=$ $\sum _{ x_{j} \in S|x_j(\lambda _i)=y} \pi ^j$ for $y \ge 0$, we have that for both $\lambda _1$ and $\lambda _2$

$$\begin{aligned} \pi _{\lambda _1}(y)=\pi _{\lambda _2}(y)=\left\{ \begin{array}{l l} \frac{1}{K+1}, &{}\quad {\text {if }}\;y \in [0,K]\\ 0, &{}{\text {otherwise}}\\ \end{array} \right. \end{aligned}$$

$\square$

4 Calculus of limit distributions of CRNs

In the previous section we have shown that CRNs are able to program any pmf on ${\mathbb {N}}$. We now define a calculus to compose and compute on pmfs. We show it is complete with respect to finite support pmfs on ${\mathbb {N}}$. The calculus we present is a left-invariant baricentric algebra (Mardare et al. 2016). Then, we define a translation of this calculus into a restricted class of CRNs. We prove the soundness of such a translation, which thus yields an abstract calculus of limit distributions of CRNs. For simplicity, in what follows we consider only pmfs with support in ${\mathbb {N}}$, but the results can be generalised to the multi-dimensional case.

Definition 8

(Syntax) The syntax of formulae of our calculus is given by

$$\begin{aligned}&P:=\,(P+P) \, | \, min(P,P) \, | \, k \cdot P \, | \, (P)_D:P \, | \, one \, | \, zero\\&D := p \,|\, p \cdot c_i + D \, \end{aligned}$$

where $k \in {\mathbb {Q}}_{\ge 0}$, $p \in {\mathbb {Q}}_{[0,1]}$ are rational and $V=\{c_1,\ldots ,$ $c_n \}$ is a set of variables with values in ${\mathbb {N}}$.

A formula P denotes a pmf that can be obtained as a sum, minimum, multiplication by a rational, or convex combination of pmfs one and zero. Given a formula P, variables $V=\{c_1,\ldots ,c_n \}$, called environmental inputs, model the influence of external factors on the probability distributions of the system. V(P) represents the variables in P. An environment $E: V \rightarrow {\mathbb {Q}}_{[0,1]}$ is a partial function which maps each input $c_i$ to its valuation normalized to [0, 1]. Given a formula P and an environment E, where $V(P)\subseteq dom(E)$, with dom(E) domain of E, we define its semantics, $[\![P]\!]_E$, as a pmf (the empty environment is denoted as $\emptyset$). D expresses a summation of valuations of inputs $c_i$ weighted by rational probabilities p, which evaluates to a rational $[\![D]\!]_E$ for a given environment. We require that, for any D, the sum of p coefficients in D is in [0, 1]. This ensures that $0\le [\![D]\!]_E \le 1$. The semantics is defined inductively as follows, where the operations on pmfs are defined in Sect. 4.1.

Definition 9

(Semantics) Given formulae P,  $P_1,$ $P_2$ and an environment E, such that $V(P)\cup V(P_1)\cup V(P_2)$ $\subseteq dom(E)$, we define

$$\begin{aligned}&[\![one]\!]_E=\pi _{one}\quad [\![ zero]\!]_E=\pi _{zero}\\&[\![P_1+P_2]\!]_E=[\![P_1]\!]_E + [\![P_2]\!]_E \\&[\! [min(P_1,P_2)]\!]_E=min([\![P_1]\!]_E,[\![P_2]\!]_E) \\&[\![k\cdot P]\!]_E=\frac{k_1\cdot ([\![P]\!]_E)}{k_2}\quad {\text {for }}k=\frac{k_1}{k_2}{\text { and }}k_1,k_2 \in {\mathbb {N}}\\&[\![(P_1)_D:(P_2)]\!]_E=([\![P_1]\!]_E)_{[\![D]\!]_E} : ([\![P_2]\!]_E) \\&[\![p]\!]_E=p \\&[\![p \cdot c_i + D]\!]_E=p \cdot E(c_i) + ([\![D]\!]_E) \end{aligned}$$

where

$$\begin{aligned}\pi _{one}(y)=\left\{ \begin{array}{l l} 1, &{} {\text {if }}\;y=1 \\ 0, &{}{\text {otherwise}}\\ \end{array} \right. , \pi _{zero}(y)=\left\{ \begin{array}{l l} 1, &{}{\text {if }}\;y=0\\ 0, &{}{\text {otherwise}}\\ \end{array} \right. . \end{aligned}$$

To illustrate the calculus, consider the Bernoulli distribution with parameter $p\in {\mathbb {Q}}_{[0,1]}$. We have $bern^p=(one)_p:zero$, where $[\![bern^p]\!]_{\emptyset }(y)=\{p {\text { if }} y=1;1-p {\text { if }}y=0;0 \,\,\, {\text {otherwise}}\}$.

The binomial distribution can be obtained as a sum of n independent Bernoulli distributions of the same parameter. Given a random variable with a binomial distribution with parameters (n, p), if n is sufficiently large and p sufficiently small then this approximates a Poisson distribution with parameter $n\cdot p$.

4.1 Operations on distributions

In this section, we define a set of operations on pmfs needed to define the semantics of the calculus. We conclude the section by showing that these operations are sufficient to represent pmfs with finite support in ${\mathbb {N}}$.

Definition 10

Let $\pi _1:{\mathbb {N}} \rightarrow [0,1]$, $\pi _2:{\mathbb {N}} \rightarrow [0,1]$ be two pmfs. Assume $p \in {\mathbb {Q}}_{[0,1]}$, $y \in {\mathbb {N}}$, $k_1 \in {\mathbb {N}}$ and $k_2 \in {\mathbb {N}}_{>0}$, then we define the following operations on pmfs:

• The sum or convolution of $\pi _1$ and $\pi _2$ is defined as

  $$\begin{aligned} (\pi _1 + \pi _2)(y)=\sum _{ (y_i,y_j) \in {{\mathbb {N}}\times {\mathbb {N}}} \, s.t. \, y_i+y_j=y} \pi _1(y_i)\pi _2(y_j). \end{aligned}$$

• The minimum of $\pi _1$ and $\pi _2$ is defined as

  $$\begin{aligned}&min(\pi _1,\pi _2)(y)\\&\quad =\sum _{(y_i,y_j) \in {\mathbb {N}}\times {\mathbb {N}}\,s.t.\,min(y_i,y_j)\,=y} \pi _{1}(y_i)\pi _{2}(y_j). \end{aligned}$$

• The multiplication of $\pi _1$ by the constant $k_1$ is defined as

  $$\begin{aligned} (k_1 \pi _1) (y)= \left\{ \begin{array}{l l} \pi _1\left( \frac{y}{k_1}\right) , &{}\quad {\text { if}} \,\frac{y}{k_1} \in {\mathbb {N}}\\ 0, &{} {\text {otherwise}}\\ \end{array} \right. \end{aligned}$$

• The division of $\pi _1$ by the constant $k_2$ is defined as

  $$\begin{aligned} \frac{\pi }{k_2}(y)=\sum _{y_i \in {\mathbb {N}}\,s.t.\, y=\lfloor y_i / k_2 \rfloor }\pi (y_i). \end{aligned}$$

• The convex combination of $\pi _1$ and $\pi _2$, for $y \in {\mathbb {N}}$, is defined as

  $$\begin{aligned} ( (\pi _1)_p : (\pi _2) )(y)= p\pi _1 (y)+ (1-p) \pi _2 (y). \end{aligned}$$

Example 4

Consider the following pmf $\pi _1: {\mathbb {N}}\rightarrow [0,1]$

$$\begin{aligned} \pi _{1}(y_1)= \left\{ \begin{array}{l l} \frac{1}{6}, &{}\quad {\text {if} }\;y_1=3\\ \frac{5}{6}, &{}\quad {\text {if }}\;y_1=0\\ 0, &{}\quad {\text {otherwise}}\\ \end{array} \right. \end{aligned}$$

and the following pmf $\pi _2: {\mathbb {N}}\rightarrow [0,1]$

$$\begin{aligned} \pi _{2}(y_2)= \left\{ \begin{array}{l l} \frac{1}{2}, &{}\quad {\text {if }}\;y_2=5\\ \frac{1}{2}, &{}\quad {\text {if }}\;y_2=1\\ 0, &{}\quad {\text {otherwise}}\\ \end{array} \right. \end{aligned}$$

Then the sum of $\pi _1$ and $\pi _2$ is:

$$\begin{aligned} (\pi _1 + \pi _2)(y)= \left\{ \begin{array}{l l} \frac{1}{12}, &{}\quad {\text {if }}\;y=8\\ \frac{5}{12}, &{}\quad {\text {if }}\;y=5\\ \frac{1}{12}, &{}\quad {\text {if }}\;y=4\\ \frac{5}{12}, &{}\quad {\text {if }}\;y=1\\ 0, &{}\quad {\text {otherwise}}\\ \end{array} \right. \end{aligned}$$

Example 5

Consider the pmfs $\pi _1$ and $\pi _2$ of Example 4 then

$$\begin{aligned} min(\pi _1,\pi _2)(y)= \left\{ \begin{array}{l l} \frac{1}{12}, &{}\quad {\text {if} }\;y=3\\ \frac{1}{12}, &{}\quad {\text {if }}\;y=1\\ \frac{5}{6}, &{}\quad {\text {if }}\;y=0\\ 0, &{}\quad {\text {otherwise}}\\ \end{array} \right. \end{aligned}$$

Example 6

Consider the pmf $\pi _2$ of Example 4, then

$$\begin{aligned} 2\pi _{2}(y)= \left\{ \begin{array}{l l} \frac{1}{2}, &{}\quad {\text {if }}\;y=10\\ \frac{1}{2}, &{}\quad {\text {if }}\;y=2\\ 0, &{}\quad {\text {otherwise}}\\ \end{array} \right. \end{aligned}$$

Example 7

Consider the following formula

$$\begin{aligned} P_1=(one)_{0.001\cdot c + 0.2}:(4\cdot one)+(2 \cdot one)_{0.4}:(3\cdot one), \end{aligned}$$

with set of environmental variables $V=\{c \}$ and an enviroment E such that $V(P_1)\subseteq dom(E)$. Then, according to Definition 10 we have that

$$\begin{aligned}{}[P_1]_E(y)=\left\{ \begin{array}{ll} (0.001\cdot [c]_E+0.2)\cdot 0.4, &{} {\text {if }}\;y=3\\ (0.001\cdot [c]_E+0.2)\cdot 0.6, &{} {\text {if }}\;y=4\\ (1-(0.001\cdot [c]_E+0.2))\cdot 0.4, &{} {\text {if }}\;y=6\\ (1-(0.001\cdot [c]_E+0.2))\cdot 0.6, &{} {\text {if }}\;y=7\\ 0, &{} {\text {otherwise}}\\ \end{array} \right. \end{aligned}$$

The convex combination operator is the only one that is not closed with respect to pmfs whose support is a single point. Lemma 1 shows the associativity of the convex distribution.

Lemma 1

Given probability mass functions $\pi _1$, $\pi _2 : {\mathbb {N}} \rightarrow [0,1]$, $p_1,p_2,p_3,p_4 \in [0,1]$ and $k \in {\mathbb {Q}}_{\ge 0}$, then the following equations hold:

• $k(({\pi _1})_p : \pi _2)=(k{\pi _1})_p : (k\pi _2)$

• $(({\pi _1})_{p_1} : \pi _2)_{p_2}: \pi _3 = ({\pi _1})_{p_3} :( ({\pi _2})_{p_4}:\pi _3)$ iff $p_3=p_1 p_2$ and $p_4=\frac{(1-p_1)p_2}{1-p_1 p_2}$

• $({\pi _1})_p : \pi _2=({\pi _2})_{1-p} : \pi _1$

• $({\pi _1})_p : \pi _1=\pi _1$.

Proof

We need to prove each statement.

Case $k(({\pi _1})_p : \pi _2)=(k{\pi _1})_p : (k\pi _2)$.

For $y \in {\mathbb {N}}$ we have that

$$\begin{aligned}&k({(\pi _1)}_p : \pi _2)(y)\\&\quad = \sum _{y_i \in {\mathbb {N}} s.t. \lfloor ky_i \rfloor =y} (p\pi _1 (y_i)+(1-p)(\pi _2(y_i)))\\&\quad =\sum _{y_i \in {\mathbb {N}}\, s.t.\, \lfloor ky_i \rfloor =y} (p\pi _1 (y_i))\\&\qquad +\sum _{y_i \in {\mathbb {N}}\, s.t. \, \lfloor ky_i \rfloor =y}((1-p)(\pi _2(y_i)))\\&\quad =p\cdot \sum _{y_i\in {\mathbb {N}}\, s.t.\, \lfloor ky_i \rfloor =y} (\pi _1 (y_i))+(1-p)\cdot \sum _{y_i \in {\mathbb {N}}\, s.t. \,\lfloor ky_i \rfloor =y}((\pi _2(y_i)))\\&\quad =(k{\pi _1})_p : (k\pi _2))(y) \end{aligned}$$

Case $(({\pi _1})_{p_1} : \pi _2)_{p_2}: \pi _3 = ({\pi _1})_{p_3} :( ({\pi _2})_{p_4}:\pi _3)$ iff $p_3=p_1 p_2$ and $p_4=\frac{(1-p_1)p_2}{1-p_1 p_2}$.

For $y \in {\mathbb {N}}$ we have that

$$\begin{aligned}&(({\pi _1}_{p_1} : \pi _2)_{p_2}:\pi _3)(y)\\&\quad =p_2(p_1 \pi _1(y)+(1-p_1)\pi _2(y))+(1-p_2)\pi _3(y)\\&({\pi _1}_{p_3} :( {\pi _2}_{p_4}:\pi _3))(y)\\&\quad =p_3\pi _1(y)+ (1-p_3)(p_4 \pi _2(y)+(1-p_4)\pi _3(y)) \end{aligned}$$

These are equal if

$$\begin{aligned}&p_1 p_2=p_3\\&p_4-p_3 p_4=p_2-p_1p_2 \\&1-p_2=(1-p_3)(1-p_4) \end{aligned}$$

and these conditions are satisfied if and only if $p_3=p_1 p_2$ and $p_4=\frac{(1-p_1)p_2}{1-p_1 p_2}$.

Case $({\pi _1})_p : \pi _2=({\pi _2})_{1-p} : \pi _1$.

For $y \in {\mathbb {N}}$ by Definition 10 it holds that

$$\begin{aligned}&(({\pi _1})_p : \pi _2)(y)=p \pi _1(y) + (1-p) \pi _2 (y)\\&\quad =(1-p)\pi _2 (y)+ p \pi _1 (y)= (({\pi _2})_{1-p} : \pi _1)(y) \end{aligned}$$

Case $({\pi _1})_p : \pi _1=\pi _1$.

For $y \in {\mathbb {N}}$ by Definition 10 it holds that

$$\begin{aligned}&(({\pi _1})_p : \pi _1)(y)=p \pi _1(y) + (1-p) \pi _1 (y)\\&\quad =(p+1-p)\pi _1 (y)= \pi _1(y) \end{aligned}$$

$\square$

Having formally defined all the operations on pmfs, we can finally state the following proposition guaranteeing that the semantics of any formula of the calculus is a pmf.

Proposition 1

Given P, a formula of the calculus defined in Definition  8, and an environment E such that $V(P)\subseteq dom(E)$, then $[\![P]\!]_E$ is a pmf.

Proof

The proof is by structural induction on the structure of P with basic cases $[\![one]\!]_E=\pi _{one}$ and $[\![zero]\!]_E=\pi _{zero}$, which are pmfs by definition for any E. $\square$

The following theorem shows that our calculus is complete with respect to finite support distributions.

Theorem 6

For any pmf $f:{\mathbb {N}}\rightarrow [0,1]$ with finite support there exists a formula P such that $[\![P]\!]_\emptyset =f$.

Proof

Given a pmf $f:{\mathbb {N}} \rightarrow [0,1]$ with finite support $J=(z_1,\ldots ,z_{|J|})$ we can define $P=(z_1\cdot one)_{f(z_1)}:((z_2\cdot one)_{\frac{f(z_2)}{1-f(z_1)}}:(\ldots :((z_i\cdot one)_{\frac{f(z_i)}{\prod _{j=1}^{i-1}(1-f(z_j))}}:\ldots :((z_n\cdot one )))))$. Then, $[\![P]\!]_\emptyset =f$. $\square$

Proof of Theorem 6 relies only on a subset of the operators, but the other operators are useful for composing previously defined pmfs.

5 CRN implementation

We show how the operators of the calculus can be realized by operators on CRSs. The resulting CRSs produce the required distributions at steady state, that is, in terms of the steady state distribution of the induced CTMC. Thus, we need to consider a restricted class of CRNs that always stabilize and that can be incrementally composed. The key idea is that each such CRN has output species that cannot act as a reactant in any reaction, and hence the counts of those species increase monotonically.¹ This implies that the optimized CRSs shown in Sect. 3.2 cannot be used compositionally.

5.1 Non-reacting output CRSs (NRO-CRSs)

Since in the calculus presented in Definition 8 we consider only finite support pmfs, in this section we are limited to finite state CTMCs. This is important because some results valid for finite state CTMCs are not valid in infinite state spaces. Moreover, any pmf with infinite support on natural numbers can always be approximated under the $L^1$ norm (see Corollary 1).

Given a CRS $C=(\varLambda ,R,x_0)$, we call the non-reacting species of C the subset of species $\varLambda _r \subseteq \varLambda$ such that given $\lambda _r \in \varLambda _r$ there does not exist $\tau \in R$ such that $r^{\lambda _r}_{\tau }>0$, where $r^{\lambda _r}_{\tau }$ is the component of the source complex of the reaction $\tau$ relative to $\lambda _r$, that is, $\lambda _r$ is not a reactant in any reaction. Given C we also define a subset of species, $\varLambda _o \subseteq \varLambda$, as the output species of C. Output species are those whose limit distribution is of interest. In general, they may or may not be non-reacting species; they depend on the observer and on what he/she is interested in observing.

Definition 11

A non-reacting output CRS (NRO-CRS) is a tuple $C=(\varLambda ,\varLambda _o,R,$ $x_0)$, where $\varLambda _o \subseteq \varLambda$ are the output species of C such that $\varLambda _o \subseteq \varLambda _r$, where $\varLambda _r$ are the non-reacting species of C.

NRO-CRNs are CRSs in which the output species are produced monotonically and cannot act as a reactant in any reaction. A consequence of Theorem 3 is the following lemma, which shows that this class of CRNs can approximate any pmf with support on natural numbers, up to an arbitrarily small error.

Lemma 2

For any probability mass function $f : {\mathbb {N}}^{m} \rightarrow [0,1]$ there exists a NRO-CRS such that the joint limit distribution of its output species approximates f with arbitrarily small error under the $L^1$ norm. The approximation is exact if the support of f is finite.

Proof

This lemma is a consequence of Theorems 3 and 4. In fact, by construction, all CRSs used in those theorems are non-reacting output. $\square$

5.1.1 NRO-CRS operators

A NRO-CRS operator is a NRO-CRS such that, given as input the output of certain NRO-CRSs, it produces as output a (set of) species that at steady state implement a given operation. We define the following NRO-CRS operators and show their correctness.

Definition 12

Let $C_1=(\varLambda _1,\varLambda _{o_1},R_1,x_{0_1})$ and $C_2=(\varLambda _2,\varLambda _{o_2},R_2,$ $x_{0_2})$ be NRO-CRSs such that $\varLambda _1 \cap \varLambda _2=\emptyset$. Then, for $\lambda _{o_1}\in \varLambda _{o_1},\lambda _{o_2}\in \varLambda _{o_2}$, $\{\lambda _{out},\lambda _z,\lambda _{r_1},$ $\lambda _{r_2}\}\cap (\varLambda _1\cup \varLambda _2)=\emptyset$, $k\in {\mathbf {N}},p\in [0,1]$, we define the following NRO-CRS operators:

$$\begin{aligned}&Sum(C_1,\lambda _{o_1},C_2,\lambda _{o_2},\lambda _{out})\\&\quad = (\varLambda _1\cup \varLambda _2\cup \{ \lambda _{out}\},\{\lambda _{out}\},R_1\cup R_2 \cup \\&\qquad \{ \lambda _{o_1} \rightarrow \lambda _{out}, \lambda _{o_2} \rightarrow \lambda _{out}\},x_0)\\&Min(C_1,\lambda _{o_1},C_2,\lambda _{o_2},\lambda _{out}) \\&\quad =(\varLambda _1\cup \varLambda _2\cup \{ \lambda _{out}\},\{\lambda _{out}\},R_1 \cup R_2 \cup \\&\qquad \{ \lambda _{o_1} + \lambda _{o_2}\rightarrow \lambda _{out}\},x_0) \\&\qquad Mul(C_1,\lambda _{o_1},k,\lambda _{out})\\&\quad = (\varLambda _1 \cup \{\lambda _{out}\},\{\lambda _{out}\},R_1\cup \\&\qquad \left\{ \lambda _{o1} \rightarrow \underbrace{\lambda _{out}+\cdots +\lambda _{out}}_{k\, times}\right\} ,x_{0}) \\&Div(C_1,\lambda _{o_1},k,\lambda _{out})\\&\quad = (\varLambda _1 \cup \{\lambda _{out}\},\{\lambda _{out}\},R_1\cup \\&\qquad \left\{ \underbrace{\lambda _{o_1}+\cdots +\lambda _{o_1}}_{k \, times}\rightarrow \lambda _{out}\right\} ,x_{0})\\&Con(C_1,\lambda _{o_1},C_2,\lambda _{o_2},p,\lambda _{out})\\&\quad = (\varLambda _1\cup \varLambda _2\cup \{\lambda _z,\lambda _{r_1},\lambda _{r_2}, \lambda _{out}\},\{\lambda _{out}\},R_1\cup R_2 \cup \\&\qquad \,\{\lambda _z \rightarrow ^{p} \lambda _{r_1},\lambda _{z} \rightarrow ^{1-p} \lambda _{r_2},\\&\qquad \lambda _{o_1} +\lambda _{r_1} \rightarrow \lambda _{r_1}+\lambda _{out}, \lambda _{o_1} +\lambda _{r_2} \rightarrow \lambda _{r_2}+\lambda _{out}\},x_0) \end{aligned}$$

where $x_{0}(\lambda )={\left\{ \begin{array}{ll} x_{0_1}(\lambda )&{}\quad {\text { if }}\lambda \in \varLambda _1\\ x_{0_2}(\lambda )&{}\quad {\text { if }}\lambda \in \varLambda _2\\ 1&{}\quad {\text { if }}\lambda = \lambda _z\\ 0&{}\quad {\text { otherwise}}\\ \end{array}\right. }$

Theorem 7

Let $C_1=(\varLambda _1,\varLambda _{o_1},R_1,x_{0_1})$ and $C_2=(\varLambda _2,$ $\varLambda _{o_2},R_2,$ $x_{0_2})$ be NRO-CRSs such that $\varLambda _1 \cap \varLambda _2=\emptyset$. Then, for $\lambda _{o_1}\in \varLambda _{o_1},\lambda _{o_2}\in \varLambda _{o_2},$ $\lambda _{out}\not \in \varLambda _1\cup \varLambda _2$, $k\in {\mathbf {N}},p\in [0,1]$ we have:

$$\begin{aligned}&\pi _{\lambda _{out}}^{Sum(C_1,\lambda _{o_1},C_2,\lambda _{o_2},\lambda _{out})}= \pi _{\lambda _{o_1}}^{C_1}+\,\pi _{\lambda _{o_2}}^{C_2}\\&\pi _{\lambda _{out}}^{Min(C_1,\lambda _{o_1},C_2,\lambda _{o_2},\lambda _{out})}=min\left( \pi _{\lambda _{o_1}}^{C_1},\pi _{\lambda _{o_2}}^{C_2}\right) \\&\pi _{\lambda _{out}}^{Mul(C_1,\lambda _{o_1},k,\lambda _{out})}= k\pi _{\lambda _{o_1}}^{C_1}\\&\pi _{\lambda _{out}}^{ Div(C_1,\lambda _{o_1},k,\lambda _{out})}= \frac{\pi _{\lambda _{o_1}}^{C_1}}{k}\\&\pi _{\lambda _{out}}^{Con(C_1,\lambda _{o_1},C_2,\lambda _{o_2},p,\lambda _{out})}= \left( \pi _{\lambda _{o_1}}^{C_1}\right) _{p}:\pi _{\lambda _{o_2}}^{C_2} \end{aligned}$$

The proof of Theorem 7 is not trivial, and is given in the next subsection. The key difficulties lie in the fact that we need to compose stochastic processes and show that the resulting process has the required properties.

Example 8

We consider the pmfs $\pi _1$ and $\pi _2$ of Example 4. Using the results of Theorem 3 we build the CRSs $C_1$ and $C_2$ such that $\lambda _{out_1}$ and $\lambda _{out_2}$, unique output species of $C_1$ and $C_2$ respectively, admit as steady state distribution exactly $\pi _1$ and $\pi _2$. $C_1=(\{\lambda _z,\lambda _1,\lambda _{1,1},\lambda _{o_1}\},$ $\{\lambda _{o_1}\},R,x_0)$ has the following reactions

$$\begin{aligned} \lambda _{z} \rightarrow ^{\frac{1}{6}} \lambda _{1,1}; \quad \lambda _{z} \rightarrow ^{\frac{5}{6}} \emptyset ; \quad \lambda _{1}+\lambda _{1,1} \rightarrow ^{1} \lambda _{1,1} + \lambda _{o_1}; \end{aligned}$$

where $\emptyset$ is the empty set and $x_0$ is such that: $x_0(\lambda _1)=3,\, x_0(\lambda _z)=1, \, x_0(\lambda _{1,1})=0,\, x_0(\lambda _{o_1})=0$.

The CRS $C_2$ has the following reactions

$$\begin{aligned}&\lambda _{z'} \rightarrow ^{\frac{1}{2}} \lambda _{1',1'};\quad \lambda _{z'} \rightarrow ^{\frac{1}{2}} \lambda _{2',2'};\\&\lambda _{1'}+\lambda _{1',1'} \rightarrow ^{1} \lambda _{1',1'} + \lambda _{o_2};\\&\lambda _{2'}+\lambda _{2',2'} \rightarrow ^{1} \lambda _{2',2'} + \lambda _{o_2}; \end{aligned}$$

with initial condition $x_0$ such that: $x_0(\lambda _{1'})=5,\, x_0(\lambda _z')=1, \, x_0(\lambda _{1',1'})=0, \,x_0(\lambda _{2',2'})=1,\,x_0(\lambda _{2'})=5,\,x_0(\lambda _{o_2})=0$. Then, applying the Sum operator circuit, we add the following reactions

$$\begin{aligned} \lambda _{o_1} \rightarrow ^{1} \lambda _{out}; \quad \lambda _{o_2} \rightarrow ^{1} \lambda _{out}; \end{aligned}$$

$Sum(C_1,\lambda _{o_1},C_2,\lambda _{o_2},\lambda _{out})$ has unique output species $\lambda _{out}$, whose limit distribution, $\pi _{\lambda _{out}}$, is equal to $\pi _{1}+\,\pi _{2}$ described in Example 4.

In what follows, we present in extended form the operator for convex combination, and introduce a new operator, which implements the convex distribution with external inputs ($ConE(\cdot )$).

Considering $C_1$ and $C_2$, as previously, then we need to derive a CRS operator $Con(C_1,\lambda _{o_1},C_2,\lambda _{o_2},p,\lambda _{out})$ such that $\pi _{\lambda _{out}}=(\pi ^{C_1}_{\lambda _{o_1}})_p : (\pi ^{C_2}_{\lambda _{o_2}})$. That is, at steady stade, $\lambda _{out}$ equals $\pi ^{C_1}_{\lambda _{o_1}}$ with probability p and $\pi ^{C_2}_{\lambda _{o_2}}$ with probability $1-p$. This can be done by using Theorem 4 to generate a bi-dimensional synthetic coin with output species $\lambda _{r_1},\lambda _{r_2}$ such that their joint limit distribution is

$$\begin{aligned} \pi _{\lambda _{r_1},\lambda _{r_2}}(y_1,y_2)= {\left\{ \begin{array}{ll} p &{} \quad {\text {if }}\;y_1=1\;{\text { and }}\;y_2=0\\ 1-p &{} \quad {\text {if }}\;y_1=0\;{\text { and }}\;y_2=1\\ 0 &{} \quad {\text {otherwise}}\\ \end{array}\right. }. \end{aligned}$$

That is, $\lambda _{r_1}$ and $\lambda _{r_2}$ are mutually exclusive at steady state. Using these species as catalysts in $\tau _3:\lambda _{o_1}+\lambda _{r_1}\rightarrow \lambda _{r_1}+ \lambda _{out}$ and $\tau _4:\lambda _{o_2}+\lambda _{r_2}\rightarrow \lambda _{r_2}+ \lambda _{out}$ we have exactly the desired result at steady state.

Example 9

Consider the following NRO-CRSs $C_1=(\{\lambda _{o_1}\},\{\lambda _{o_1}\},\{\},x_{0_1} )$ and $C_2=(\{\lambda _{o_2}\},\{\lambda _{o_2}\},\{\},x_{0_2} )$, with initial condition $x_{0_1}(\lambda _{o_1})=10$ and $x_{0_2}(\lambda _{o_2})=20$. Then, the operator $Con(C_1,\lambda _{o_1},C_2,\lambda _{o_2},0.3,\lambda _{out})$ implements the operation $\pi _{\lambda _{out}}=(\pi ^{C_1}_{\lambda _{o_1}})_{0.3}( \pi ^{C_2}_{\lambda _{o_2}})$ and it is given by the following reactions:

$$\begin{aligned}&\lambda _z \rightarrow ^{0.3} \lambda _{r_1};\quad \lambda _z \rightarrow ^{0.7} \lambda _{r_2};\\&\lambda _{r_1} + \lambda _{o_1} \rightarrow \lambda _{r_1}+\lambda _{out};\quad \lambda _{r_2} + \lambda _{o_2} \rightarrow \lambda _{r_2}+\lambda _{out} \end{aligned}$$

with initial condition $x_0$ such that $x_0(\lambda _z)=1$, $x_0(\lambda _{r_1})=x_0(\lambda _{r_2})=x_0(\lambda _{out})=0.$

Let $C_1,C_2$ be as above and $f=p_0+p_1\cdot c_1+\cdots +p_n\cdot c_n$ with $p_1,\ldots ,p_n \in {\mathbb {Q}}_{[0,1]}$, $V=\{c_1,\ldots ,c_n\}$ a set of environmental variables, and E, an environment such that $V\subseteq dom(E)$. Then, computing a CRS operator $ConE(C_1,\lambda _{o_1},C_2,\lambda _{o_2},$ $f(E(V)),\lambda _{out})$ such that $\pi _{\lambda _{out}}=(\pi ^{C_1}_{\lambda _{o_1}})_{f(E(V))} : (\pi ^{C_2}_{\lambda _{o_2}})$ is a matter of extending the previous circuit. First of all, we can derive the CRS to compute f(E(V)) and $1-f(E(V))$ and memorize them in some species. This can be done as f(E(V)) is semi-linear (Chen et al. 2014). Then, as $f(E(V))\le 1$ by assumption, we can use these species as catalysts to determine the output value of $\lambda _{out}$, as in the previous case. As shown in Sect. 5.2, this circuit, in the case of external inputs, introduces an arbitrarily small, but non-zero, error, due to the fact that there is no way to know when the computation of f(E(V)) terminates.

Example 10

Consider the following NRO-CRSs $C_1=(\{\lambda _{o_1}\},\{\lambda _{o_1} \},\{\},x_{0_1} )$ and $C_2=(\{\lambda _{o_2}\},\{\lambda _{o_2}\},\{\},x_{0_2} )$, with initial condition $x_{0_1}(\lambda _{o_1})=10$ and $x_{0_2}(\lambda _{o_2})=20$. Then, consider the following functions $f(E(c))=E(c)$, where E is a partial function assigning values to c, and it is assumed $0.001\le E(c)\le 1$ and that $E(c)\cdot 1000 \in {\mathbb {N}}$. Then, the operator $ConE(C_1,\lambda _{o_1},C_2,\lambda _{o_2},f,\lambda _{out})$, implements the operation $\pi _{\lambda _{out}}=(\pi ^{C_1}_{\lambda _{o_1}})_{E(c)}(\pi ^{C_2}_{\lambda _{o_2}})$ and it is given by the following reactions:

$$\begin{aligned}&\tau _{1}:\lambda _{c} \rightarrow ^{k_1} \lambda _{\mathrm {Cat}_1}+\lambda _{\mathrm {Cat}_2};\, \tau _{2}:\lambda _{Tot}+\lambda _{\mathrm {Cat}_2}\rightarrow ^{k_1} \emptyset \\&\tau _{3}:\lambda _{z}+\lambda _{\mathrm {Cat}_1} \rightarrow ^{k_2} \lambda _{1};\, \tau _{4}:\lambda _{z}+\lambda _{Tot}\rightarrow ^{k_2} \lambda _2\\&\tau _{5}:\lambda _{o_1}+\lambda _{1} \rightarrow ^{k_2} \lambda _{1}+\lambda _{out};\, \tau _{6}:\lambda _{o_2}+\lambda _{2} \rightarrow ^{k_2} \lambda _{2}+\lambda _{out} \end{aligned}$$

where $\lambda _{c},\lambda _{\mathrm {Cat}_1},\lambda _{\mathrm {Cat}_2},\lambda _{z},\lambda _{1}$ and $\lambda _{2}$ are auxiliary species with initial condition $x_0$ such that $x_0(\lambda _{\mathrm {Cat}_1})=x_0(\lambda _{\mathrm {Cat}_2})$ $=x_0(\lambda _{1})=x_0(\lambda _{2})=0,$ $x_0(\lambda _{Tot})=1000, x_0(\lambda _{z})=1$, $x_0(\lambda _{c})=E(c)\cdot 1000$ and $k_1 \gg k_2$. Reactions $\tau _1,\tau _2$ implement f(E(c)) and $1-f(E(c))$ and store these values in $\lambda _{Cat_1}$ and $\lambda _{Tot}$. These are used in reactions $\tau _3$ and $\tau _4$ to determine the probability that the steady state value of $\lambda _{out}$ is going to be determined by reaction $\tau _5$ or $\tau _6$.

5.2 Correctness of the CRS-operators

We prove the correctness of Theorem 7. For the sake of simplicity, we consider only the Sum operator, as other operators have similar proofs. The key idea of the proof is to make use of Eq. (1) to show that the resulting CRS implements the desired operation at steady state.

Proposition 2

Let $C_1=(\varLambda _1,\varLambda _{o_1},R_1,x_{0_1}), \, C_2=(\varLambda _2,$ $\varLambda _{o_2},R_2,x_{0_2})$ be NRO-CRSs such that $\varLambda _1 \cap \varLambda _2=\emptyset$ and $\{\lambda _{out}\}\cap (\varLambda _1 \cup \varLambda _2)=\emptyset$. Then for $\lambda _{o_1} \in \varLambda _{o_1}$ and $\lambda _{o_2} \in \varLambda _{o_2}$ the CRS $Sum(C_1,\lambda _{o_1},C_2,\lambda _{o_2},\lambda _{out})=C_c$ is such that $\pi ^{C_c}_{\lambda _{out}}= \pi _{\lambda _{o_1}}^{C_1}+\,\pi _{\lambda _{o_2}}^{C_2}$.

Proof

Consider the counting processes $J_{\lambda _{o_1}}^{C_c}$ and $J_{\lambda _{o_2}}^{C_c}$, acording to the stocahstic model introduced in (1), which give the number of molecules of $\lambda _{o_1}$ and $\lambda _{o_2}$ produced until time t in $C_c$. Using Eq. (1) we have

$$\begin{aligned}&J_{\lambda _{o_1}}^{C_c}(t)=\sum _{\tau \in R_1\cup R_2 \cup \{\tau _{s_1}, \tau _{s_2}\}}p^{\lambda _{o_1}}_{\tau }Y_{\tau }\left( \int _0^t \! \alpha _{\tau }(X^{C_c}(s)) \, {\mathrm {d}}s\right) \\&J_{\lambda _{o_2}}^{C_c}(t)=\sum _{\tau \in R_1\cup R_2 \cup \{\tau _{s_1}, \tau _{s_2}\}}p^{\lambda _{o_2}}_{\tau }Y_{\tau }\left( \int _0^t \! \alpha _{\tau }(X^{C_c}(s)) \, {\mathrm {d}}s\right) \end{aligned}$$

where $p^{\lambda _{o_1}}_{\tau }$ and $p^{\lambda _{o_2}}_{\tau }$ represent the number of molecules of $\lambda _{o_1}$ and $\lambda _{o_2}$ produced by the occurrence of reaction $\tau$. Recall that $\tau _{s_1}$ and $\tau _{s_2}$ are such that $\tau _{s_1}: \lambda _{o_1}\rightarrow \lambda _{out}$ and $\tau _{s_2}: \lambda _{s_2}\rightarrow \lambda _{out}$ and $\varLambda _1 \cap \varLambda _2=\emptyset$. As a consequence, $p^{\lambda _{o_1}}_{\tau _{s_1}}=p^{\lambda _{o_1}}_{\tau _{s_2}}=p^{\lambda _{o_2}}_{\tau _{s_1}}=p^{\lambda _{o_2}}_{\tau _{s_2}}=0$ and we can write

$$\begin{aligned} J_{\lambda _{o_1}}^{C_c}(t)=&\sum _{\tau \in R_1\cup R_2 \cup \{\tau _{s_1}, \tau _{s_2}\}}p^{\lambda _{o_1}}_{\tau }Y_{\tau }\left( \int _0^t \! \alpha _{\tau }(X^{C_c}(s)) \, {\mathrm {d}}s\right) \\ =&\sum _{\tau \in R_1}p^{\lambda _{o_1}}_{\tau }Y_{\tau }\left( \int _0^t \! \alpha _{\tau }(X^{C_c}(s)) \, {\mathrm {d}}s\right) \end{aligned}$$

and

$$\begin{aligned}&J_{\lambda _{o_2}}^{C_c}(t) \\&\quad =\sum _{\tau \in R_1\cup R_2 \cup \{\tau _{s_1}, \tau _{s_2}\}}p^{\lambda _{o_2}}_{\tau }Y_{\tau }\left( \int _0^t \! \alpha _{\tau }(X^{C_c}(s)) \, {\mathrm {d}}s\right) \\&\quad =\sum _{\tau \in R_2 }p^{\lambda _{o_2}}_{\tau }Y_{\tau }\left( \int _0^t \! \alpha _{\tau }(X^{C_c}(s)) \, {\mathrm {d}}s\right) \end{aligned}$$

Moreover, $r_{\tau _{s_1}}^{\lambda }=p_{\tau _{s_1}}^{\lambda }=r_{\tau _{s_2}}^{\lambda }=p_{\tau _{s_2}}^{\lambda }=0$ for any $\lambda \in \varLambda -\{\lambda _{out},\lambda _{o_1},\lambda _{o_2}\}$, that is, $\tau _{s_1}$ and $\tau _{s_2}$ do not produce or consume any species in $\varLambda -\{\lambda _{out},\lambda _{o_1},\lambda _{o_2}\}$. As a consequence, because $x_0(\lambda )=x_{0_1}(\lambda )$ for all $\lambda \in \varLambda _1-\{\lambda _{o_1}\}$, we have $\int _0^t \! \alpha _{\tau }(X^{C_c}(s)) \, {\mathrm {d}}s=\int _0^t \! \alpha _{\tau }(X^{C_1}(s)) \, {\mathrm {d}}s$ for all $\tau \in R_1$ . In exactly the same way, it is possible to show that the same relation holds for $\lambda _{o_2}$ with respect to $X^{C_2}$, and as a consequence it is also true that $\int _0^t \! \alpha _{\tau }(X^{C_c}(s)) \, {\mathrm {d}}s=\int _0^t \! \alpha _{\tau }(X^{C_2}(s)) \, {\mathrm {d}}s$ for all $\tau \in R_2$. As a result:

$$\begin{aligned}&\sum _{\tau \in R_1}p^{\lambda _{o_1}}_{\tau }Y_{\tau }\left( \int _0^t \! \alpha _{\tau }(X^{C_c}(s)) \, {\mathrm {d}}s\right) =\sum _{\tau \in R_1}p^{\lambda _{o_1}}_{\tau }Y_{\tau }\left( \int _0^t \! \alpha _{\tau }(X^{C_1}(s)) \, {\mathrm {d}}s\right) \\&\sum _{\tau \in R_2 }p^{\lambda _{o_2}}_{\tau }Y_{\tau }\left( \int _0^t \! \alpha _{\tau }(X^{C_c}(s)) \, {\mathrm {d}}s\right) =\sum _{\tau \in R_2 }p^{\lambda _{o_2}}_{\tau }Y_{\tau }\left( \int _0^t \! \alpha _{\tau }(X^{C_2}(s)) \, {\mathrm {d}}s\right) \end{aligned}$$

Considering that $\lambda _{o_1}$ is an output species in $C_1$ and $\lambda _{o_2}$ is an output species in $C_2$, that is, NRO-CRSs, then for any $\tau \in R_1$ we have that $\upsilon _{\tau }^{\lambda _{o_1}}=p_{\tau }^{\lambda _{o_1}}$ and for any $\tau \in R_2$ $\upsilon _{\tau }^{\lambda _{o_2}}=p_{\tau }^{\lambda _{o_2}}$. As a consequence:

$$\begin{aligned} X^{C_1}_{\lambda _{o_1}}(t)=&X^{C_s}_{\lambda _{o_1}}(0)+\sum _{\tau \in R_1}{p}^{\lambda _{o_1}}_{\tau } Y_{\tau }\left( \int _0^t \! \alpha _{\tau }(X^{C_1}(s)) \, {\mathrm {d}}s \right) \\ =&X^{C_1}_{\lambda _{o_1}}(0)+J^{C_c}_{\lambda _{o_1}}(t)\\ X^{C_2}_{\lambda _{o_2}}(t)=&X^{C_2}_{\lambda _{o_2}}(0)+\sum _{\tau \in R_2}{p}^{\lambda _{o_2}}_{\tau } Y_{\tau }\left( \int _0^t \! \alpha _{\tau }(X^{C_2}(s)) \, {\mathrm {d}}s \right) \\ =&X^{C_2}_{\lambda _{o_2}}(0)+J^{C_c}_{\lambda _{o_2}}(t) \end{aligned}$$

According to the fact that in the composed NRO-CRS $\lambda _{out}$ is produced only by $\tau _{s_1}$ and $\tau _{s_2}$ such that $p_{\tau _{s_1}}^{\lambda _{out}}=p_{\tau _{s_2}}^{\lambda _{out}}=1$, and that $\lambda _{out}$ is not consumed in any reaction, and its initial molecular count is 0. Then, it is possible to write:

$$\begin{aligned} X^{C_c}_{\lambda _{out}}(t)=0 + Y_{\tau _{s_1}}\left( \int _{0}^t \! \alpha _{\tau _c}(X^{C_c}(s))\! {\mathrm {d}}s\right) +Y_{\tau _{s_2}}\left( \int _{0}^t \! \alpha _{\tau _c}(X^{C_c}(s)\right) \! {\mathrm {d}}s) \end{aligned}$$

In the same way we can define the stochastic model for the number of molecules of $\lambda _{o_1}$ or $\lambda _{o_2}$ present in $C_c$ at a given time, as given by the number of molecules produced minus the number of molecules consumed. As $\lambda _{o_1}$ and $\lambda _{o_2}$ are consumed only by $\tau _{s_1}$ and $\tau _{s_2}$, and they are not reactant in any other reaction, we have:

$$\begin{aligned} X^{C_c}_{\lambda _{o_1}+\lambda _{o_2}}(t)=&X^{C_c}_{\lambda _{o_1}}(0) + X^{C_c}_{\lambda _{o_2}}(0) \\&+\sum _{\tau \in R_1}{p}^{\lambda _{o_1}}_{\tau } Y_{\tau }\left( \int _0^t \! \alpha _{\tau }(X^{C_1}(s)) \, {\mathrm {d}}s \right) \\&+\sum _{\tau \in R_2}{p}^{\lambda _{o_2}}_{\tau } Y_{\tau }\left( \int _0^t \! \alpha _{\tau }(X^{C_2}(s)) \, {\mathrm {d}}s \right) \\&- Y_{\tau _{s_1}}\left( \int _{0}^t \! \alpha _{\tau _c}(X^{C_c}(s))\! {\mathrm {d}}s\right) \\&-Y_{\tau _{s_2}}\left( \int _{0}^t \! \alpha _{\tau _c}(X^{C_c}(s))\! {\mathrm {d}}s\right) \\ =&X^{C_1}_{\lambda _{o_1}}(t)+ X^{C_2}_{\lambda _{o_2}}(t)-X^{C_c}_{\lambda _{out}}(t) \end{aligned}$$

because $X^{C_c}_{\lambda _{o_1}}(0)=X^{C_1}_{\lambda _{o_1}}(0)$ and $X^{C_c}_{\lambda _{o_2}}(0)=X^{C_1}_{\lambda _{o_2}}(0)$ by assumption.

The set of reachable states from $x_0$ in $X^{C_c}$ is finite because the set of reachable states from $x_{0_1}$ in $X^{C_1}$ and from $x_{0_2}$ in $X^{C_2}$ are finite by assumption and $\tau _c$, in a finite time, can fire only a finite number of times. This implies that $X^{C_c}(t)$ for $t \rightarrow \infty$ will reach a bottom strongly connected component (BSCC) of the underlying graph of the state space, with probability 1 in finite time, because of a well known result of CTMC theory (Kwiatkowska et al. 2007). In a BSCC, any pair of configurations $x_1$ and $x_2$ are such that $x_1 \rightarrow ^* x_2$ and $x_2 \rightarrow ^* x_1$. Therefore, any configuration x in any BSCC reachable by $X^{C_c}$ from $x_0$ is such that $x(\lambda _{o_1})=x_0(\lambda _{o_2})=0$, because in a configuration $x_i$ where $x_i(\lambda _{o_1})>0$ or $x_i(\lambda _{o_2})>0$ it is always possible to reach a configuration $x_j$ where $x_j(\lambda _{o_1})=x_i(\lambda _{o_1})-1$ or $x_j(\lambda _{o_2})=x_i(\lambda _{o_2})-1$ and $x_j(\lambda _{out})=x_i(\lambda _{out})+1$, but then there is no way to reach $x_i$ from $x_j$ because $\lambda _{out}$ is not reactant in any reaction in $R_1\cup R_2 \cup \{\tau _{s_1},\tau _{s_2}\}$. Therefore

$$\begin{aligned}&\lim _{t \rightarrow \infty }Prob(X^{C_c}_{\lambda _{o_1}+\lambda _{o_2}}(t)=\,0|X^{C_c}(0)=x_0)=1 \implies \\&\lim _{t \rightarrow \infty } Prob( X^{C_1}_{\lambda _{o_1}}(t)+X^{C_2}_{\lambda _{o_2}}(t)-X^{C_c}_{\lambda _{out}}(t) = 0|\\&X^{C_c}(0)=x_0,X^{C_1}(0)=x_{0_1},\\&X^{C_2}(0)=x_{0_2})=1 \implies \\&\lim _{t \rightarrow \infty }Prob(X^{C_c}_{\lambda _{out}}(t)=X^{C_1}_{\lambda _{o_1}}(t)+X^{C_2}_{\lambda _{o_2}}(t)|\\&X^{C_c}(0)=x_0,X^{C_1}(0)=x_{0_1},\\&X^{C_2}(0)=x_{0_2})=1 \end{aligned}$$

This concludes the proof. $\square$

5.3 Compiling into the class of NRO-CRSs

Given a formula P as defined in Definition 8, then $[\![P]\!]_E$ associates to P and an environment E a pmf. We now define a translation of P, T(P), into the class of NRO-CRSs that guarantees that the unique output species of T(P), at steady state, approximates $[\![P]\!]_E$ with arbitrarily small error for any environment E such that $V(P)\subseteq dom(E)$. In order to define such a translation we need the following renaming operator.

Definition 13

Given a CRS $C=(\varLambda ,R,x_0)$, for $\lambda _t \in \varLambda$ and $\lambda _1 \not \in \varLambda$ we define the renaming operator $C\{\lambda _1 \leftarrow \lambda _t \}=C_c$ such that $C_c=((\varLambda -\{\lambda _t\})\cup \{ \lambda _1\}, R\{\lambda _1 \leftarrow \lambda _t \},x_0')$, where $R\{\lambda _1 \leftarrow \lambda _t \}$ substitutes any occurrence of $\lambda _t$ with an occurrence of $\lambda _1$ for any $\tau \in R$ and $x_0'(\lambda )=\{x_0(\lambda )$ if $\lambda \ne \lambda _t; \, x_0(\lambda _t)$ if $\lambda =\lambda _1 \}$.

This operator produces a new CRS where any occurrence of a species is substituted with an occurrence of another species previously not present.

Definition 14

(Translation into NRO-CRSs) Define the mapping T by induction on syntax of formulae P:

$$\begin{aligned} T (one)=&(\{\lambda _{out}\},\{\lambda _{out}\},\emptyset ,x_0)\quad {\text {with }}x_0(\lambda _{out})=1;\\ \quad T (zero)=&(\{\lambda _{out}\},\{\lambda _{out}\},\emptyset ,x_0)\,\,\, {\text {with} } x_0(\lambda _{out})=\,0;\\ T (P_1+P_2)=&Sum(T(P_1)\{\lambda _{o_1} \leftarrow \lambda _{out}\},\\&\lambda _{o_1},T(P_2)\{\lambda _{o_2} \leftarrow \lambda _{out}\},\lambda _{o_2},\lambda _{out});\\ T (k\cdot P)=&Div(Mul(T(P)\{\lambda _{o} \leftarrow \lambda _{out}\},\\&\quad \lambda _{o},k_1,\lambda _{out})\{\lambda _{o'} \leftarrow \lambda _{out}\}),\lambda _{o'},k_2,\lambda _{out});&\nonumber \\ T (min(P_1,P_2)=&Min(T(P_1)\{\lambda _{o_1} \leftarrow \lambda _{out}\},\\&\lambda _{o_1},T(P_2)\{\lambda _{o_2} \leftarrow \lambda _{out}\},\lambda _{o_2},\lambda _{out});&\nonumber \\ T ((P_1)_D:P_2)=&\left\{ \begin{array}{ll} Con(T(P_1)\{\lambda _{o_1} \leftarrow \lambda _{out}\},\lambda _{o_1},T(P_2)\{\lambda _{o_2} \leftarrow \lambda _{out}\},\\ \quad \quad \quad \lambda _{o_2},D,\lambda _{out}), \quad \quad {\text {if }}D=p\\ ConE(T(P_1)\{\lambda _{o_1} \leftarrow \lambda _{out}\},\lambda _{o_1},T(P_2)\{\lambda _{o_2} \leftarrow \lambda _{out}\},\\ \quad \quad \quad \lambda _{o_2},D,\lambda _{out}), \quad \quad {\text { if }}D=p+\sum _{i=1}^{m}p_i\cdot c_i\\ \end{array} \right. \end{aligned}$$

for $m>1$, $k \in {\mathbb {Q}}_{>0}$, $k_1,k_2 \in {\mathbb {N}}$ such that $k=\frac{k_1}{k_2}$ and formulae $P_1,P_2$, which are assumed to not contain species $\lambda _{o_1},\lambda _{o_2}$.

Example 11

Consider the formula $P_1=(one)_{0.001\cdot c + 0.2}$ $(4\cdot one)+(2 \cdot one)_{0.4}(3\cdot one)$ of Example 7, and an environment E such that $0.000125 \le E(c) \le 1$ and suppose $E(c) \cdot 800 \in {\mathbb {N}}$. We show how the translation defined in Definition 14 produces a NRO-CRS C with output species $\lambda _{out}$ such that $\pi _{\lambda _{out}}=[\![P_1]\!]_E$. Consider the following NRO-CRSs $C_1,C_2,C_3,C_4$ defined as $C_1=(\{\lambda _{c_1}\},\{\lambda _{c_1}\},\{\},x'_0)$ with $x_0(\lambda _{c_1})=1$, $C_2=(\{\lambda _{c_2}\},$ $\{\lambda _{c_2}\},\{\},$ $x_0)$ with $x_0(\lambda _{c_2})=1$, $C_3=(\{\lambda _{c_3}\},\{\lambda _{c_3}\},\{\},$ $x_0)$ with $x_0(\lambda _{c_3})=1$, and $C_4=(\{\lambda _{c_4}\},\{\lambda _{c_4}\},\{\},$ $x_0)$ with $x_0(\lambda _{c_2})=1$. Then, we have that :

$$\begin{aligned} C^c_1=&ConE(C_1,\lambda _{c_1},Mul(C_2,\lambda _{c_2},4,\lambda _{out})\{\lambda _{o_2}\leftarrow \lambda _{out}\},\\&\lambda _{o_2},0.001\cdot c+0.2,\lambda _{out_1})\\ C^c_2=&Con(Mul(C_3,\lambda _{c_3},2,\lambda _{out})\{\lambda _{o_3}\leftarrow \lambda _{out}\},\lambda _{o_3},\\ {}&Mul(C_4,\lambda _{c_4},3,\lambda _{out})\{\lambda _{o_4}\leftarrow \lambda _{out}\},\lambda _{o_4},0.4,\lambda _{out_2}) \end{aligned}$$

are such that $\pi _{\lambda _{out_1}}=\left\{ \begin{array}{ll} (0.001\cdot [\![c]\!]_E+0.2), &{} {\text {if }}\;y=1\\ 1-(0.001\cdot [\![c]\!]_E+0.2),&{}{\text {if }}\;y=4 \\ 0,&{}{\text {otherwise}}\\ \end{array} \right.$, and   $\pi _{\lambda _{out_2}}=\left\{ \begin{array}{ll} 0.4, &{}{\text {if }}\;y=2\\ 0.6, &{}{\text {if }}\;y=3\\ 0, &{}{\text {otherwise}}\\ \end{array} \right.$. Then, consider the CRS $C=Sum(C^c_1\{\lambda _{t_1\leftarrow \lambda _{out_1}}\},\lambda _{t_1},$ $C^c_2\{\lambda _{t_2\leftarrow \lambda _{out_2}}\},$ $\lambda _{t_2},\lambda _{out})$ and we have $\pi _{\lambda _{out}}=[\![P_1]\!]_E$ with arbitrarily small error. The reactions of C are shown below

$$\begin{aligned}&Mul{\text { on inputs}} \{\tau _1: \lambda _{C_2}\rightarrow 4\lambda _{o_1};\quad \tau _2: \lambda _{C_3}\rightarrow 2\lambda _{o_2};\\&\tau _3: \lambda _{C_4}\rightarrow 3\lambda _{o_3}.\\&C^c_1 \left\{ \begin{array}{l} \tau _4: \lambda _{env}\rightarrow ^{k} \lambda _{cat_1}+\lambda _{cat_2};\\ \tau _5: \lambda _{cat_1} + \lambda _{z} \rightarrow \lambda _{1} \\ \tau _6: \lambda _{cat_2} + \lambda _{tot} \rightarrow ^{k} \emptyset ;\\ \tau _7: \lambda _{tot} + \lambda _{z} \rightarrow \lambda _{2}\\ \tau _8: \lambda _{1} + \lambda _{o_1} \rightarrow \lambda _{o_1} + \lambda _{out_1};\\ \tau _9: \lambda _{2} + \lambda _{o_2} \rightarrow \lambda _{o_2}+ \lambda _{out_1} \end{array} \right. \\&C^c_2 \left\{ \begin{array}{l} \tau _{10}: \lambda _{z_1}\rightarrow ^{0.6} \lambda _{r_1};\\ \tau _{11}: \lambda _{z_1} \rightarrow ^{0.4} \lambda _{r_2}\\ \tau _{12}: \lambda _{r_1} + \lambda _{o_3}; \rightarrow \lambda _{r_1}+\lambda _{out_2};\\ \tau _{13}:\lambda _{r_2} + \lambda _{o_4} \rightarrow \lambda _{r_2}+\lambda _{out_2} \end{array} \right. \\&Sum \left\{ \tau _{14}:\lambda _{out_1}\rightarrow \lambda _{out};\quad \tau _{15}:\lambda _{out_2}\rightarrow \lambda _{out} \right. \end{aligned}$$

for $k\gg 1$ and initial condition such that $x_0(\lambda _{env})=E(c)\cdot 800$, $x_0(\lambda _{tot})=800$, $x_0(\lambda _{z})=x_0(\lambda _{z_1})=x_0(\lambda _{z_2})=1=x_0(\lambda _{c_1})=x_0(\lambda _{c_2})=x_0(\lambda _{c_3})=x_0(\lambda _{c_4})=1$, and all other species initialized with 0 molecules.

Proposition 3

For any formula P we have that T(P) is a NRO-CRS.

Proof

The proof is by structural induction. The base cases are T(zero) and T(one), which are NRO-CRSs by definition. Assuming $T(P_1)$ and $T(P_2)$ are NRO-CRNs then application of operators of sum, Mul, Div, Min, Con and ConE on these CRSs produces a NRO-CRNs by definition of the operators. $\square$

Given a formula P and an environment E such that $V(P)\subseteq dom(E)$, the following theorem guarantees the soundness of T(P) with respect to $[\![P]\!]_E$. In order to prove the soundness of our translation we consider the measure of the multiplicative error between two pmfs $f_1$ and $f_2$ with values in ${\mathbb {N}}^m$, $m>0$ as $e_m(f_1,f_2)=\max _{n \in {\mathbb {N}}^m}\min (\frac{f_1(n)}{f_2(n)},\frac{f_2(n)}{f_1(n)})$.

Theorem 8

(Soundness) Given a formula P and $\lambda _{out}$, unique output species of T(P), then, for an environment E such that $V(P)\subseteq dom(E)$, it holds that $\pi ^{T(P)}_{\lambda _{out}} = [\![P]\!]_E$ with arbitrarily small error under multiplicative error measure.

The proof follows by structural induction.

Remark 3

A formula P is finite by definition, so Theorem 8 is valid because the only production rule which can introduce an error is $(P_1)_D:(P_2)$ in the case $D\ne p_0$, and we can always find reaction rates to make the total probability of error arbitrarily small. Note that, by using the results of Soloveichik et al. (2008), it would also be possible to show that the total error can be kept arbitrarily small, even if a formula is composed from an unbounded number of production rules. This requires small modifications to the ConE operator following ideas in Soloveichik et al. (2008).

Observe that compositional translation, as defined in Definition 14, generally produces more compact CRNs with respect to the direct translation in Theorem 3, and in both cases the output is non-reacting, so the resulting CRN can be used for composition. For a distribution with support J direct translation yields a CRN with 2|J| reactions, whereas, for instance, the support of the sum pmf has the cardinality of the Cartesian product of the supports of the input pmfs.

6 Discussion

Our goal was to explore the capacity of CRNs to compute with distributions. This is an important goal because, when molecular interactions are in low number, as is common in various experimental scenarios (Qian and Winfree 2014), deterministic methods are not accurate, and stochasticity is essential for cellular circuits. Moreover, there is a large body of literature in biology where stochasticity has been shown to be essential and not only a nuisance (Eldar and Elowitz 2010). Our work is a step forward towards better understanding of molecular computation. In this paper we focused on error-free computation for distributions. It would be interesting to understand and characterize what would happen when relaxing this constraint. That is, if we admit a probabilistically (arbitrarily) small error, does the ability of CRNs to compute on distributions increase? Another interesting topic to investigate is whether we can relax the constraint that the output species are produced monotonically. In fact, this is a constraint that is generally not present in natural systems where species undergo production and degradation reactions. More specifically, we require that a CRN will reach a state where no reactions can happen. In terms of sampling from the distribution, this would require sampling an ensemble of cells since sampling a single cell would yield a single state. Also, we would like to address the problem if it is possible to implement distributions in CRNs without leaders (species being present with initial number of molecules equal to 1) and without knowing the precise initial number of molecules for each species. Our constructions, except for the uniform distribution, crucially rely on these assumptions, though may be challenging to obtain in technologies such as DNA strand displacement (Soloveichik et al. 2010). As a consequence, DNA implementation would become easier if these constraints can be removed. However, it is worth noting that, in a practical scenario, leaders can be thought of as single genes or localized structures (Qian and Winfree 2014), and there exist CRN techniques to produce given concentrations independently of initial conditions (Shinar and Feinberg 2010).