Computing and Bounding Equilibrium Concentrations in Athermic Chemical Systems

In general, chemical equilibria of complex chemical systems are not analytically solvable and numerical tools are required for analysis. Such tools include NUPACK [17] for thermodynamic analysis of nucleic-acid systems, as well as more abstract platforms that support domain-level abstraction and free energy specification [12] including via rule-based modeling [11], and software for computing steady-state concentrations of chemical reaction networks [7, 10]. However, engineers often want a deeper understanding of the equilibrium than what analytically opaque numerical calculations can provide. Moreover, we often seek to understand infinite classes of designs such as logic circuits constructed from gate modules or parameterized constructions. For example, there is a growing body of work on leak in DNA-based systems, exhibiting a family of schemes parameterized by a “redundancy parameter” meant to decrease leak arbitrarily at the cost of additional system components [13, 5, 15, 1, 16, 14]. It is proven that producing off-target species necessarily decreases the overall number of separate complexes – the change controlled by the redundancy parameter – and thus incurs an entropic penalty. However, due to system complexity, the relationship between this rigorously proven thermodynamic unfavorability and the actual concentrations of off-target species often remains implicit.

While properties of equilibria of abstract coupled chemical reactions have been extensively explored in the chemical reaction network theory literature (e.g., [6], Chapter 14, for detailed-balance equilibria), and full explicit-parameter schemes (e.g., [9]) can characterize the entire space of equilibria, these analytical approaches have severe limitations. For large systems the resulting explicit formulas are unwieldy and offer no effective guidance on how to choose the parameters so that off-target concentrations remain below a desired bound.

In this paper we are interested in the following problem: Our chemical species are complexes (termed polymers) made up of indivisible units called monomers. We assume there are finitely many possible polymers, and among these we are given a set of on-target species that we want to have desired equilibrium concentrations, and all other off-target species to have some sufficiently low equilibrium concentrations. Our task is to determine a consistent (detailed-balance) equilibrium, and therefore the concentrations of the monomers that would lead to this equilibrium. Note that our interest in setting equilibrium concentrations rather than solving for them (based on initial concentrations or total monomer concentrations, for example) is misaligned with most computational approaches. For instance, it can be seen as the reverse problem to NUPACK 3’s concentrations tool which takes the concentrations of the monomers (strands) and returns the corresponding equilibrium concentrations of the polymers (complexes).

While having the flexibility to set monomer concentrations may appear to simplify the problem in the same way that finding some detailed-balance equilibrium is easier than computing the one consistent with input monomer concentrations, the complexity comes from simultaneously ensuring that off-target polymers remain in low concentration. Consider the natural approach of taking logarithms of all concentrations, thereby converting the detailed-balance equations (balancing each reaction) into a linear system amenable to standard linear solvers. When we fix the concentrations of on-target polymers, the remaining (off-target) concentrations are typically under-determined. The linear system then describes an unbounded affine subspace, and there is no obvious way to extract upper bounds on off-target species’ concentrations without leaving the linear framework. Our solution is an iterative algorithm that assigns concentrations to off-target polymers in decreasing order of concentration, ensuring that each off-target species remains below a desired threshold concentration when possible. Importantly, terminating the algorithm at any iteration still provides valid upper bounds for all remaining off-target polymers.

In the most general formulation of the monomer-polymer equilibrium concentration problem, the polymer free energies can be assigned arbitrarily incorporating binding strength, geometric constraints, etc. In this work, we focus on the simpler athermic case rather than tackling the problem in its full generality. Our allowed polymers are such that all possible reactions between them are enthalpy-neutral. This model is consistent with systems of strong fully-complementary DNA domains in which domain-level bonds can only switch binding partners but not de-hybridize. Thermodynamic Binding Network (TBN) [5, 2] saturated configurations capture this condition, but our setting is more general without a built-in notion of domains (binding sites).

The main result of this paper is Algorithm 1 and Theorem 5.4 showing how starting with desired concentrations of on-target polymers (already in detailed balance), we can set the concentrations of off-target polymers to satisfy detailed balance and thus thermodynamic equilibrium. In Section 4 we explain the apparent difficulties in balancing reactions which our approach needed to overcome.

If, rather than computing exact concentrations of off-target polymers, it is sufficient to bound them, then we refer the reader to Section 6. In Section 7 we apply our framework to the analysis of systems in the TBN model specifically, connecting the combinatorial notions of stability and entropy loss in the TBN model to equilibrium concentrations. In Section 8, we show applications of our method to the analysis of a simple TBN AND gate, as well as a parameterized family of signal propagation systems (translator cascades) from prior work. For the translator cascade, we argue that tuning concentrations of on-target polymers according to our framework is essential for leak to decrease exponentially with the redundancy parameter. We conclude with a discussion of future work (Section 9), including a formulation of new combinatorial conditions in the TBN model to make our framework easily applicable.

Let $\mathbb{N}$ denote the set of nonnegative integers. Given a finite set $𝒜$, we define ${\mathbb{N}}^{𝒜}$ as the set of functions $f:𝒜\to \mathbb{N}$.

A multiset $ℳ$ over the finite set $𝒜$ is described by its counting function $f_{ℳ}\in {\mathbb{N}}^{𝒜}$, where for each element $a\in 𝒜$, the value $f_{ℳ}(a)$ indicates how many times $a$ appears in $ℳ$. We often write $ℳ\in {\mathbb{N}}^{𝒜}$ to mean that $ℳ$ is a multiset over $𝒜$, and we denote the count of $a\in 𝒜$ in $ℳ$ by $ℳ[𝒶]$. The notation $a\in ℳ$ means that $ℳ[𝒶]\ge \mathcal{1}$. The cardinality of a multiset $ℳ$, denoted $|ℳ|$, is the total number of elements in the multiset $|ℳ|={\sum }_{𝒶\in 𝒜}ℳ[𝒶]$. For two multisets $ℳ$ and ${ℳ}^{\prime }$ over $𝒜$, we define their union $ℳ+{ℳ}^{\prime }$ as the multiset whose counting function is the pointwise sum $(ℳ+{ℳ}^{\prime })(𝒶)=ℳ[𝒶]+{ℳ}^{\prime }[𝒶]\text{ for all }𝒶\in 𝒜$. For example, let $ℳ=\{𝒶,𝒶,𝒷,𝒸\}$. Then, $ℳ[𝒶]=\mathcal{2}$, $ℳ[𝒷]=\mathcal{1}$, $ℳ[𝒸]=\mathcal{1}$, and $ℳ[𝒹]=\mathcal{0}$ for all $d\notin \{a,b,c\}$. Also, the cardinality of $ℳ$ is $|ℳ|=\mathcal{2}+\mathcal{1}+\mathcal{1}=\mathcal{4}$. Finally, note that $ℳ$ could also be written as the union $\{a,b\}+\{a,c\}$, for example.

The linear combination of multisets with nonnegative integers is defined analogously: For multisets ${ℳ}_{\mathcal{1}},\mathrm{\dots },{ℳ}_{𝓃}$ and $a_1,\mathrm{\dots },a_n\in \mathbb{N}$, ${\sum }_{i=1}^n{a}_i\cdot {ℳ}_{𝒾}$ corresponds to the counting function ${\sum }_{i=1}^n{a}_i\cdot {ℳ}_{𝒾}[𝒶]$. Let ${ℳ}_{\mathcal{1}},{ℳ}_{\mathcal{2}}\in {\mathbb{N}}^{𝒮}$ be two multisets over the same set $𝒮$. The difference ${ℳ}_{\mathcal{1}}-{ℳ}_{\mathcal{2}}$ is defined as the multiset $ℳ\in {\mathbb{N}}^{𝒮}$ such that for every $a\in 𝒮$, $({ℳ}_{\mathcal{1}}-{ℳ}_{\mathcal{2}})[𝒶]={ℳ}_{\mathcal{1}}[𝒶]-{ℳ}_{\mathcal{2}}[𝒶]$ provided that ${ℳ}_{\mathcal{1}}[𝒶]\ge {ℳ}_{\mathcal{2}}[𝒶]$ for all $a\in 𝒮$.

We also define the intersection of a multiset with a set. Given a multiset $ℳ$ over $𝒜$ and a subset $𝒮\subset 𝒜$, the intersection $ℳ\cap 𝒮$ is a multiset over $𝒜$ defined by $(ℳ\cap 𝒮)[𝒶]=ℳ[𝒶]$ if $a\in 𝒮$, otherwise $(ℳ\cap 𝒮)[𝒶]=\mathcal{0}$. For instance, $\{a,a,b,c\}\cap \{a,c\}=\{a,a,c\}$.

The main object of this paper is an abstract model of monomers and polymers, motivated by systems in DNA nanotechnology. This model captures how simple indivisible molecules (monomers) combine to form complexes (polymers) under specific physical and chemical constraints.

Let ${𝐱}^{\mathrm{𝟎}}\in {(0,1)}^{{𝚿}^{\mathrm{𝟎}}}$ represent the vector of concentrations for all monomers, and let $𝐱\in {(0,1)}^{𝚿}$ represent the vector of concentrations for all polymers (also called configuration). The relationship between monomer and polymer concentrations is governed by mass conservation. Specifically, we require

where $𝐀\in {\mathbb{N}}^{|{𝚿}^{\mathrm{𝟎}}|\times |𝚿|}$ is a matrix such that each entry $A_{ij}$ specifies the number of monomers of type $i$ in polymer $j$.

For example, the polymer $P=\{m_1,m_1,m_2,m_3\}$ contains two copies of $m_1$, one of $m_2$, and one of $m_3$. Note that we will be interested in cases where the set of polymers of interest $𝚿$ is a finite (proper) subset of all possible polymers over ${𝚿}^{\mathrm{𝟎}}$.

In DNA nanotechnology the monomers are typically DNA strands with different sequences. Polymers are analogous to a multistranded DNA structure composed of multiple DNA strands. We use the term polymer rather than “complex” in order to better emphasize their composition from monomers and to be consistent with the TBN literature.

To model the equilibrium behavior of such systems, we use the free energy formulation in the notation of Dirks et al. [4], where the equilibrium concentrations are obtained by minimizing the following free energy function (corresponding to the pseudo-Helmholtz free energy used throughout chemical reaction network theory literature [8, 6]):

where $x_P$ denotes the concentration of polymer $P$, and ${\mathrm{\Omega }}_P$ is its partition function corresponding to the exponential of the polymer’s negative free energy. The minimization is subject to the mass conservation constraint given in Equation 1.

In this work, we focus on athermic systems where all interactions are equally favored enthalpically. Thus we assume that ${\sum }_{P\in 𝚿}{x}_P\cdot \log {\mathrm{\Omega }}_P$ is constant for all configurations $𝐱$ satisfying mass conservation (Equation 1), yielding a simplified cost function that is entirely entropic:

This function serves as the objective of our optimization problem. Its minimizer under the constraint of Equation 1 corresponds to the equilibrium concentration of polymers.

Understanding the equilibrium requires us to formalize how polymers can transform into one another. We do so by defining reconfigurations and the reactions they induce.

We occasionally use Greek alphabets to denote the reactions. Note that while $\cong$ is a binary relation, a reaction $\alpha :{ℳ}_{\mathcal{1}}\to {ℳ}_{\mathcal{2}}$ is an ordered pair representing a directed transformation between the multisets. Our arguments will refer to reactants (left-hand side) and products (right-hand side) of particular reactions.¹¹1All reactions are of course reversible, but we treat each direction as a separate reaction.

An important property of the minimum of Equations 2 and 3 is that it satisfies detailed balance over reactions [8, 6]. For us (Equation 3), for any single reaction $\alpha :{ℳ}_{\mathcal{1}}\to {ℳ}_{\mathcal{2}}$, the equilibrium concentrations satisfy:

We say that a reaction $\alpha$ is balanced when this equality holds. This notion of balance will be central to our characterization of equilibrium concentrations formed by on-target polymers, which we explore next.

The fact that balance leads to the minimum of the free energy is well-established, and we present its sketch in Appendix B for the completeness of the paper.

We now define the set $𝒮$ of on-target polymers – intuitively, these are the high-concentration polymers whose equilibrium concentration we set as input to our algorithm.

While presented in terms of restricting $𝒮$, another (perhaps more apropos) interpretation of the first condition is that we are restricting $𝚿$ to be only those polymers that can be obtained via a reconfiguration of polymers over $𝒮$. The interpretation of the second condition is that every reaction among polymers in $𝒮$ is in detailed balance. More precisely:

If $𝒮=𝚿$ then we are done, being guaranteed that all reactions are in detailed balance. However, the case of interest is where we do not know the concentration exponents of all polymers – making it the goal of this paper to compute them or bound them to ensure that the equilibrium concentration of off-target polymers is small.

We now define canonical reactions, which are the reactions with reactants from $𝒮$. These reactions can generate polymers outside of $𝒮$ from polymers in $𝒮$. Each canonical reaction has key quantities we call imbalance and novelty which will play a key role in our algorithm.

Given the definition of multiset subtraction, the novelty $l(\alpha )$ is always non-negative.

Intuitively, the imbalance of the reaction represents how far from detailed balance it is if we consider only the polymers whose concentrations are known. The novelty is the number of off-target (new) polymers produced by the same reaction.

Intuitively, a large $k(\alpha )$ means a larger bias of the reaction toward its reactants, i.e., elements of $𝒮$, prior to assigning concentration exponents to off-target polymers. This is desired since it implies less pressure to make off-target polymers and gives us more room to maneuver in assigning concentrations to them.

Intuitively, a larger novelty $l(\alpha )$ means that the canonical reaction is more entropically favored to its products. Note that the term ${\sum }_{P\in 𝚿}{x}_P\log x_P$ in Equation 3 is minimized when the concentration is more “spread out” over the different species (like Shannon entropy), and thus there is an entropic force to generate new polymers. A large $l(\alpha )$ is thus undesired. As justified in the subsequent results, the ratio $k(\alpha )/l(\alpha )$ captures the effective tradeoff between imbalance and novelty.

Note that summing canonical reactions gives another canonical reaction. While on the one hand this leads to infinitely many canonical reactions, on the other hand, this additive property allows us to analyze the set of all canonical reactions using a Hilbert basis of finite size. We explain how to employ the Hilbert basis in our algorithm to avoid the infinity of canonical reactions in Appendix A.

The notion of stability of $𝒮$, defined below, captures the idea that on-target polymers are in higher concentration compared to off-target polymers. Recall that concentration exponents $\mu$ are $\le 1$ for polymers in $𝒮$ (Definition 3.1). As shown later by Theorem 5.4, the following definition ensures that all concentration exponents of polymers outside of $𝒮$ are greater than $1$ by our algorithm. Since the base concentration , this implies that off-target polymers are in smaller concentration.

While our formalism allows different concentration exponents for different polymers in $𝒮$, nearly all insight can be obtained from the simple case of uniform on-target concentration exponents:

Our main goal is to ensure that in equilibrium the concentration of on-target polymers is much higher than the others. This section presents two examples illustrating why this is nontrivial.

Consider a uniform set of on-target polymers $𝒮=\{A,B,C,D\}$, aiming at the equilibrium concentration $x_A=x_B=x_C=x_D=c$. For two off-target polymers $P_1$ and $P_2$, suppose that we figured out a canonical reaction

This reaction is in detailed balance if and only if $x_{P_1}\cdot x_{P_2}=x_A^2\cdot x_B\cdot x_C=c^4.$ At this point, it is unclear how to balance $x_{P_1}$ and $x_{P_2}$ without inspecting the other reactions: For example, the reaction $\alpha$ is balanced by assigning $x_{P_1}=x_{P_2}=c^2$, but another canonical reaction $A+2B+2D\to 2P_1$ would not be if it exists.

The following example exhibits a different potential problem. While Remark 3.2 shows that reactions entirely within $𝒮$ are balanced, it is not clear that given our choice of concentration exponents $\mu$ for $𝒮$, the reactions involving polymers outside of $𝒮$ could be balanced at all. Suppose there are two reactions

These reactions intuitively say that the off-target polymer $P_3$ is non-interacting to other polymers in these reactions so that the above problem of balancing the concentrations over off-target polymers is absent.⁴⁴4Looking ahead, these non-interacting polymers are indeed easier to assign concentrations to as shown in Section 5.1. However, balancing the reaction $\beta$ suggests the concentration $x_{P_3}=c^2$ at equilibrium, but the reaction $\gamma$ gives $x_{P_3}=c^{1.5}$. In other words, it is unclear if the concentrations of each of the polymers in $𝒮$ can be the same or even close to each other.

Despite these difficulties, our main result guarantees that the configuration on on-target polymers we demand is in equilibrium without the above issue, and can be extended to the configuration over all polymers at equilibrium, and ensures that the off-target polymers remain in very low concentration at equilibrium. Intuitively, we prove that there exists a canonical reaction for which we can assign concentrations to product polymers without creating conflicts with other reactions, and we provide a method to find this reaction. Moreover, as we will see later, these assignments occur in order of decreasing concentration allowing us to restrict the concentration of off-target species.

For the remainder, we fix ${𝚿}^{\mathrm{𝟎}}$, $𝚿$, and a particular set of on-target polymers $𝒮$ with concentration exponents $\mu$.

To systematically assign concentration exponents to all polymers, we will organize them into hierarchical groups called levels. The process begins with a designated set of polymers $𝒮$, the on-target polymers, which are assigned initial concentration exponents via $\mu :𝒮\to {ℝ}^{+}$.

All other polymers, called off-target polymers, will be partitioned into level sets ${𝒮}_1,{𝒮}_2,\mathrm{\dots }$, constructed inductively based on how these polymers appear in certain canonical reactions. At each level $i$, we will compute a scalar value ${\mu }_i$ that serves as the concentration exponent for all polymers newly added at that level.

In this way, we gradually extend the initial function $\mu$ to a global function $\overline{\mu }:𝚿\to {ℝ}^{+}$ that assigns a concentration exponent to every polymer in the system. This inductive process and the precise definitions of ${\mu }_i$, ${𝒮}_i$, and $\overline{\mu }$ will be described in detail below.

While $\widehat{{ℳ}_{\mathcal{2}}}$ technically depends on the level index $i$, we suppress this dependency in the notation for simplicity; it will be clear from the context that it is updated at each level. Additionally, by the definition of a canonical reaction, we have $\overline{\mu }({ℳ}_{\mathcal{1}})=\mu ({ℳ}_{\mathcal{1}})$ in the expression above.

At the heart of this inductive construction is the requirement that each $i$-levelizing reaction maintains balance with respect to the assigned concentration exponents. That is, for every reaction $\alpha :{ℳ}_{\mathcal{1}}\to {ℳ}_{\mathcal{2}}$ used to define level ${𝒮}_i$, the assigned exponents ensure the reaction remains detailed balanced, $c^{\overline{\mu }({ℳ}_{\mathcal{1}})}=c^{\overline{\mu }({ℳ}_{\mathcal{2}})}$. This holds because each polymer in ${ℳ}_{\mathcal{2}}$ that already appeared in previous levels contributes its already-defined exponent, while newly introduced polymers in ${𝒮}_i$ all receive the same value ${\mu }_i$. Decomposing ${ℳ}_{\mathcal{2}}$ into previously levelized polymers ${\widehat{ℳ}}_{\mathcal{2}}$ and new polymers in ${𝒮}_i$, we obtain:

ensuring that the total concentration exponent on both sides of the reaction is equal. This confirms that the assignment of ${\mu }_i$ for polymers in level ${𝒮}_i$ is consistent with the balance requirements dictated by the underlying reactions.

The full procedure for determining the level sets ${𝒮}_i$ and the corresponding concentration exponents ${\mu }_i$ is formally specified in Definitions 5.1 and 5.2 and carried out through the step-by-step process described in Algorithm 1. This algorithm inductively builds the extended concentration exponent function $\overline{\mu }$ by identifying polymers that can be levelized at each step and assigning them an appropriate exponent value based on their participation in canonical reactions.

While the algorithm itself does not explicitly state a stopping condition, its termination is ensured by the structure of the level construction process. Each time a new level ${𝒮}_i$ is defined, it includes at least one polymer not present in any previous level. Since the set of all polymers $𝚿$ is finite by assumption, only finitely many levels can be introduced. As a result, the inductive process must terminate after assigning a level and concentration exponent to each polymer, thus completing the definition of the extended function $\overline{\mu }$.

For us, canonical reactions represent the simplest meaningful interactions and serve as building blocks for more complex behavior. Their central role for our results is motivated by the following lemma:

Often we are only interested in an upper bound on the concentration of an off-target undesired polymer. In this case, rather than generating its exact equilibrium concentration via Algorithm 1, we can more efficiently compute an upper bound by narrowing our focus to reactions that directly produce $P$ instead of examining the full set of canonical reactions. This leads to a simpler surrogate quantity that approximates $\overline{\mu }(P)$ from below.

Since $P$ has not yet been levelized by construction, any such reaction $\alpha$ producing $P$, obviously involves at least one unlevelized polymer, ensuring $l_i(\alpha )\ne 0$.

Intuitively, ${\tilde{\mu }}_i(P)$ captures a conservative estimate of the exponent with which polymer $P$ can grow in concentration at level $i$, based only on reactions that actually produce $P$. Since it is defined as a minimum over a subset of possible reactions, it is easier to compute than the full concentration exponent $\overline{\mu }(P)$, yet is still meaningful for analysis:

While the TBN model is combinatorial in nature, quantifying over discrete (saturated) configurations, at the end we are most often interested in determining real-valued concentrations, which are accessible to (bulk) wet-lab experiment. The framework developed in this paper helps to bridge this gap between combinatorics of discrete configurations and concentrations.⁸⁸8Of course, much of the heavy lifting in bridging this gap is done by the derivation [4] of the cost function $g(𝐱)$, but our work expands on it beyond numerical simulation.

Consider a saturated (i.e., maximally bonded) configuration $ℳ$ in the TBN model. Quantified over all saturated reconfigurations ${ℳ}^{\prime }\cong ℳ$, the key quantity of interest in the TBN model is the “entropy” of ${ℳ}^{\prime }$, defined as the number of polymers in ${ℳ}^{\prime }$ (i.e., $|{ℳ}^{\prime }|$). The TBN model defines the “stable” configurations to be those that have maximum entropy among all saturated configurations.

Corresponding to the TBN literature [5], a multiset of polymers $ℳ$ is called TBN-stable if any reaction $ℳ\to {ℳ}^{\prime }$ has $|ℳ|\ge |{ℳ}^{\prime }|$. We say that a set $𝒮$ of polymers is TBN-stability closed if every multiset $ℳ\in {\mathbb{N}}^{𝒮}$ is TBN-stable and further any reaction $ℳ\to {ℳ}^{\prime }$ where ${ℳ}^{\prime }$ contains a polymer outside of $𝒮$ (i.e., $P\in {ℳ}^{\prime }$ for some $P\notin 𝒮$) satisfies $|ℳ|>|{ℳ}^{\prime }|$. In other words, producing a polymer outside of a TBN-stability closed set $𝒮$ costs entropy. The following lemma connects our notion of stability of on-target polymers (Definition 3.4) to TBN-stability.

In this section we show several applications of our mathematical tools in the analysis of existing systems of interest in the DNA molecular programming literature. We base our arguments on previous results proving the entropy loss (quantity $e(\alpha )$ of Corollary 7.2) of the systems in producing off-target (leak) species.

We note that while existing literature succeeds in characterizing entropy loss via TBN-like combinatorial arguments, more work is needed to develop similarly rigorous combinatorial arguments on novelty (quantity $l(\alpha )$); see also the discussion in Section 9. In the examples to follow, we claim to have identified the worst-case canonical reactions – i.e., having the least $e(\alpha )/l(\alpha )$ or $k(\alpha )/l(\alpha )$ ratio – without proof.

We first consider the TBN AND gate introduced in prior work [5, 2] and recently experimentally realized [14]. Figure 1 (left) shows the desired functionality of the AND gate in which inputs $A$ and $B$ cooperate to produce $C$. We are interested in bounding the concentration of $C$ (leak) that can be produced in the absence of inputs $A$ and $B$. Phrased in our terminology, combinatorial TBN arguments have shown that $𝒮=\{X,Y,Z\}$ is TBN-stability closed, as well as $\{X,Y,Z,A\}$ and $\{X,Y,Z,B\}$ where one of the two inputs is present. This implies that any canonical reaction $\alpha$ producing $C$ has entropy loss $e(\alpha )\ge 1$. However, such arguments do not directly connect entropy loss to leak concentrations, justifying the need for a tool like our Corollary 7.2.

To apply our framework to the TBN AND gate, we use Corollary 7.2. In the absence of both inputs, take $𝒮=\{X,Y,Z\}$. We claim that the worst-case canonical reaction $\alpha$ producing $C$ is $X+Y+Z\to G_1+C$ shown in Figure 1 (middle). This reaction has entropy loss $e(\alpha )=1$ and novelty $l(\alpha )=2$ since $G_1$ and $C$ are outside of $𝒮$. Thus for any base concentration , there is an equilibrium with $[X]=[Y]=[Z]=c$ where the leak concentration is $[C]\le c^{1.5}$.

Similarly, consider the case of having input $B$ but no input $A$. We take $𝒮=\{X,Y,Z,B\}$ and claim⁹⁹9While we do not provide a formal proof that this reaction is worst-case, our confidence is based on the observation that we cannot increase the novelty further without increasing the entropy loss in proportion. For example, we can combine two copies of reaction $\alpha$ to yield reaction ${\alpha }^{\prime }$ but then $e({\alpha }^{\prime })=2e(\alpha )$ and $l({\alpha }^{\prime })=2l(\alpha )$, maintaining the same ratio. that the worst-case canonical reaction $\alpha$ producing $C$ is $X+Y+Z+B\to G_2+G_3+C$ as shown in Figure 1 (right). This reaction has entropy loss $e(\alpha )=1$ and novelty $l(\alpha )=3$ since $G_2$, $G_3$, and $C$ are outside of $𝒮$, yielding the equilibrium leak concentration of $[C]\le c^{1.33}$. Our analysis thus provides concrete polynomial upper bounds on leak, consistent with the qualitative expectation of the TBN model that leak should become comparatively negligible in the limit of decreasing $c$. The bound is smaller when both inputs are absent than when $B$ is present.

As the next example, we consider the “leakless” DNA strand displacement system previously theoretically and experimentally studied [15]. That work focuses on a family of “translator” modules that convert an input strand to an output strand of independent sequence. The family is parameterized by the redundancy parameter $N$ defined as the number of bound domains in each fuel polymer $F_i$, such that the number of domains in each signal $X_i$ is $N+1$. Leak is expected to decrease with decreasing overall concentration (as for the AND gate), as well as with increasing redundancy $N$. The decrease of leak with increasing $N$ was confirmed by experiment, at least for small $N$. The reactions of the translator with $N=3$ are shown in Figure 2.

To apply our framework, we consider the system without toeholds, driven solely by entropy; with long domains alone, we are in the enthalpy neutral (athermic) regime. Prior work focused on the case where the system was prepared in a state with only fuel polymers ($F_i$), all at equal concentration, and zero initial concentration of waste polymers ($W_i$). Let on-target set be $𝒮=\{F_i,W_i\mid i=1,\mathrm{\dots },n\}$ with uniform $\mu =1$. Rephrasing in our terminology, ref. [13] proved that any canonical reaction $\alpha :ℳ\to {ℳ}^{\prime }$ where $X_i\in {ℳ}^{\prime }$ has entropy loss $e(\alpha )=N-1$, and this fact was used to argue that by increasing $N$ we can make leak arbitrarily small. We now show that considering novelty in addition to entropy loss makes this argument problematic, and suggest an alternative parameter setting to arbitrarily decrease leak.

Consider a cascade of two translators. Importantly, the number of fuels for a single translator module is $N+1$; i.e., $X_i$ and $X_j$ overlap in sequence for but are completely sequence-independent for $j=i+N+1$. We claim that the worst-case leak pathway $\alpha$ is where the first translator leaks resulting in the triggering of the second translator. This pathway generates $l(\alpha )=N+3$ new off-target polymers: $1$ for the leaked upstream translator (all fuels coming together), $N+1$ for the triggered fuels of the second translator, and $1$ for the output. Thus $e(\alpha )/l(\alpha )=(N-1)/(N+3)$. By Corollary 7.2, the concentration of the leak product in the uniform setting is bounded above by $c^{4/3}$ for $N=3$. The bound tightens to $c^2$ with increasing $N$; however, it does not get arbitrarily small. This suggests that we do not decrease equilibrium leak concentration arbitrarily by increasing “redundancy” $N$ because while the entropy loss increases, the novelty increases as well.

To arbitrarily decrease leak, we propose to use positive initial concentrations of the waste polymers $W_i$. The fuel polymers are denoted $F_1,\mathrm{\dots },F_n$, and the waste polymers produced after translator triggering are $W_1,\mathrm{\dots },W_n$. Let $X_0$ represent the input and $X_n$ represent the output. We have $n=2(N+1)$ for two translators composed together. This cascade of length $n$ proceeds through the following reactions described in Figure 2:

for $i=1,\mathrm{\dots },n$.

We first consider the triggered system (with-input) showing at least a constant fraction of the signal is propagated to the end as $X_n$. We then focus on the case without input and bound the leak. Note that this analysis utilizes Theorem 5.4 directly rather than Corollary 7.2 because we will have different concentration exponents $\mu$ for $F_i$ and $W_i$ in $𝒮$.

For the triggered (with-input) system we define our on-target set as $𝒮=\{F_i,W_i,X_0,X_i\mid i=1,\mathrm{\dots },n\}$. We assign concentrations to the on-target polymers as follows: All fuel concentrations are equal to $2c$, and all waste concentrations are equal to $c$. These concentrations correspond to concentration exponents $\mu (F_i)=1+{\log }_{c}2$ and $\mu (W_i)=1$. Let the concentration of the output in the final layer $X_n$ be $y$ and assign balancing concentrations of the other $X_i$. Since $[F_i]/[W_i]=2$, we have , meaning that more than half of the total signal ($X_i$) is at the output layer.

Now, we investigate the system without input. Let $𝒮=\{F_i,W_i\mid i=1,\mathrm{\dots },n\}$. For the situation to properly correspond to the with-input case, we need to ensure that all monomer concentrations are the same between the two cases, except removing the monomer corresponding to input $X_0$. Rather than thinking about specific monomers, we start in the with-input case and conceptually run reactions $X_i+W_i\to X_{i-1}+F_i$ to completely push all $X_i$ to $X_0$, and then remove $X_0$ from the system.¹⁰¹⁰10More precisely, this ensures that the concentrations of monomers making up on-target polymers are the same between the with-input and without-input case (other than $X_0$). Since the total amount of all $X_i$ is less than $2y$ in the with-input case, this results in: and $[W_i]>c-2y$. These correspond to $\mu (F_i)>{\log }_c(2c+2y)$ and .

Recall that redundancy $N$ results in entropy penalty $N-1$. We claim that the reaction with the smallest imbalance-novelty ratio (i.e., worst-case) is reaction $\beta$:

where $G$ is the “large polymer” formed after all $F_1,\mathrm{\dots },F_N$ displace the top strand from $F_{N+1}$, and $X_n$ is the leak output. The imbalance of this reaction is:

We can ensure that $k(\beta )$ is at least a constant fraction of $n$ for small enough $c$. For example, if we let $y\le c/4$, then $k(\beta )\ge n/4$ for any $c\le 0.0064$. The novelty is independent of $n$: $l(\beta )=2$ since $G$ and $X_n$ are not in $𝒮$. Therefore, the imbalance-novelty ratio $k(\beta )/l(\beta )$ of the worst case reaction is at least $n/8$, which increases linearly with $N$ (recall $n=2(N+1)$ for two translators composed together). Applying Theorem 5.4 leads to leak concentration of at most $c^{n/8}=c^{1/4+N/4}$. This upper bound¹¹¹¹11Note that this upper bound is loose because of inequalities such as Equation 4. implies smaller-than $c$ concentration of leak for $N\ge 4$, with the leak exponentially decreasing for larger $N$.

To summarize, by increasing redundancy $N$ in the appropriate regime, we maintain the property that a constant fraction of the input is converted to output in the with-input case, while arbitrarily (exponentially in $N$) decreasing leak in the without-input case.

Our results suggest a few important directions for future work. Given the central role of worst-case canonical reactions – i.e., canonical reactions with the lowest imbalance-novelty ratio (for Algorithm 1 and Theorem 5.4) or entropy loss-novelty ratio (Corollary 7.2) – it is important to develop formal techniques to prove that a given canonical reaction is indeed worst-case overall, at least or for a particular off-target polymer (for Section 6). Note that while combinatorial techniques in prior work in the TBN model have focused on proving entropy loss, more work is needed to study the ratio directly, making our framework more easily applicable. While we believe the canonical reactions highlighted in Section 8 are indeed worst-case, the argument is informal.

Another promising avenue of research is to establish a more direct link between a polymer’s monomer composition and its equilibrium concentration. Our current framework is effectively reaction-centric, inferring concentrations based on how polymers transform into one another. An alternative approach could be to derive concentration bounds directly from the structural properties of the off-target polymers, such as their size (monomer count) or degree of overlap with one another (multiset difference). Nonetheless, we hope that a variety of structure-based results could be proven based on a reduction to our canonical reaction framework.

Finally, this work has focused exclusively on the athermic case, where all molecular interactions are enthalpically neutral. While this is a reasonable and useful abstraction for systems with strong, saturated bonds, many real-world molecular systems, including many popular in DNA molecular programming, involve a range of binding strengths and enthalpic effects (e.g., from toehold binding). Extending our algorithmic framework to incorporate user-specified $\mathrm{\Delta }G$’s for each polymer could significantly broaden its applicability, and although this would complicate our algorithm, we do not anticipate any insurmountable difficulties.