{'fraction_non_alphanumeric': 0.061736301457591514, 'fraction_numerical': 0.02346491789840408, 'mean_word_length': 4.354272390618571, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 17, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Noise Control for DNA Computing

Synthetic biology is a growing interdisciplinary field, with far-reaching applications, which aims to design biochemical systems that behave in a desired manner. With the advancement of strand-displacement DNA computing, a large class of abstract biochemical networks may be physically realized using DNA molecules. Methods for systematic design of the abstract systems with prescribed behaviors have been predominantly developed at the (less-detailed) deterministic level. However, stochastic effects, neglected at the deterministic level, are increasingly found to play an important role in biochemistry. In such circumstances, methods for controlling the intrinsic noise in the system are necessary for a successful network design at the (more-detailed) stochastic level. To bridge the gap, the noise-control algorithm for designing biochemical networks is developed in this paper. The algorithm structurally modifies any given reaction network under mass-action kinetics, in such a way that (i) controllable state-dependent noise is introduced into the stochastic dynamics, while (ii) the deterministic dynamics are preserved. The capabilities of the algorithm are demonstrated on a production-decay reaction system, and on an exotic system displaying bistability. For the production-decay system, it is shown that the algorithm may be used to redesign the network to achieve noise-induced multistability. For the exotic system, the algorithm is used to redesign the network to control the stochastic switching, and achieve noise-induced oscillations.


Introduction
Synthetic biology is an interdisciplinary field of science and engineering that aims to construct biochemical systems with prescribed behaviors [1,2]. At the theoretical level, the synthetic systems may significantly enhance our understanding of biology. At the practical level, they may have broad applications, e.g. in medicine [3,4,5,6,7,8], industry [9,10], and nanotechnology [11,12]. The systems may also be of interest to NASA for optimizing extraterrestrial explorations [13]. A proof-of-concept for synthetic biology is a synthetic oscillator called the repressilator, which was implemented in vivo [14]. The experimental advances since the repressilator range from isolated synthetic biochemical networks, to microorganisms containing partially, or even fully, synthetic DNA molecules (synthetic life) [15,16,17,18]. Examples include microorganisms containing a synthetic bistable switch [19], and a cell-density controlling quorum sensor [20], microorganisms producing antimalarial drugs [7,8], and synthetic systems designed for tumor detection, diagnosis and adaptive drugresponse [4,5].
The construction of biochemical networks in synthetic biology may be broken down into two steps: firstly, an abstract system is constructed, displaying prescribed properties, and taking the form of a chemical reaction network [21,22,23]. Secondly, the abstract network is mapped to a suitable physical network, which may then be integrated into a desired environment (e.g. a test-tube, or a living cell) [24].
In the first step of network construction, the goal is to obtain an abstract network with desired dynamics. In this paper, we consider two dynamical models of reaction networks under mass-action kinetics [25,23]: the deterministic model, and the stochastic model (see Methods for more details). The deterministic model takes the form of the reaction rate equations, which are ordinary-differential equations governing the time-evolution of the species concentrations [23,25]. The stochastic model takes the form of a Markov chain, which may be simulated using the Gillespie stochastic simulation algorithm [26]. The Gillespie algorithm generates noisy copy-number time-series, with the copy-number distribution matching that obtained from the underlying chemical master equation [23,25,26,27]. The stochastic model is more-detailed, taking into an account the discreteness of the species counts, and the stochastic nature of the dynamics, which may be particularly important in biochemistry, where reaction networks may contain lowabundance species [31,32,14,19,22,28,29,30]. On the other hand, the deterministic model is less-detailed, and more appropriate when the species are in high-abundance, and the discreteness and stochasticity are negligible [33].
In the second step of network construction, the goal is to engineer a physical network whose dynamics match well the dynamics of a given abstract network, over a suitable timeinterval. Engineering an appropriate physical network may proceed indirectly, by altering a preexisting physical network, or directly, by engineering a network, which involves a given set of physical species, from scratch. The advantage of the former approach is that a preexisting network may display (partially) desirable dynamical properties. However, such a network may involve DNA and RNA molecules, proteins, and metabolites [2], some of which may have complex biophysical properties. Consequently, the disadvantage is that the structure (and, thus, the dynamics) of such a network cannot generally be modified in an arbitrary manner. In the latter approach, one may choose the physical species, at the expense of having to build a network from scratch. In the subfield of DNA computing, the latter approach is followed, and physical networks are engineered with chemical species consisting exclusively of DNA molecules, interacting via the toeholdmediated DNA strand-displacement mechanism [24]. DNA production is systematic and cost-effective, and, due to the fact that DNA biophysics is relatively well-understood, one has more freedom in controlling the structure of corresponding physical networks. More precisely, an abstract network under mass-action kinetics may be mapped to a DNA-based physical network provided it consists of up to second-order reactions, with rate coefficients varying over up to six orders of magnitude. The resulting physical network has identical deterministic dynamics as the abstract network (in the asymptotic limit of some of the kinetic parameters [24]), up to a scaling of the dependent variables. A proof-of-concept for DNA computing is a synthetic oscillator called the displacillator, which was implemented in vitro [34].
While the deterministic model of reaction networks is lessdetailed, it is also simpler than the stochastic model, making it attractive for guiding the construction of networks, predicting accurately their mean-field behavior [24,14,19,21,22,23]. However, when noise is an important part of the dynamics, the stochastic model has to be considered. The intrinsic noise, often arising in biochemistry, may be controlled in two ways: it may be decreased (e.g. as in [32]), in order to reduce the differences between the stochastic and deterministic dynamics. On the other hand, it may be increased, in a state-dependent manner, in order to favorably change the stochastic dynamics. In the language of molecular computing, the latter approach corresponds to exploiting the proven computational power of the stochastic reaction networks [35], by reprogramming the underlying intrinsic noise. Let us note that exploitations of the noise for enhancing biological functions have been reported in applications [31,30]. In this paper, we follow the latter approach, and present the noise-control algorithm (given as Algorithm 1) which maps an input reaction network to output networks whose stochastic dynamics have an additional controllable state-dependent noise. Importantly, the input and output networks have identical deterministic model in appropriate limits of some of the parameters introduced by the algorithm. The algorithm may play a significant role in the biochemical network synthesis, allowing for a deterministic-stochastic hybrid approach. More precisely, when constructing abstract and physical networks, one may use the deterministic model to guide the construction, and then apply the algorithm to favorably modify the intrinsic noise in the stochastic model, while preserving the desired deterministic dynamics. The algorithm may also be used to adjust the intrinsic noise to favorably interact with environment-induced effects (e.g. extrinsic noise).
The rest of the paper is organized as follows: we introduce Algorithm 1 by applying it to the test network (1), which at the deterministic level displays a globally attracting equilibrium point. We show that the algorithm can favorably modify the stationary probability distribution underlying (1) at arbitrary points of the state-space, without influencing the deterministic dynamics. For example, it is shown that the algorithm may be used to redesign (1) to achieve noise-induced multimodality (multistability). We then apply Algorithm 1 to the exotic network (11), which at the deterministic level displays a bistability involving an equilibrium point and a limit cycle. The algorithm is used to redesign (11) to increase the stochastic switching between the two attractors, and to achieve noise-induced oscillations.

A One-species Regular System
Consider the one-species production-decay reaction network R(s), given by (1).
(2) Species s from network (1) reacts according to the two reactions with rate coefficients k 1 , k 2 ∈ R ≥ , where R ≥ is the set of nonnegative real numbers, and ∅ is the zero-species (denoting species which are not of interest). In this paper, we assume reaction networks are under mass-action kinetics, with the reactions taking place in unit-volume reactors. Let us denote the concentration of species s from (1) at time t ∈ R ≥ byx =x(t) ∈ R ≥ . The initial value problem for the deterministic model (also called the drift) for network (1) is given by system (2), withx 0 ≥ 0 (see also Methods). Since the deterministic model (2) has a globally attracting equilibrium point, given by k 1 /k 2 , network (1) is said to be regular [23].
Let us denote the copy-number of species s from (1) at time t ≥ 0 byX(t) ∈ N 0 , where N 0 is the set of integers. Under the stochastic model,X(t) is modelled as a continuoustime, discrete-space Markov chain (see also Methods), which can be generated by using the Gillespie stochastic simulation algorithm [26]. GivenX(t), there will be a mean interevent time until one of the reactions from (1) fires. The mean interevent time is given by 1/α(X(t)), and when the event takes place, the probability that the i-th reaction from (1) fires is equal toα i (X(t))/α(X(t)), for i ∈ {1, 2}. Here,α 1 = k 1 , andα 2 (x) = k 2 x, are the so-called propensity functions of the first, and second, reactions from (1), respectively. Function α(x) = k 1 + k 2 x is the total propensity function of network (1), i.e. the sum of propensity functions of all the underlying reactions.
We now wish to structurally modify network (1) in such a way that the deterministic model from (2) is preserved, while an arbitrary nonnegative function, defined on a bounded discrete domain, is added to the total propensity function of (1). The latter requirement implies that the interevent time would be controllably decreased in a state-dependent manner. Equivalently, the two requirements imply that a controllable state-dependent noise would be introduced into the stochastic dynamics. We have designed a three-step algorithm, given as Algorithm 1, which achieves such goals for arbitrary reaction networks under mass-action kinetics. Let us describe properties of the algorithm by applying it on network (1).
Firstly, we wish to introduce an additional speciess into network (1), in such a way that species s ands satisfy a pairwise stoichiometric conservation law. Secondly, we require that the enlarged network has the same deterministic model as network (1), despite the added speciess, which may be achieved by adding another auxiliary species. More precisely, let us consider networkR 1 (s,s) ∪ R 2 1 (s), given by: Species s,s, I 1 from (3) react according to the four reactions with rate coefficients k 1 , k 2 , 1/µ ∈ R ≥ . NetworkR 1 =R 1 (s,s), given in (3), is obtained from networkR =R(s), given by (1), in the following way: since the first reaction inR increases copy-number of s by one,s and I 1 are added to the reactants of the reaction, and I 1 is added to the products, leading to the first reaction inR 1 . Since the second reaction inR decreases copy-number of s by one,s is added to the products, leading to the second reaction inR 1 . This ensures that the desired conservation law holds. The superscript in I 1 indicates that species I 1 is involved as a catalyst in a reaction ofR 1 in which s is increased by one. The subscript in R 2 1 = R 2 1 (s) indicates that the network describes production and decay of I 1 .
The initial value problem for the deterministic model of (3) is given by where x = x(t) ∈ [0, c] ∩ R ≥ , and y = y(t) ∈ R ≥ , are the concentrations of species s, and I 1 , from (3), respectively, with x 0 , y 0 , c ∈ R ≥ . We have used the kinetic conservation law x(t) = c − x(t), wherex(t) is the concentration of speciess, and c < ∞ is a time-independent conservation constant. Note that the conservation law truncates x-state-space. Let us now describe relationships between systems (2) and (4), starting with the weak statement: for c > k 1 /k 2 , and for any µ ≥ 0, solutions of (2) and (4) are the same in the long-time limit t → ∞. More precisely, the x-component of the equilibrium point of (4) is identical to the equilibrium point of (2), and both are stable. In Supplementary Information (SI) Text, we justify the strong statement: for sufficiently large c, and for µ 1, solutions of (2) and (4), with the same initial conditions, are approximately the same at each time t ≥ 0. For these reasons, we call R 2 1 a drift-corrector network.

Zero-Drift Network
Having completed the first two steps, let us focus on the third (and final) step, in which we introduce arbitrary noise into the stochastic model of (3), without influencing the deterministic model (4). Let us start our consideration by embedding into (3) network R 3 1,1 = R 3 1,1 (s,s), which is given by The subscript in R 3 1,1 indicates that the underlying reactions have one molecule of s, and one ofs, as reactants. The two reactions in (5) preserve the conservation law from (3). Furthermore, they fire with the same rates, with the first reaction leading to a unit-production, while the second to a unit-decay, of species s. Consequently, embedding R 3 1,1 into (3) does not affect the underlying deterministic model (4), and we call R 3 1,1 a zero-drift network. However, R 3 1,1 does affect the underlying stochastic model [36,37,38,23]. To illustrate this, let us consider network R 3 1,1 in isolation: the reactions from (5) fire when X(t) ∈ (0,C), but not when X(t) ∈ {0,C}, so that R 3 1,1 in isolation fires until X(t) takes one of the extreme values {0,C}. Here, X(t) ∈ N 0 , and C ∈ N, C < ∞, are the copynumber of species s appearing in (3) and (5) at time t ≥ 0, and the conservation constant, respectively. Let us note that a possible biologically-relevant realization of network (5), aside from DNA strand-displacement mechanism, is a dimer version of the bifunctional histidine kinase/phosphatase reported in [39].
In SI Text, we derive equation (SI7) which describes the effective behavior of the Markov chain X(t) from network R 1 ∪ R 2 1 ∪ R 3 1,1 in the limit µ → 0, and it follows that the effective total propensity function of the network, denoted α(x), satisfies Functionα(x) has the form of the total propensity of network (1), and K 1,1 β 1,1 (x) is the propensity function of reactions in (5), with the scaled factors given by Figure 1(a), where one can notice its parabolic shape, arising from the underlying conservation law X(t) +X(t) = C, which holds for all t ≥ 0, wherē X(t) ∈ N 0 is the copy-number ofs at time t ≥ 0. Comparing (6) and (7), it follows that, as µ → 0, the mean interevent time for X(t), from networkR 1 ∪R 2 1 ∪R 3 1,1 , is lower than that ofX(t), from network (1), in the regions of the common statespace where β 1,1 (x) = 0, i.e. for x ∈ (0,C). Coefficient K 1,1 controls by how much the interevent time is reduced. Equivalently, β 1,1 (x), and K 1,1 , determine the support, and magnitude, respectively, of the state-dependent intrinsic noise which network (5) introduces into the dynamics of network (3).
To study this further, in SI Text we derive the following two equations (given as (SI9), and (SI13), respectively) where p(x) is the stationary probability mass function (PMF) corresponding to networkR 1 ∪ R 2 1 ∪ R 3 1,1 in the limit µ → 0, i.e. the probability that there are x molecules of species s as µ → 0 in the long-time limit t → ∞. Let us interpret analytical results (9) and (10), and compare them with the numerically obtained counterparts. In Figure 1(b), we display numerically obtained stationary x-marginal PMFs for different values of K 1,1 , with the rest of the (dimensionless) parameters fixed to k 1 = 2.5, k 2 = 0.5, µ = 10 −3 , and C = 15. It can be seen that, for K 1,1 = 0, i.e. when the zero-drift network R 3 1,1 does not fire, the PMF matches that of network (1), i.e. it is a Poissonian, as predicted by (9). Let us note that the matching of the PMFs of networks (1) andR 1 ∪ R 2 1 ∪ R 3 1,1 relies on choosing sufficiently large rate coefficients 1/µ in the drift-corrector network R 2 1 . When K 1,1 = 5, the PMF appears closer to a uniform distribution, than does the PMF when K 1,1 = 0. Finally, for the larger value K 1,1 = 10 5 , i.e. when zero-drift network R 3 1,1 fires much faster than networkR 1 , the PMF redistributes across the domain, accumulating at the boundary, and becoming bimodal. This is in qualitative agreement with (6), and in quantitative agreement with (10), which predicts p(0) ≈ 0.7 and p(15) ≈ 0.3. In Figure 1(c), a representative sample path is shown, obtained by applying the Gillespie algorithm on networkR 1 ∪ R 2 1 ∪ R 3 1,1 , when K 1,1 = 10 5 . Also shown is a trajectory obtained by numerically solving the deterministic model (4). Consistent with Figure 1(b), the sample path switches between the boundary of the state-space, with a bias towards the left boundary point x = 0. This is in contrast to the deterministic trajectories, which are globally attracted to the equilibrium point x = 5.

General Zero-Drift Networks R 3 n,n
Zero-drift network R 3 1,1 (s,s), given by (5), involves a single molecule of s ands as reactants, and adds the noise at x ∈ [1,C − 1], i.e. in the interior of the state-space. Similar networks may be used to add the noise at any point in the state-space, without influencing the deterministic dynamics. In particular, in (19) and (20), we present general zero-drift networks R 3 n,n (s,s), which involve n molecules of s, andn ofs, as reactants, and add the noise at x ∈ [n,C −n], where n,n ∈ N 0 , and (n +n) ≤ C (see also SI Text). Embedding a union of such networks, ∪ (n,n) R 3 n,n (s,s), into (3), we arrive at the result similar to (6), with K 1,1 β 1,1 (x) replaced by the linear combination ∑ (n,n) K n,n β n,n (x). The scaled rate coefficient K n,n , and function β n,n (x), are given as (S14), and (S15), respectively, in SI Text, where we also justify that an arbitrary nonnegative function, with compact support, may be approximated by a suitable sum ∑ (n,n) K n,n β n,n (x). To illustrate general zero-drift networks, let us start with embedding into network (3) zero-drift network R 3 5,10 (s,s), satisfying (19) with n = 5 andn = 10. In Figure 1(d), we show propensity function β 5,10 (x), which is nonzero only at x = 5. In (e), we show the numerically approximated stationary x-marginal PMFs underlying networkR 1 ∪ R 2 1 ∪ R 3 5,10 for different values of K 5,10 , with the rest of the parameters as in Figure 1(b). One can notice that, under the action of network R 3 5,10 , the PMF is gradually decreased to nearly zero at x = 5 (the deterministic equilibrium), and becomes bimodal, with the two noise-induced maxima at x = 4 and x = 6. In (f), we show a corresponding representative sample path.
In general, noise-induced multimodality may be achieved by a suitable combination of zero-drift networks. For example, let us synthetize noise such that the stationary PMF is trimodal, and nearly zero everywhere, except at x ∈ {1, 7, 11}. Such a task may always be achieved by a suitable combination of the basis zero-drift networks, i.e. those zero-networks that induce noise only at a single point in the state-space (e.g. subnetwork R 3 5,10 with propensity function shown in Figure 1(d), see also SI Text). In the present case, one could construct the thirteen basis zero-drift networks which add large enough noise at x ∈ [0, 15] \ {1, 7, 11}. Here, for simplicity, we achieve the task with only four zero-drift networks. In Figures 1(g) . We denote β (x) ≡ β 0,15 (x) + β 2,9 (x) + β 8,5 (x) + β 12,0 (x), and, for simplicity, take K ≡ K 0,15 = K 2,9 = K 8,5 = K 12,0 . The resultant propensity function β (x) is shown in (g), while in (h) it can be seen that the PMF becomes trimodal for sufficiently large K, with the maxima at x = {1, 7, 11}. This is consistent with the corresponding representative sample path shown in blue in panel (i), which display tristability. Let us note that, while the stochastic dynamics display multistability in (c), (f) and (i), the corresponding deterministic dynamics, also shown in the plots, remain monostable.

A Two-species Exotic System
Consider the two-species networkR(s 1 , s 2 ), given bỹ where species s 1 and s 2 react according to the eleven reactions with rate coefficients k 1 , k 2 , . . . , k 11 ≥ 0. We denote the copynumbers of species s 1 , and s 2 , at time t by X 1 (t), and X 2 (t), respectively. It was established in [21] that, for particular choices of the rate coefficients, the deterministic model of reaction network (11), given as equation (SI17) in SI Text, exhibits exotic dynamics: it undergoes a homoclinic bifurcation, and displays a bistability involving a limit cycle and an equilibrium point. On the other hand, it is demonstrated in [22] that the stochastic model of (11) is not necessarily sensitive to the deterministic bifurcation, and may effectively behave in a monostable manner. The latter point is demonstrated in Figure 2(c), where we show in red numerically approximated x 1 -solutions of (SI17), one initiated in the region of attraction of the equilibrium point, while the other of the limit cycle. For a comparison, we also show in blue a representative sample path generated by applying the Gillespie algorithm on (11). It can be seen that the stochastic solution spends significantly more time near the deterministic equilibrium point. To gain a clearer picture, we display in Figures 2(a), and (b), the joint, and the x 1 -marginal, stationary PMFs, respectively, underlying network (11), which have been obtained numerically for the same parameter values as in Figure 2(c). In (b), one can notice that the PMF is bimodal, but the left peak, corresponding to the limit cycle, is significantly smaller than the right peak, which corresponds to the stable equilibrium point. We now apply Algorithm 1 on network (11) to achieve two goals. Firstly, we balance the sizes of the two peaks of the stationary PMF from Figure 2(b), thereby forcing the stochastic system to spend comparable amounts of time at the two deterministic attractors. Secondly, we reverse the situation shown in Figure 2(b), by making the left PMF peak significantly larger than the right one, thereby forcing the stochastic system to spend most of the time near the limit cycle. We could achieve the goals by introducing speciess 1 ,s 2 into (11), and using suitable basis zero-drift networks. We take a simpler approach, by mapping (11) toR 1 (s 1 , s 2 ,s 2 ) ∪ R 2 1 (s 2 ) ∪ (R 3 0,C 2 −10 (s 2 ,s 2 ) ∪ R 3 30,0 (s 2 ,s 2 )), which is given by equation (SI18) in SI Text. For our purposes, only one ofs 1 ,s 2 is sufficient, since the stochastic dynamics of s 1 and s 2 are coupled. We have chosens 2 for convenience, since x 2 -state-space may be truncated at a lower value, C 2 = 180, than x 1 -state-space (see also Figure 2 (a)). The x 2 -component of the deterministic limit cycle satisfies x 2 ∈ (10, 30). Correspondingly, we introduce two zero-drift networks: R 3 0,C 2 −10 (s 2 ,s 2 ), and R 3 30,0 (s 2 ,s 2 ), which redistribute the PMF from x 2 ∈ [0, 10], and from x 2 ∈ [30,C 2 ], respectively, to the limit cycle region, x 2 ∈ (10, 30). We fix the scaled rate coefficient K 2 0,C 2 −10 to a large value (so that the PMF is nearly zero for x 2 ∈ [0, 10]), and vary the coefficient K 2 30,0 , which redistributes the PMF from the deterministic equilibrium point to the limit cycle. Network R 2 1 (s 2 ) is necessary for the preservation of the deterministic dynamics of (11) under the application of Algorithm 1.
In Figures 2(d), and (e), we show the joint, and x 1 -marginal, stationary PMFs for an intermediate value of K 2 30,0 , when the PMF is partially redistributed from x 2 ∈ [30,C 2 ] to x 2 ∈ (10, 30), so that the two peaks in (e) are of comparable sizes.
In Figure 2(f), we show a representative sample path, obtained by applying the Gillespie algorithm on network (SI18) from SI Text, together with the deterministic trajectories obtained by solving (SI17). One can notice that the stochastic system now spends significantly more time near the limit cycle, when compared to (c). In Figures 2(f)-(g), we show analogous plots, but for a sufficiently large value of K 2 30,0 , when the PMF is almost completely redistributed from x 2 ∈ [30,C 2 ] to x 2 ∈ (10, 30). Now, in contrast to Figures 2(a)-(c), the PMF becomes essentially unimodal, and concentrated around the limit cycle. Let us note that the red trajectories from Figures 2(f) and (i) were generated by numerically solving the deterministic model of network (11), given by (SI17). For our purposes, it is not necessary to solve the corresponding (stiff) deterministic model of network (SI18). The reason is that Algorithm 1 does not influence the deterministic equilibrium points of a given reaction network, regardless of the choice of the kinetic algorithm parameters. For example, while the deterministic limit cycle is not necessarily preserved for the algorithm parameters chosen in Figure 2(i), the enclosed deterministic unstable focus is necessarily preserved. Thus, the blue sample path corresponds to noise-induced oscillations either near a deterministic limit cycle, or near a deterministic unstable focus.

Summary
In this paper, we have presented the noise-control algorithm, which is given as Algorithm 1. The algorithm maps an input chemical reaction network to output networks, all under massaction kinetics, by introducing appropriate additional species and reactions, such that the output networks satisfy the following two properties. Firstly, the output networks have the same deterministic model as the input network, in appropriate limits of some of the parameters (rate coefficients) introduced by the algorithm. Secondly, controllable state-dependent noise is introduced into the stochastic model of the output networks. Thus, Algorithm 1 may be used to control the intrinsic noise of a given reaction network under mass-action kinetics, while preserving the deterministic dynamics. Let us note that the asymptotic conditions for the algorithm parameters are necessary for preservation of the time-dependent deterministic solutions. However, the time-independent deterministic solutions (the deterministic equilibrium points), which capture important features of the deterministic dynamics, are preserved under the algorithm even if the asymptotic conditions are not satisfied.
The algorithm has been applied to a test problem, taking the form of the one-species production-decay system given by (1). Using analytical and numerical methods, we have shown that the additional intrinsic noise, introduced by the algorithm, may be used to favorably modify the stationary probability mass function at arbitrary points in the state-space, as demonstrated in Figure 1. For example, in Figure 1(b), the noise is added to the whole interior of the state-space, while in (e) only at a single point, in both cases resulting in noise-induce bimodality. On the other hand, in Figure 1(h), by adding the noise to specific points in the state-space, the network is redesigned to display noise-induced trimodality. As shown in Figures 1(c), (f), (i), the blue stochastic trajectories display multistability, while the red deterministic ones remain monostable.
The algorithm has also been applied to a more challenging problem, taking the form of the two-species system given by (11), which, for the parameters taken in this paper, at the deterministic level displays a bistability involving an equilibrium point and a limit cycle [21,22]. At the stochastic level, the system is significantly more likely to be found near the equilibrium point, as demonstrated in Figures 2(a)-(c). We have used the algorithm to redesign network (11), so that the stochastic system spends comparable amounts of time near the two attractors, as demonstrated in Figures 2(d)-(f). The network was also redesigned to display noise-induced oscillations, which is shown in Figures 2(g)-(i).
The controllable state-dependent noise is generated by Algorithm 1 using the zero-drift networks (19) and (20). Any nonnegative function, defined on a bounded discrete domain, may be represented by a linear combination of propensity functions induced by an appropriate union of the zero-drift networks. Thus, choosing suitable zero-drift networks, the algorithm may control the intrinsic noise at arbitrary points in the state-space of the stochastic dynamics of reaction networks. The cost of such a precision in nose-control is a larger number of reactants in the underlying zero-drift networks. However, while the high-molecular reactions introduced by the algorithm are more expensive to synthetize, they do not limit applicability of Algorithm 1 to synthetic biology. The reason for this is that such reactions may always be broken down into sets of up-to bi-molecular reactions, with asymptotically equivalent deterministic and stochastic dynamics [40,41]. In particular, a zero-drift network, involving reactions of order (n +n), may be broken down into 2(n +n) − 2 reactions of upto second-order, which may be readily mapped to DNA-based physical networks. Algorithm 1 may constitute a qualitatively novel finding which will facilitate the progress of DNA computing [24]. In particular, a hybrid approach for constructing DNA-based reaction networks may be used: the deterministic model may be used to guide the construction of reaction networks, and then Algorithm 1 may be applied to favorably reprogram the intrinsic noise in the stochastic model, while preserving the mean-field behavior. The algorithm may be of critical importance when the synthetic networks involve species at low copy-numbers, since then the stochastic effects may play a significant role [31,32,14,19,22,28,29,30], uncontrollably contaminating the performance of the synthetic networks. In such circumstances, Algorithm 1 may be used for controlling the stochastic effects, enriching the DNA-based synthetic systems with novel, noise-induced functionalities.

Methods
Let us consider the mass-action reaction network R given by where s 1 , . . . , s N are the reacting species, k j the reaction rate coefficients, and c i j , c i j the stoichiometric coefficients. Let us denote by c j , c j ∈ N N 0 the vectors of the stoichiometric coefficients of reaction j, and ∆x j = c j − c j . The deterministic model of reaction network (12) is given by the following system of ordinary-differential equations (ODEs), also known as the reaction rate equations [23,25]: Here, x = x(t) ∈ R N ≥ is the vector of species concentrations, i.e. x i (t) is the concentration of species s i at time t, and x c j ≡ ∏ N l=1 x c l j l , with the convention that 0 0 ≡ 1. The stochastic model of reaction network (12) is given by the following system of difference-differential equations, also known as the chemical master equation (CME) [23,25,27]: Here, p(x,t) is the probability mass function (PMF), i.e. the probability that the vector of copy-numbers X = X(t) ∈ N N 0 of species s 1 , . . . , s N at time t is given by x. Linear operator L is called the forward operator, and step operator E is the propensity function [23,25] of the j-th reaction from (12), and is given by where x c l j l denotes a falling factorial of x l , i.e. x c l j l ≡ x l (x l − 1) . . . (x l − c l j + 1).

Supplementary Information (SI) Text
The Deterministic Dynamics of NetworkR 1 ∪ R 2 1 in the Limit µ → 0 Let us analyse system (4) in the asymptotic limit µ → 0. It follows from the Tikhonov theorem [42] that the ODE for y, given by second equation in (4), reduces to the algebraic equation y = (c − x) −1 as µ → 0. Substituting the algebraic equation into (4) results in Initial value problems (2) and (SI1) have the same form, and let us denote their solutions byx(t;x 0 ) and x(t; x 0 ), respectively. Then, choosing c ≥ max t≥0x (t;x 0 ) < ∞, and x 0 =x 0 , ensures that concentration of auxiliary speciess is nonnegative,x(t) = c − x(t) ≥ 0, and that the solutions of (2) and (4) are asymptotically equivalent in the limit µ → 0.
with i ≥ 2. Substituting (SI4) into (SI2), and equating terms of equal powers in µ, the following system of equations is obtained: Order 1/µ equation. A suitable form of the zero-order approximation of the PMF follows from the Bayes theorem: p 0 (x, y,t) = p 0 (y|x)p 0 (x,t), where p 0 (y|x) is the stationary PMF of y conditional on x, while p 0 (x,t) is the marginal PMF of x. Substituting p 0 (x, y,t) = p 0 (y|x)p 0 (x,t) into the first equation in (SI5), with t, x fixed, leads to −L 2 1 p 0 (y|x) = 0. It follows that p 0 (y|x) is a Poisson distribution with parameter (C − x) −1 , so that the zero-order PMF is given by Order 1 equation. Substituting (SI6) into the second equation in (SI5), summing over all the possible states y ∈ N 0 , using (SI3), and equalities ∑ y yp 0 (y|x) = (C − x) −1 and ∑ y p 0 (y|x) = 1, one obtains the effective CME, given by where L is the forward operator corresponding to network (1), and has the following form Setting the left-hand side (LHS) to zero, and taking K 1,1 = 0 in (SI7), and assuming C is fixed to a sufficiently large value, it follows that the stationary PMF is a Poisson distribution with parameter k 1 /k 2 [27]: otherwise. (SI9) Let us substitute the power-series expansion with i ≥ 2, into (SI7) with the LHS set to zero, and consider the limit K 1,1 → ∞. Then, equating terms of equal powers in 1/K 1,1 , one obtains: Order 1 equation. The solution to the first equation in (SI11) is given by where a ∈ R ≥ is an arbitrary constant. Order 1/K 1,1 equation. Multiplying the second equation in (SI11) by x, and summing over x ∈ N 0 , with the convention that f 0 (x) = 0 and β 1,1 (x) = 0 for x / ∈ [0,C], one obtains the solvability condition 0 = ∑ ∞ x=0 xL f 0 (x), which implies a = k 1 /k 2 . Substituting a into (SI12) leads to the zero-order approximation of the stationary PMF: The propensity function of reactions underlying R 3 n,n (s,s), n,n ∈ N 0 , and (n +n) ≤ C, is given by K n,n β n,n : [0,C] → R ≥ , with K n,n = M n,n k n,n , and β n,n (x) = (M n,n ) −1 where the scaling factor M n,n is introduced to approximately normalize β n,n (x), and is given by Here, we take the convention ∏ N i=0 f (i) = 1 if N < 0, where f (i) is an arbitrary function of i. Function β n,n (x) is nonzero on the interval [n,C −n], with the single maximum approximately at Cn/(n +n).
Interior zero-drift networks. Zero-drift network R 3 n,n (s,s), with n,n = 0, satisfies (19), and the propensity function of its reactions, which is proportional to (SI15), is nonzero only in the interior of the state-space. Since the propensity function of R 3 n,n (s,s), with n,n = 0, attains its maximum in the interior of the domain, we call the network an interior zero-drift network.
Boundary zero-drift networks. Network R 3 0,n (s,s), satisfying (20), is a zero-drift network in the limit µ 0,n → 0. Furthermore, in the same limit, the first two reactions from (20) have the same propensity function, which is proportional to (SI15) with n = 0, and which is nonzero at the left boundary point, x = 0. Similarly, network R 3 n,0 = R 3 0,n (s, s;B, k n,0 , µ n,0 ) is a zero-drift network as µ n,0 → 0, and its first two reactions have the same propensity function, which is nonzero at the right boundary point, x = C. Since networks with n = 0 (respectively,n = 0) generate propensity functions with the maximum values at the left (respectively, right) boundary point, we call such networks left (respectively, right) boundary zero-drift networks.
Basis zero-drift networks. Stoichiometric coefficients n,n control the support of the intrinsic noise, which network R 3 n,n introduces into the stochastic dynamics, via the control of support of the compact function (SI15). The larger the sum (n +n) is, with (n +n) ≤ C, the smaller the support of (SI15), and, hence, one obtains a more precise noise-control. In the special case when n +n = C, the propensity function (SI15) is nonzero only at a single point in the state-space, x = n. We call networks R 3 n,n (s,s), with n +n = C, basis zero-drift networks, and the corresponding propensity functions basis propensity functions. Any nonnegative function, defined on a bounded discrete domain, may be represented by a suitable linear combination of the basis propensity functions.

The Deterministic Model for NetworkR
The deterministic model of network (11) is given by where x 1 = x 1 (t), x 2 = x 2 (t) are the concentrations of species s 1 , s 2 , respectively, at time t.
where s 1 , . . . , s N , are the species, k j the reaction rate coefficients, and c i j , c i j the stoichiometric coefficients. (1) Step: Reaction networkR, given by (16), is mapped to a pairwise conservative networkR 1 given bŷ R 1 (s 1 , . . . , s N ,s 1 , . . . ,s N ) : Here,s i , I ∆x i j i are additional species, ∆x i j = (c i j − c i j ), and 1 N (·) is the indicator function of the natural numbers. (2) Step: For each species I ∆x i j i , a drift-corrector network is constructed, R 2 ∆x i j (s i ) = R 2 ∆x i j (s i ; I ∆x i j i , µ), given by where 0 ≤ µ 1. (3) Step: For each speciess i , a union of zero-drift networks may be constructed. Let n,n ∈ N 0 , and (n +n) ≤ C i . Network R 3 n,n (s i ,s i ) = R 3 n,n (s i ,s i ; k i n,n ), with n,n = 0, is given by Network R 3 0,n (s i ,s i ) = R 3 0,n (s i ,s i ; B i , k i 0,n , µ 0,n ), withn = 0, is given by where 0 ≤ µ 0,n 1, and B i is an additional species. Network R 3 n,0 = R 3 0,n (s i , s i ;B i , k i n,0 , µ n,0 ).
Output: An output reaction network R is given by where R 2 = ∪ i ∪ ∆x i j R 2 ∆x i j (s i ), and R 3 = ∪ i ∪ (n,n) R 3 n,n (s i ,s i ).