Computing Real Numbers with Large-Population Protocols Having a Continuum of Equilibria

The computation of real numbers provides a theoretical benchmark for an abstract model’s computational power. In 1936, Turing [32] considered computable numbers while researching a discrete model now known as the Turing machine. Then, Shannon [30] introduced the general-purpose analog computer (GPAC) to the world. The GPAC is a mathematical abstraction of Bush’s differential analyzer. Bush [8] developed this machine to obtain numerical solutions of ordinary differential equations (ODEs). After much revision [29, 27, 18, 17], we can now characterize GPACs by polynomial initial value problems (PIVPs):

where $\mathbf{p}$ is a vector of multivariate polynomials, $\mathbf{x}$ is a vector of variables, and $\mathbf{x}_{0}$ is its initial value. Bournez et al. [7] have shown that GPACs are provably equivalent to Turing machines in terms of computability and complexity. Chemical reaction networks (CRNs), as a restricted form of GPACs, are widely used in molecular programming applications. The deterministic semantic of CRNs usually assumes well-mixed solutions and mass-action kinetics. CRN dynamics can be modeled by polynomial systems with extra constraints [21] (see the preliminary section for more details). Although the constraints could seemingly weaken the computational power of CRNs, they are able to simulate GPACs. Hence, the CRN model is also Turing complete [13]. The key idea is to encode every variable $x$ (possibly having negative values) in a GPAC by the difference between a pair of positive variables $x^{+}$ and $x^{-}$ in the copycat CRN. That is, $x(t)=x^{+}(t)-x^{-}(t)$ for all $t\geq 0$.

A straightforward idea to “compute” a number $\alpha$ by a GPAC/CRN is to designate a variable $x$ such that $x(t)\to\alpha$ as $t\to\infty$. Further down this line, Huang et al. defined the class of real-time computable real numbers by chemical reaction networks [24], wherein the real-time computability notion essentially requires $x(t)$ to converge exponentially fast to $\alpha$. Huang et al. then demonstrate that GPACs and CRNs compute the same class of real numbers in real-time [23]. Notably, famous transcendental numbers like $e$ and $\pi$ are among this class. Later, Huang [22] added Euler’s $\gamma$ and Dottie number (the unique real root of the equation $\cos(x)=x$) to the list. Note that in order to get a meaningful measure of time, in Huang et al.’s work, all variables in ODEs associated with GPACs or CRNs are bounded. Hence, the notion of time aligns with the time parameter $t$ of the ODEs since no more than a linear speedup is allowed. This notion prevents the so-called Zeno phenomena [19] caused by inordinate speedup. Another, perhaps equivalent, metric is the arc length of $\mathbf{x}(t)$. A complexity theory [7] of analog computation can then be built upon this metric.

Also, note that the notion Huang et al. used has a descriptive nature: the initial values are all zero, and rate constants are integrals (rationals). Hence, one can only encode a finite amount of information in the initial value and rate constant. In this sense, a system that computes a number $\alpha$ also gives a finite description of $\alpha$. This idea will prove crucial for later discussions.

Population protocols (PPs) [2] can be treated as special CRNs with two inputs, two outputs, and a unit rate constant per reaction. Observe that population protocols conserve the population since interactions do not change agent counts. This property imposes a vital constraint on the population: the sum of all agents remains constant. After normalization, we can say that this sum is one without loss of generality.

Classical PPs are concerned with computing predicates over their input configuration in a discrete, stochastic setting. Therefore, computing real numbers does not make sense under this setting. In a series of papers [3, 5, 6], Bournez et al. developed a setting they termed a Large-Population Protocol (LPP). LPPs are a new analytical setting rather than a new computational model. They are still population protocols, but their behaviors are analyzed as the population grows to infinity. Bournez et al. aimed for asymptotic results independent of the initial configuration. Our approach departs from theirs in that dependence on initial values is considered. Note that we should approach with caution when defining what is meant by “initial configuration” since LPPs need to talk about different population sizes along the way.

The LPP setting provides a playground for computing real numbers by population protocols: a number $\nu\in[0,1]$ is said to be computable by an LPP if the proportion of agents in a set of marked states converges to $\nu$ over time as the population grows to infinity. An additional requirement is that the ODE induced by the balance equation of the LPP must have finitely many equilibria. This condition enforces (exponential) stability and detaches the dependence on initial values.

Bournez et al. specifically asked[6] “… is it possible to compute solutions of trigonometric equations? E.g., can we design a protocol insuring that, asymptotically, a ratio $\frac{\pi}{4}$ of the nodes are in some prescribed state?” Their answer is negative. They proved that, under their definition, LPPs compute exactly the algebraic numbers in [0,1]. Closely related is a result by Fletcher et al. [15], where a class of numbers, called Lyapunov numbers, are found to be the set of algebraic numbers. It is not a coincidence: both works involve ODEs with either finitely many or isolated equilibria. Hence, Tarski’s quantifier elimination over real closed fields can be performed. Alternatively, one can find more modern elimination algorithms involving Gröbner bases in standard algebraic geometry texts [10, 31].

Why is there a gap between LPPs (cannot compute transcendentals) and GPAC/CRN (can compute $e$ and $\pi$)? The following example helps explore the question.

It is evident that with initial values $F(0)=\frac{1}{2}$, $E(0)=\frac{1}{2}$, and $G(0)=0$, the solution of $F(t)\to\frac{1}{2}e^{-1}$ as $t\to\infty$. Correspondingly, in the LPP setting, if we let $F(0)=\lfloor\frac{N}{2}\rfloor$, $E(0)=\lfloor\frac{N}{2}\rfloor$, and $G(0)=0$ for a population of size $N$, we can observe from the experimental results in Figure 1 that the proportion of $F$ gets very close to $\frac{1}{2}e^{-1}$ when $N$ becomes larger. In fact, the so-called Kurtz’s Theorem [26], which says ODE solution $F(t)$ equals the proportion random variable $F^{(N)}(t)/N$ almost surely when $N$ goes to infinity, guarantees the limit is $\frac{1}{2}e^{-1}$. Although the CRN in the above example is not technically a PP, we can translate it into a probabilistic PP (Construction 1, see also [11], Theorem 1). We maintain the protocol in the current form for simplicity.

So the above LPP can “compute” $\frac{1}{2}e^{-1}$, a transcendental number! Is it a counter-example to Bournez et al.’s algebraic result? Not really. Since the system of ODEs in the example has the whole $F$-$G$ plane (with $E=0$) as equilibria, the system has a continuum of equilibria. Note that the key reason GPACs/CRNs in [23] can compute transcendentals is that they are systems with a continuum of equilibria.

The Motivating Example invites us to extend the notion of LPP-computable numbers: We can drop the finitary requirement on equilibria and therefore include systems with a continuum of equilibria. The seminal work [4] by Sanjay and Bernstein provides an in-depth study of this type of system.

The move has several consequences. An immediate one is that computation now depends on initial values: since the equilibria “stick together”, we need the initial values to determine which equilibrium the system converges to. Another one is stability: we lose asymptotic stability. As pointed out in [4], “… Since every neighborhood of a non-isolated equilibrium contains another equilibrium, a non-isolated equilibrium cannot be asymptotically stable. Thus asymptotic stability is not the appropriate notion of stability for systems having a continuum of equilibria.” However, it is inaccurate to immediately conclude that such systems must be fragile and sensitive to perturbations. Indeed, Sanjay and Bernstein investigated semistability, a more appropriate notion of stability for these systems.

We now state our new finding under the extended notion of LPP-computable numbers.

Our proof is constructive: We give an explicit algorithm to translate a bounded GPAC or CRN into an LPP. A critical step is the construction of a cubic form (homogeneously degree 3 polynomial) system that conserves the population. Then, we transform the system from the cubic form into the quadratic form. Finally, the system can be rewritten as a probabilistic population protocol.

The key feature of our construction is that we always guarantee the system is either CRN-implementable or PP-implementable. Bournez et al. constructed a system in [6] that claimed to compute any algebraic number. While their system is a GPAC, it is not a PP in general. As a by-product, our construction fixes this flaw.

We have discussed that GPACs can be characterized by polynomial initial value problems (PIVPs) in the previous section. To align with the ODE descriptions of GPACs, we will use the ODE characterization of CRNs and PPs throughout this paper. This choice clarifies the discussion of our construction and proof in the coming sections.

We start with a few notations and conventions in this paper.

A typical definition of a CRN or PP specifies a pair $(S,R)$, with $S$ being the set of species or states and $R$ being the set of reactions or interactions between states. In general, a reaction for an abstract CRN looks like $X+Y\xrightarrow{\mathmakebox{k}}X+W+Z$, where $k$ is the reaction constant. There is no restriction on the number of reactants (input) or products (output). Reactions in population protocols, however, are strictly limited to those with two inputs and two outputs: $X+Y\to W+Z.$

The so-called deterministic mass-action semantic of a CRN defines a polynomial differential system, which governs the CRN’s dynamics. One can find descriptions of this semantic in many works of literature (e.g., [23]). We will use these polynomial differential systems directly throughout the paper.

Now, we introduce notation and terminology for polynomials.

Homogeneous polynomials play a fundamental role in our construction. We call them forms for short.

We say a function $f(t)$ is implementable (generable) by a GPAC (CRN, LPP) if there exist a GPAC (CRN, LPP) and a variable (species, state) $x$ in it, such that $f(t)=x(t)$ for all $t$. We choose the word implementable over generable (as in [18]) to reflect the idea that we program a GPAC (CRN, LPP) to implement a function.

Note that if a variable $x$ has initial value $a$, we can introduce a new variable $y=x-a$ to make $y(0)=0$. Therefore, without loss of generality, we often assume $\mathbf{x}(0)=\mathbf{0}$ in this paper.

The restriction is rooted in reaction modeling: a negative term in $x^{\prime}$ means to destroy some molecular $x$ at some rate. However, a reaction cannot destroy a non-reactant, so $x$ must appear as a reactant in the reaction. Such restriction extends to Population Protocols together with some extra constraints.

In this paper, we will also call a CRN containing only unimolecular reactions of the form $X\to Y$ a unimolecular population protocol (UPP). Similarly, a CRN containing only termolecular reactions of the form $X+Y+Z\to U+V+W$ is called a termolecular population protocol (TPP). In general, a $k$-PP is a CRN that contains only reactions of the form $X_{1}+\cdots+X_{k}\to Y_{1}+\cdots+Y_{k}$. Generally, we have the following characterization.

From above, we see that turning a GPAC into a PP means having to dance with more and more chains. We must carefully appease all restrictions, which makes the translation process difficult.

We first revisit some useful concepts regarding PPs.

The balance equation is a description of the system’s dynamics. Intuitively, it says whenever $q_{1}$ and $q_{2}$ bump into each other, if there is a reaction $(q_{1},q_{2})\to(q_{3},q_{4})$ in the PP, then the pair $(q_{1},q_{2})$ turns into $(q_{3},q_{4})$; otherwise, we interpret $(q_{1},q_{2})$ as a null reaction and the agents remain unchanged. We adopt an asymmetric interpretation of reactions in PPs throughout the paper. That is, the ordered pairs $(q_{i},q_{j})$ and $(q_{j},q_{i})$, when used as reactants, do not necessarily result in the same products. The assumption is not essential: the choice is to align with the existing literature (e.g., [6]) for convenience.

The above balance equation is for deterministic PPs, where $\delta_{q_{1},q_{2},q_{3},q_{4}}$ is either 0 or 1. It is more convenient to use probabilistic transition rules in many scenarios. We introduce probabilistic LPPs (PLPPs).

The balance equation remains mostly the same as Equation (4). To simplify the notations, we often take $Q=[n]=\{1,2,\cdots,n\}$ and will still use $q_{i}$ or use terms like “state $i$” when referring to a state.

The ODE associated with a PLPP can be written as

where $\mathbf{x}=(x_{1}(t),\cdots,x_{n}(t))\in\mathbb{R}^{n}$ is a variable vector and its $i$-th component, $x_{i}$, tracks the proportion of the total population in state $i$ over time $t$. We will often call $\mathbf{x}=(x_{1},x_{2},\cdots,x_{n})$ or $\mathbf{x}$ an LPP (PLPP).

A simple observation results directly from the definition of LPP or PLPP.

Although obvious, the above is the single most important observation in this paper. Many constructions and proofs rely on it. This observation is frequently used to re-write the constant term in an ODE. We will rewrite a constant $c$ in either of the following ways, depending on our needs:

• •

  $c=c\cdot 1=c(x_{1}+\cdots+x_{n})$, or

• •

  $c=c\cdot 1\cdot 1=c(x_{1}+\cdots+x_{n})(x_{1}+\cdots+x_{n})$.

An intuitive way to think about equation 7:

• •

  The term $-2x_{r}$ represents all reactions that consume $x_{r}$. A complete set of rules for a PP must include all combinations of pairs $(x_{r},x_{k})$, or $(x_{k},x_{r})$, where $k$ iterates through all $[n]$. These combinations sum to $-2x_{r}$.

• •

  The $f(x)$ represents all reactions that produce $x_{r}$. This observation is critical in the constructive proof of our main theorem.

We now give a formal definition of LPP-computable numbers.

We extend the definition from [6]. Note that the definition also applies to large-population unimolecular protocols (LPUPs).

A major change is our removal of the finitary requirement on equilibria. A note originally in [6] discusses the rationale for the requirement of finitely many equilibria. This assumption is needed to “avoid pathological cases, in particular the case of idle systems $q\ q^{\prime}\rightarrow q\ q^{\prime}$ for all $q$ and $q^{\prime}$. Indeed, in idle systems, all initial states are equilibria, and such a system would compute any real of [0,1], depending on the initial configuration.”

This pathological behavior could occur, if one is allowed to initialize the system with any real numbers. However, we can prevent this case by restricting the initial configuration to rational numbers. Also note that by forcing the initial counts of a species $X$ to be a fraction of the total population $N$, we exclude those systems (e.g.,[12, 1]) that crucially depend on small counts in a very large population. These systems are not modeled correctly by ODEs. For instance, one can not initialize $X$ to have 100 molecules. Any state initialized to a constant amount will be a fraction that approaches zero, as $N$ grows to infinity. In our definition, such $X$ will be initialized to have zero molecules for a fixed $N$.

For the rest of the paper, whenever we mention LPP computable number, we mean the extended notion defined above, unless specified otherwise.

We will prove our main theorem in this section. It suffices to show that LPPs can simulate GPACs. Since one can transform a bounded GPAC into a bounded CRN with zeroed initial values [23], it is sufficient to consider only this type of CRN. We develop a four-stage algorithm that translates CRNs into LPPs.

• •

  Input: A CRN that computes a target number $\alpha$. We now enter a cascade of transformations:

• •

  Stage 1: into a CRN-implementable quadratic system,

• •

  Stage 2: into a TPP-implementable cubic form system,

• •

  Stage 3: into a PP-implementable quadratic form system, and

• •

  Stage 4: into a PLPP.

First, we discuss some basic operations as building blocks for our constructive proof. The notations used are as described in the preliminary section. The naming of operations “$\varepsilon$-trick” and “$\lambda$-trick” (the use of the Greek letters $\varepsilon$ and $\lambda$) originates in the proof of Lemma 5 in [6].

The operation is essentially a change of variable $t\to\varepsilon t$ and a direct result of the chain rule. The behavior of the ODE system will not change; it only alters the flow of time. Indeed, $t$ in the new system corresponds to $\varepsilon t$ of the original system. In our construction, we will use this operation to shrink the coefficient of an ODE system. The operation can be seen in [11, 6].

We can do a more general change of variable (composition.)

A time dilating operation, as above, controls the flow of time with the function $g(t)$. More accurately, $t$ in the new system corresponds to $G(t)=\int_{0}^{t}g(s)ds$ of the original system. If we want the new system to mimic the old system’s limiting behavior (i.e., to achieve the value $\lim_{t\to\infty}x_{1}(t)$), we must guarantee $G(\infty)=\infty$. Thus, time in the new system is nonterminating and can go to infinity in the eyes of the old system. Otherwise, the new system would terminate somewhere in the middle, relative to the old system. To obtain nonterminating time, it is enough to consider functions that converge slowly, e.g., $g(t)=\frac{1}{1+t}$. Alternatively, one could ensure $g$ is bounded from below, i.e., $g(t)>c>0$, for some constant $c$.

Time dilation for CRNs can be found in [25] (page 15, Theorem 3.7).

Next, we introduce an operation that shrinks variables to [0,1].

The biggest challenge in translating a CRN into an LPP is: the CRN does not necessarily preserve population, while the LPP requires it. In fact, our input CRN initializes at zero. However, $x_{1}$ needs to compute some number $\alpha$ as all other components remain positive throughout time; evidently, not preserving the population. A straightforward idea to “balance” the system, is to introduce a variable $x_{0}$, such that its derivative cancels all changes to the population size by other variables. That is, $x^{\prime}_{0}=-\sum_{i=1}^{n}x_{i}^{\prime}$.

However, the resulting $x_{0}$ is not CRN-implementable in general. As a result, it is not transformable to an LPP. The rationale is clear in the derivative of $x_{i}$. In general, $x_{i}^{\prime}$ has the form

A typical $p_{i}$ would contain at least one monomial $m$. Evidently, $x_{0}\not\in m$ since the variable $x_{0}$ is newly introduced. Then $m$ will contribute to a negative term in $x_{0}^{\prime}$ and $x_{0}\not\in m$, making $x_{0}$ not CRN-implementable and, in turn, not PP-implementable. That is why Bournez et al.’s construction in [6] (Lemma 2) fails to produce a population protocol in general. To address the this issue, we purposefully multiply $x_{0}$ with every $x_{i}^{\prime}$ and introduce the following operation.

What makes the operation different from a normal time dilation operation is, the function $x_{0}$ is not known in advance: It needs to adjust/customize according to the system it is dilating. In our construction, we also need to combine the operation with the $\lambda$-trick to achieve nonterminating time. This operation also turns a non-conservative system into a conservative CRN-implmentable one.

We will present our construction in the following subsections.

The following is a CRN version of [9] (Theorem 1).

With this result, we can assume a CRN’s associated polynomial system is quadratic with out loss of generality.