Quasi-robust control of biochemical reaction networks via stochastic morphing

One of the main objectives of synthetic biology is the development of molecular controllers that can manipulate the dynamics of a given biochemical network that is at most partially known. When integrated into smaller compartments, such as living or synthetic cells, controllers have to be calibrated to factor in the intrinsic noise. In this context, biochemical controllers put forward in the literature have focused on manipulating the mean (first moment) and reducing the variance (second moment) of the target molecular species. However, many critical biochemical processes are realized via higher-order moments, particularly the number and configuration of the probability distribution modes (maxima). To bridge the gap, we put forward the stochastic morpher controller that can, under suitable timescale separations, morph the probability distribution of the target molecular species into a predefined form. The morphing can be performed at a lower-resolution, allowing one to achieve desired multi-modality/multi-stability, and at a higher-resolution, allowing one to achieve arbitrary probability distributions. Properties of the controller, such as robustness and convergence, are rigorously established, and demonstrated on various examples. Also proposed is a blueprint for an experimental implementation of stochastic morpher.

Manipulating the dynamics of cells in vivo, and designing their synthetic counterparts in vitro, are complementary goals of synthetic biology, both involving overcoming nonlinear, non-modular and stochastic nature of biochemical networks [2]. In particular, when biochemical systems are integrated into smaller-volume compartments, such as living or synthetic cells, the lower copy-numbers of some of the underlying species give rise to intrinsic noise [9][10][11][22][23][24][25][26]. The induced stochasticity requires theoretical and experimental methods (see section S1 in the electronic supplementary material and §5, respectively) which are more involved than their larger-volume (deterministic) counterparts [14,15,27]. In this context, a central problem is the development of a controller network that, when embedded into a given input network, ensures that the resulting output network executes a predefined stochastic dynamics in a stable and accurate manner, see figure 1. Importantly, the structure and dynamics of the input network are either unknown (black-box) or only partially known (grey-box). Stochastic control can be sought over probability distributions of the desired biochemical species (weak control), or at the level of the underlying sample paths (strong control). See also section S2 in the electronic supplementary material, where we rigorously formulate concepts for stochastic biochemical control.
Molecular controllers must satisfy a set of constraints arising from biochemistry in order to be experimentally realizable. At the structural level, controllers should consist of up to second-order (bi-molecular) reactions [28]; at the kinetic level, the underlying reaction rates should be tunable over a sufficiently large range. At the dynamical level, assuming stability, a controller is said to be robust if the dynamics of desired species in the output network does not explicitly depend on the initial conditions, nor on the rate coefficients from the input network; if these two conditions are satisfied only up to an arbitrarily small error under a suitable timescale separation, then the controller is said to be quasi-robust, see section S2 in the electronic supplementary material for more details. (Quasi-) robustness ensures that the output network accurately traces the predefined dynamics even when the initial conditions are unknown or experimentally difficult to manipulate. For example, the division of a living cell induces extrinsic noise into the initial conditions of the daughter cells [22]; analogously, designing synthetic cells can involve dividing biochemical systems from test-tubes into a larger number of cell-like vesicles with a limited accuracy [9][10][11]. (Quasi-)robustness also ensures that the predefined dynamics can be achieved by fine-tuning the rate coefficients in the controller without having to know the unknown and uncontrollable rate coefficients from a black-box input network. Dynamics achieved under (quasi-) robust controllers are said to display (quasi-)robust adaptation, and plays important roles in biology, e.g. in cell signalling, glycolysis and chemokinesis [29][30][31][32][33][34]; the same is true for timescale separations (slow-fast dynamics), which underpin biochemical multi-stability, oscillations and bifurcations [22][23][24][25]27,[35][36][37][38].
(Quasi-)robust biochemical controllers developed in the literature have been predominantly focused on manipulating the stationary mean (first-moment) of the target species [39][40][41][42] and reducing their stationary variance (secondmoment) [43,44]. This approach is a step forward from controlling the deterministic dynamics to which the underlying stochastic dynamics converges in the thermodynamic limit [45]. However, many important biochemical phenomena, such as cellular differentiation and memory, quorum sensing, bacterial chemokinesis and antibiotic resistance, are realized via higher-order moments of the underlying probability distributions [24,25,35,37,[46][47][48]. For example, as a consequence of containing a gene that can stochastically and transiently (reversibly) switch between two different states, some cells can produce a key regulatory protein whose abundance follows a probability distribution with two modes (maxima), with each of the modes giving rise to a distinct transient cell phenotype. The genetic bi-modality thus gives rise to a biphenotypic population of cells that can have a higher chance of surviving a changing environment [46]; for example, this phenomenon allows some bacteria to persistently survive antibiotic treatments [48]. In this context, particularly important is the number and configuration of the modes present in the probability distributions of the molecular species, and the timing and pattern of stochastic switching in the underlying sample paths. Such dynamically exotic and biochemically important phenomena cannot be ensured using controllers that target only the mean and variance.
To bridge the gap, we put forward a quasi-robust controller called stochastic morpher, presented algorithmically in section S3 in the electronic supplementary material. Stochastic morpher consists of a suitable faster interfacing sub-network de-signed to override the underlying black-box input network and, together with a suitable slower core sub-network, gradually transform (morph) the probability distribution of the desired species into a predened form. This control architecture is based on the phenomenon called noise-induced mixing [25] that some gene-regulatory networks utilize in vivo [37]. Control may be achieved at two different levels: the lowerresolution stochastic morpher incorporates simple production and degradation of the target species as its faster sub-network that, in combination with noise-induced mixing, allows one to achieve multi-modal probability distributions with desired number and configuration of the modes (weak control), and with controlled average timing and mode-switching pattern in the underlying multi-stable sample paths (strong control). The higher-resolution stochastic morpher involves a more complicated faster interfacing network, and can achieve arbitrary probability distributions.
The rest of the paper is organized as follows. In §2, we apply stochastic morpher on the one-species test network (2.1). In §3, we focus on the lower-resolution control in greater detail, by Figure 1. Schematic representation of biochemical control. An input network, R a ¼ R a (X ), is shown in black. The input species X ¼ {X 1 , X 2 , . . . , X N } are divided into two mutually exclusive sets: the target species X t ¼ {X 1 , X 2 , . . . , X n }, that can be explicitly (directly) targeted by the controller and are shown in yellow, and the residual species . , X N }, that can be only implicitly (indirectly) affected by the controller and are shown in white. The coupling between the input species is unknown (respectively, is only partially known) for a black-box (respectively, grey-box) input network. A known coupling between X 1 and X n+2 is depicted as a white dashed double-arrow. A controller is shown, consisting of the sub-network , called the core, and R g ¼ R g (X t , Y), called the interface, displayed in green and red, respectively. The controller core R b specifies how the controlling species . . , Y M }, introduced by the controller and shown in purple, interact among themselves, and is depicted as the green doublearrows. On the other hand, the controller interface R g specifies how the controlling species are coupled with the target species, which is displayed as the red double-arrows. Embedding the controller royalsocietypublishing.org/journal/rsif J. R. Soc. Interface 18: 20200985 explicitly controlling two target species from the three-species test network (3.1). In §4, we apply stochastic morpher on the gene-expression network (4.1), and demonstrate how implicit control can be achieved. In particular, we explicitly influence the mRNA (target species) in a suitable way, ensuring that the translated protein (residual species) is implicitly controlled. In §5, we put forward a blueprint for an experimental realization of stochastic morpher using DNA strand-displacement nanotechnology in order to achieve a bi-phenotypic synthetic cell. Finally, we conclude by presenting a summary and discussion in §6. The notation and background theory utilized in the paper are introduced as needed and are summarized in sections S1-S2 in the electronic supplementary material. General properties of stochastic morpher, outlined via specific examples in §2-4, are rigorously established in Sections S4-S6 in the electronic supplementary material.

Production-degradation input network
Consider the one-species uni-molecular input network where we adopt the convention of denoting two irreversible reactions (e.g. ! X and X ! ) jointly as a single reversible reaction (e.g. O X). In this paper, biochemical species, and their copy-numbers as a function of time t, are represented with upper-case letters (such as X, and X(t), respectively), while the copy-number values are denoted by the corresponding lower-case letters (such as x). Symbol denotes biochemical species that are not explicitly taken into an account. Furthermore, we assume reaction networks are under massaction kinetics [49], with the positive dimensionless rate coefficients displayed above/below the reaction arrows.
In what follows, we fix the (dimensionless) rate coefficients of the input network R 1 a (X) to α = (α 1 , α 2 ) = (1, 1/15). The stationary probability-mass function (PMF) of (2.1), describing the long-time dynamics of the input network and denoted by p ∞ (x), is given by the Poisson distribution with mean (centred at) α 1 /α 2 [26], denoted by p 1 (x) ¼ P(x; a 1 =a 2 ). For the particular choice of the rate coefficients, the Poisson PMF is centred at x = α 1 /α 2 = 15 and is shown as the black dots in figures 2 and 3, which are interpolated with solid black lines for visual clarity. In the rest of this section, we embed different forms of stochastic morpher into the input network (2.1) in order to desirably influence the dynamics of the species X and showcase the capabilities of the controller. Network (2.1) can be interpreted as a simplified model of genetic transcription, with X representing an mRNA species being transcribed and degraded, see §4 for a more-detailed model.

Uni-modality
Consider the stochastic morpher R b < R P g ¼ R b (Y 1 ) < R P g (X; Y 1 ), given by where X is the target species and Y 1 is the controlling species; see figure 1 for a general set-up. Controller (2.2) consists of two sub-networks: the core network R b (Y 1 ), describing a bi-molecular degradation of Y 1 , and the interfacing network R P g (X; Y 1 ) ¼ R 1 g 0 (X; ) < R 1 g 1 (X; Y 1 ), describing a degradation of X, and a production of X catalysed by Y 1 . To emphasize the catalytic role of Y 1 in R 1 g 1 , we write R 1 g 1 ¼ R 1 g 1 (X; Y 1 ); since R 1 g 0 is not catalysed by Y 1 , we write R 1 g 0 ¼ R 1 g 0 (X; ). Let us note that reaction R 1 g 0 is implicitly assumed to be of the form Y 0 + X → Y 0 , where the abundance of an additional controlling species Y 0 is absorbed into an effective rate coefficient γ 0 /ε; see also remarks at the end of this section. The super-script P from R P g stands for the Poisson distribution, as motivated shortly.
In what follows, we analyse the output network R 1 a < R b < R P g , obtained by embedding the stochastic morpher (2.2) into the input network (2.1), which we denote by (2.1) < (2.2). Assuming the copy-number of Y 1 is initially non-zero, the purpose of network R b (Y 1 ) is to fire until the unit copy-number y 1 = 1 is reached; the stationary marginal PMF of the target species X from (2.1)-(2.2) is then given by p 1 (x) ¼ P(x; (g 1 þ 1a 1 )=(g 0 þ 1a 2 )), implying that , as 1 ! 1, Therefore, as the interfacing network R P g fires faster, the input network R 1 a is over-ridden, and the stationary x-marginal PMF of the output network R 1 a < R b < R P g is gradually transformed (morphed) from the Poisson PMF centred at x = α 1 /α 2 to the Poisson PMF centred at x = γ 1 /γ 0 . Note that this uni-modal morphing controls the first-moment (mean) of the output network. In figure 2a, we display the long-time x-marginal PMFs for different values of ε, with the coefficients from R b (Y 1 ) and R P g (X; Y 1 ) fixed to β 1,1 = 1 and γ = (γ 0 , γ 1 ) = (1, 30), respectively. In particular, the long-time PMF is a Poisson distribution centred at x = 24 when ε = 10, shown as the purple squares which, in accordance with (2.3), converges close to the Poisson distribution centred at x = γ 1 /γ 0 = 30 when ε = 10 −2 , shown as the cyan histogram in figure 2a. A representative sample path corresponding to the histogram is displayed in figure 2b, together with the mean which is shown as a red dashed line. The sample path is displayed over a shorter time interval, allowing for the timescale of the underlying fluctuations around the mean to be more discernible.

Bi-modality
Consider now the stochastic morpher R b (Y 1 , Y 2 ) < R P g (X; Y 1 , Y 2 ), given by Sub-network R b (Y 1 , Y 2 ) describes first-order conversion between the two controlling species Y 1 and Y 2 , with the reaction 2Y 1 → Y 1 ensuring that the species Y 1 and Y 2 satisfy the conservation law (y 1 + y 2 ) = 1 in the long run. On the other hand, sub-network R P g (X; Y 1 , Y 2 ) involves two production royalsocietypublishing.org/journal/rsif J. R. Soc. Interface 18: 20200985 royalsocietypublishing.org/journal/rsif J. R. Soc. Interface 18: 20200985 reactions for the species X, one catalysed by Y 1 and the other by Y 2 . Ignoring the reaction 2Y 1 → Y 1 , note that (2.4) can be interpreted as describing a gene switching between two different states Y 1 and Y 2 , and producing an mRNA species X at different rates [25,37]. When ( y 1 , y 2 ) = (1, 0), reaction R 1 g 2 from (2.4) cannot fire, and the remaining faster reactions generate the Poisson PMF centred at x = γ 1 /γ 0 ; analogously, when (y 1 , y 2 ) = (0, 1), the Poisson PMF centred at x = γ 2 /γ 0 is induced. As the controlling species Y 1 and Y 2 convert between themselves, they mix the two Poisson PMFs and over-ride the input network (2.1); the resulting long-time x-marginal PMF of the output see theorem S5.1 in section S5 in the electronic supplementary material for a general result. Let us stress that bi-modality in (2.5) is achieved by exploiting the stochasticity and discreteness of the controlling species Y 1 and Y 2 ; see [25] for more details on this noise-induced mixing phenomenon. In particular, the output networks (2.1) < (2.4) displays a drastically different behaviour at the deterministic level (reaction-rate equations [49]): as a consequence of the deterministic dynamics and continuous abundances, the target species X is uni-stable with the unique equilibrium approaching zero as ε → 0. PMF (2.5) is a linear combination of two Poisson distributions (bases), whose modes are controlled with the rate coefficients γ from the interfacing network R P g (X; Y 1 , Y 2 ), while the PMF values at the modes (weights in (2.5)) are controlled with the rate coefficients β from the controller core R b (Y 1 , Y 2 ). In particular, PMF (2.5) is independent of the asymptotic parameter ε, and depends on β 1,2 and β 2,1 only via the ratio β 1,2 /β 2,1 , which determines the PMF values at the two modes, which in turn depend on the ratios γ 1 /γ 0 and γ 2 /γ 0 . However, the underlying sample paths do depend on ε, which determines the timescale of the fluctuations near each of the two modes, and on the parameters β 1,2 and β 2,1 , which influence the sample paths independently and not only via their ratio. In particular, for a fixed ratio β 1,2 /β 2,1 , the value of β 1,2 determines the timescale of stochastic switching between the two modes: given the system has started near γ 1 /γ 0 , it deterministically moves to a neighbourhood of γ 2 /γ 0 after a random holding time which is exponentially distributed with mean 1/β 1,2 . These observations are instances of the fact that PMFs do not uniquely capture time parametrizations of the underlying sample paths, which we exploit for gaining strong control. Given a fixed ratio β 1,2 /β 2,1 , the precise values of the coefficients can be fixed with the constraint β 1,2 + β 2,1 = c > 0; unless otherwise stated, we set c = 1.

Tri-modality
Stochastic morpher can be utilized to achieve multi-modality beyond bi-modality at the PMF level, and a controlled switching pattern at the underlying sample path level. For example, let us morph the PMF of the input network (2.1) into a trimodal one, with the modes x ∈ {5, 30, 55}. Furthermore, let the underlying sample paths spend on average 3 time-units in the neighbourhood of each of the modes, with the switching order 5 → 55 → 30, i.e. after being close to the mode x = 5, the system should jump near x = 55, then close to x = 30, before finally returning back to x = 5. To this end, consider Analogous to figure 2a-f, in figure 2g,h we display the longtime x-marginal PMF, and a representative sample path, of the output networks (2.1) < (2.6), with γ = (γ 0 , γ 1 , γ 2 , γ 3 ) = (1, 5, 55, 30) and (β 1,1 , β 1,2 , β 2,3 , β 3,1 ) = (1, 1/3, 1/3, 1/3). In the asymptotic limit ε → 0, the PMF approaches the linear combination of three Poisson distributions centred at x = γ 1 /γ 0 = 5, x = γ 2 /γ 0 = 55 and x = γ 3 /γ 0 = 30, each with equal weights (see theorem S5.1 in section S5 in the electronic supplementary material), which is in excellent agreement with the histogram from figure 2g, where parameter ε is two orders of magnitude larger than the rate coefficients from R 1 a (X) and R b (Y 1 , Y 2 , Y 3 ). Note that the switching order of the sample path from figure 2h mirrors the periodic conversion ). Let us stress that, despite having drastically different shapes and functionalities, the uni-modal and tri-modal PMFs displayed in figure 2a,g, respectively, have identical limiting mean (first moment) given by 30, which we show as red dashed lines in figure 2b,h. From the perspective of controllers that can manipulate only the first moment, control achieved in figure 2a,g is indistinguishable. On the other hand, using the stochastic morpher, a significantly finer control can be exerted by manipulating the higher-order moments and, thus, achieving greater biochemical functionality.

Higher-resolution control
In §2.1, we have applied the lower-resolution control consisting of the networks R b and R P g , given generally by (S13) and (S14) in section S3 in the electronic supplementary material, respectively. In this section, we replace the uni-molecular lower-resolution (Poisson distribution) interfacing network R P g with its bi-molecular higher-resolution (Kronecker-delta distribution) counterpart R d g , given by (S15) in section S3 in the electronic supplementary material.

Kronecker-delta distribution
Consider the higher-resolution stochastic morpher consists of two subnetworks: R m,1,s g 0 describes a production of X, a reversible conversion of X into a controlling species Z 1 , and a reversible conversion of X and Z 1 into another controlling species Z 2 . On the other hand, R m,1,s g 1 describes an irreversible conversion of Z 2 into Z 1 , catalysed by Y 1 . Under suitable conditions on the rate coefficients (in particular, with μ 2 γ 0,1 γ 0,2 γ 1 = (σε) −1 and 0 < μ ≪ ε, σ ≪ 1), ensuring that the controlling species Z ¼ {Z 1 , Z 2 } are sufficiently fast, network R d g from (2.7) reduces to

Hybrid control
The lower-and higher-resolution networks R P g and R d g from section S3 in the electronic supplementary material, respectively, can be combined into a composite hybrid scheme for biochemical control. For example, one may wish to obtain a more-detailed control over regions of the statespace where the target species are in lower copy numbers, while a less-detailed control may be sought over the higher copy-numbers state space. Such a hybrid approach may be experimentally desirable, as biochemical realizations of the Kronecker-delta PMFs centred at lower copy-numbers are less expensive to engineer, since a smaller number of the controlling species Z is required. For example, in figure 3c, we display morphing of the input PMF into a mixture of the Kronecker-delta distribution at x = 1 and the Poisson distribution at x = 30 via the controller (S57) from section S5 in the electronic supplementary material, while figure 3d shows a corresponding sample path.

Remarks about the stochastic morpher
The lower-resolution stochastic morpher is a negativefeedback controller; in particular, it introduces degradation reactions, and hence negative self-feedback loops, involving the target species X t . The degradation reactions, such as X ! in (2.2), can be seen as approximations of Y 0 + X → Y 0 , where Y 0 is an auxiliary exclusive catalyst, or Y 0 þ X ! , where Y 0 is a higher-concentration buffer. Alternatively, one can replace X ! with Y 1 + X → Y 1 , without changing the conclusions made in this section, but at a possible experimental cost, see §5. Analogous conclusions hold for the zero-order production reactions in the higher-resolution controller. Note that the higher-resolution controller contains both negative and positive feedback loops involving the species X t and Z. Let us also note that we require the total initial copy number of the controlling species to be non-zero; such a condition is unavoidable, as any network, once implemented, depends on the presence of suitable (buffer) species.

Bi-stable input network
In this section, we apply stochastic morpher to a more complicated reaction network, highlighting how multiple target species can be jointly controlled. To this end, consider the three-species bi-molecular input network R 2 a ¼ R 2 a (X 1 , X 2 , X 3 ), given by (3:1) It is assumed that we can explicitly control the target species X t ¼ {X 1 , X 2 }, with the remaining (residual) species being X r ¼ {X 3 }. For a particular choice of the input rate coefficients α, the long-time (x 1 , x 2 )-marginal PMF of network (3.1) is shown in figure 4a, while the underlying sample paths for X 1 and X 2 are shown in cyan and red in figure 4b, respectively. The (x 1 , x 2 )-marginal PMF is bi-modal, with the approximate modes (x 1 , x 2 ) = (10, 40) and (x 1 , x 2 ) = (40, 10), and the species X 1 and X 2 are negatively correlated.
In figure 4e,f, we set ε = 10 −2 , and one can notice an excellent match with the asymptotic prediction (3.3).
To gain more quantitative information about the convergence, in figure 4g we display the distance (error) between the long-time PMF of the output networks (3.1) < (3.2), denoted by p ε (x 1 , x 2 ), and the target PMF (3.3), as a function of the asymptotic parameter ε. Measuring the error with the l 1 -norm, kp 1 À qk l1 ¼ P x1,x2 jp 1 (x 1 , x 2 ) À q(x 1 , x 2 )j, one can notice that kp 1 À qk l1 ¼ O(1) for sufficiently small ε. In fact, more generally, assuming suitable stability of the output network, convergence of the time-dependent output PMF to the target PMF is royalsocietypublishing.org/journal/rsif J. R. Soc. Interface 18: 20200985 exponential in time t and linear in parameter ε, see theorem S5.1 in section S5 in the electronic supplementary material.
Note that the explicit control of the target species X t ¼ {X 1 , X 2 } implicitly influences the residual species X r ¼ {X 3 }. We demonstrate this influence in figure 4(h), where the long-time x 3 -marginal PMFs from the input network (3.1) is displayed as the solid black curve, while from the output networks (3.1) < (3.2) with ε = 10 −2 as the cyan histogram, showing that the PMF is redistributed across two approximately fixed modes. Also shown, as the dotted red curve, is a prediction of the x 3 -marginal PMF in the limit ε → 0, based on the results from section S6 in the electronic supplementary material; in particular, see theorem S6.1 and example S6.1 in the electronic supplementary material.

Implicit control: gene expression input network
In §3, we focused on explicitly controlling the target species, while ignoring the underlying residual effects. In this section, we shift our focus to an implicit control of the residual species via appropriate explicit manipulation of the target species.
In particular, we show that such an implicit control can be achieved via the lower-resolution stochastic morpher by sufficiently slowing down the controller core R b , in addition to sufficiently speeding up the controller interface R g . To this end, consider the two-species uni-molecular reaction network R 3 a ¼ R 3 a (X 1 , X 2 ), given by Network (4.1) can be seen as an extension of (2.1), describing a simplified model for an intracellular gene expression [37]: X 1 represents an mRNA, transcribed from a gene and translated into a protein X 2 , with each of the two species being degradable. In figure 5a,b, we display the uni-modal stationary PMFs of the target and residual species X 1 and X 2 from (4.1), respectively, for a particular choice of the input coefficients α. The goal in this section is to implicitly induce bi-modality into the protein (residual) species X 2 (and thereby create cells with two phenotypes) by explicitly influencing the mRNA (target) species X 1 . Such a setting is suitable when RNA-based controllers are royalsocietypublishing.org/journal/rsif J. R. Soc. Interface 18: 20200985 utilized which, due to their biophysical properties, are generally more readily interfaced with mRNAs than with proteins [7,17]. As opposed to explicit control, considered in § §2 and 3, where black-box input networks have been considered, implicit control can be applied only to grey-box networks, as one requires at least partial knowledge about the underlying structure and dynamics, see also section S6 in the electronic supplementary material. In what follows, we assume that the reactions which change the copy-number of X 2 are known and given by X 1 → X 1 + X 2 and X 2 ! ; we allow the underlying rate coefficients (α 3 , α 4 ) to remain unknown. To implicitly induce bi-modality into the PMF of X 2 , let us consider the stochastic morpher In the limit ε → 0, the effective behaviour of species X 2 is described by the so-called residual network , obtained by averaging out the equilibrated target species X 1 , and given by R 3 a : X 2 À ! a4 , see theorem S6.1 in section S6 in the electronic supplementary material. By further letting η → 0 in (4.2), one finds that the x 2 -marginal PMF converges close to see theorem S6.2 in section S6 in the electronic supplementary material. In contrast to the explicitly achieved PMFs in § §2 and 3, the implicitly achieved PMF (4.4) is not robust, as its modes depend explicitly on the input rate coefficients α 3 and α 4 ; however, note that the weights in (4.4) are robust. Furthermore, for a fixed ratio (β 1,2 /β 2,1 ), the free parameter (β 1,2 + β 2,1 ) is utilized in § §2 and 3 to achieve strong control of the underlying sample paths; in this section, we utilize this degree of freedom to sufficiently slow down the network Y 1 O Y 2 , i.e. we take η(β 1,2 + β 2,1 ) ≪ 1 to achieve implicit weak control. While the modes from (4.4) are not robust, their ratio, (γ 2,1 / γ 1,1 ), is robust, see corollary S6.1 in section S6 in the electronic supplementary material for a more general result. Such robustness is inherited from the target species; in particular, the ratio between the modes of the implicitly controlled species (protein) is identical to the ratio between the modes of the explicitly controlled species (mRNA). Taking (γ 2,1 /γ 1,1 ) sufficiently large ensures that the two Poisson distributions from (4.4) are wellseparated, and bi-modality achieved. Assume we want the larger mode to be three times further away than the smaller mode, with the two Poisson distributions having equal weights. Such constraints are achievable by the output networks (4.1) < (4.2) with γ = (γ 0,1 , γ 1,1 , γ 2,1 ) = (1, 1, 3) and (β 1,1 , ηβ 1,2 , ηβ 2,1 ) = (1, η/2, η/2), with 0 < η ≪ 1. In figure 5a, the cyan histogram displays the long-time x 1 -marginal PMF of (4.1)<(4.2) when ε = 10 −2 . On the other hand, in figure 5b, we display the long-time x 2 -marginal PMF of the output networks (4.1) < (4.2) with ε = 10 −2 for two different values of η, and one can notice convergence close to the target (4.4). Note that the output x 1 -marginal PMF is uni-modal, as the underlying two Poisson distributions are not well-separated; on the other hand, the limiting x 2 -marginal PMF is bi-modal, as the modes are magnified by a factor of α 3 /α 4 = 10. Let us also note that, if the ratio α 3 /α 4 is experimentally measured for the grey-box input network (4.1), then each of the protein modes, and not only their ratio, can be implicitly controlled; in fact, α 3 /α 4 can be deduced from the action of the stochastic morpher.

Proposed experimental implementation: synthetic cells
In this section, we put forward a blueprint for an in vitro experimental implementation of stochastic morpher via dynamic nucleic acid nanotechnology based on strand- royalsocietypublishing.org/journal/rsif J. R. Soc. Interface 18: 20200985 displacement reactions. Structurally, stochastic morpher involves up to bi-molecular reactions, some of which are catalytic (non-elementary), while its dynamical operation depends on time-separated kinetics. DNA/RNA strand-displacement can in principle meet both of these requirements; in particular, up to bi-molecular reactions, with arbitrary product and reactant species compositions, including catalytic reactions, are readily mapped to elementary reactions when compiled into nucleic-acid-based physical networks [14,[51][52][53]. On the other hand, the rate coefficients can be varied over at least six orders of magnitude under strand displacement [16,17,54], allowing for timescale separations. The lower-resolution control, consisting of the networks (S13)-(S14) in section S3 in the electronic supplementary material, is biochemically less costly, as it involves only two timescales and only one bimolecular reaction; its higher-resolution counterpart, given by (S13)-(S15), is biochemically more costly, since it contains a larger number of timescales, bi-molecular reactions and auxiliary species. Let us stress that quasi-robust control can be achieved with stochastic morpher by tuning ratios and orders of magnitude of the underlying rate coefficients, and not their precise values, as outlined in § §2-4, which is a task achievable within strand-displacement DNA/RNA computing. Furthermore, stochastic morpher transforms PMFs through a sequence of intermediate probability distributions that increasingly resemble the desired target distribution; therefore, satisfactory partial control may be achievable even under weaker timescale separations, e.g. see the purple intermediate PMFs from figure 2, which already display some desirable properties.

DNA Holliday junction and vesicle encapsulation
As a proof-of-concept, we propose an experimental realization of the lower-resolution stochastic morpher (2.4) from §2, which is capable of achieving bi-modality. One way to realize (2.4) is via a biochemical network satisfying the following two properties: (i) it contains two isomeric molecular species which can inter-convert (realizing the controlling species Y 1 and Y 2 , and the reaction Y 1 O Y 2 ), each triggering a catalytic production of the target species X at generally different rates (realizing the reactions Y 1 → Y 1 + X and Y 2 → Y 2 + X), and (ii) the network is integrated into an environment with exactly one copy number of the two isomeric species (realizing the conservation law ( y 1 + y 2 ) = 1, and eliminating the reaction 2Y 1 → Y 1 ). We propose to achieve these conditions using a DNA complex, known as the Holliday junction molecule, encapsulated in a nanoscale compartment, known as the small unilamellar vesicle (SUV) [9][10][11][12]55], forming a bi-phenotypic synthetic cell schematically displayed in figure 6a. The Holliday junction consists of four double-stranded arms crossing at a branch point, designed to be fixed (non-migratory) for our purposes.
In the presence of magnesium ions, the Holliday junction can adopt two distinct orientations, known as stacked conformational isomers [19,56], realizing the reversible reaction On the other hand, the SUV encapsulation is an experimentally demonstrated method for isolating and observing the dynamics of individual molecules with minimal effect from the external environment [12]. Under a sufficiently low concentration of Y 1 and Y 2 during vesicle assembly, one can utilize single-molecule fluorescence [57] to identify SUVs containing exactly one copy number of the Holliday junction.

Strand-displacement reactions
In order to control the rate coefficients of the reactions PEG layer glass g 0 e e Figure 6. Proposed experimental scheme for the stochastic morpher (2.4). (a) Displays an SUV immobilized onto a PEG or BSA passivated surface, containing an input network R a , enclosed in a dashed rounded rectangle, which is coupled to the controller R b,g , included in the red rounded rectangle. The controller includes a single DNA Holliday junction molecule, which switches between two distinct orientations, denoted by Y 1 and Y 2 , and catalytically produces the target species X, which is then degraded by suitable strand-displacement reactions. (b) Displays in a greater detail the underlying strand-displacement reactions that produce X, and which are enclosed in the red rounded rectangle in (a). The single-stranded DNA overhangs on the Holliday junction associate pairwise, forming a toehold and a branch-migration domain, shown in grey/black and blue on Y 1 and Y 2 , respectively, which are separated by a duplex. The association-activated toeholds bind to an auxiliary double-stranded DNA species X, initiating the release of the target molecule X; Y 1 and Y 2 are then recovered via strand-displacement reactions initiated by a suitable recovery strand.
royalsocietypublishing.org/journal/rsif J. R. Soc. Interface 18: 20200985 suitable DNA overhangs involved in the associative toehold activation [58]. More precisely, we propose to extend (tag) three arms of the four-armed Holliday junction with distinct single-stranded DNA molecules (overhangs), which can associate in a pairwise manner, thereby activating distinct toeholds shown as the grey/black regions attached to Y 1 and Y 2 in figure 6b. Furthermore, the overhangs also partially hybridize with each other, forming duplexes next to the toeholds. By controlling the length, and therefore the binding energy, of the duplexes in the associated DNA overhangs, one can experimentally tune the kinetics of reaction On the other hand, by controlling the length of the association-activated toeholds, one can independently tune the rates of Y 1 → Y 1 + X and Y 2 → Y 2 + X, by influencing the two subsequent strand-displacement reactions which produce the target species X from a suitable precursor molecule, denoted by X in figure 6b. Once a molecule of X is produced, Y 1 and Y 2 can be recovered via strand-displacement reactions initiated by a suitable single-stranded DNA, called a recovery strand, ensuring an effective catalytic role of Y 1 and Y 2 in the overall reaction cascade. Let us note that, despite a similarity in the domains of X and the recovery strand, the two strands have distinct toeholds, which we show in light blue and black in figure 6b, respectively; therefore, X and the recovery strand are biochemically distinct. Hence, one can design downstream strand-displacement reactions that would be activated by the unique toehold on X and subsequently monitored. For example, X could invade the desired DNA duplex labelled with a fluorophore and a quencher; the displaced fluorophore-tagged DNA strand can then be detected through the increase in fluorescence.
The operational timescale of the stochastic morpher is limited by the abundance of the DNA complexes that act as substrates for activation and degradation of X, and for recovery of Y 1 and Y 2 . When a significant fraction of the resources is consumed, the underlying effective rate coefficients slow down and the action of the controller weakens. This time limitation can be relaxed by replenishing the resources either externally by using semi-permeable SUVs with porous membranes [59], or internally by encapsulating suitable DNA buffers [60]. In fact, the proposed gene-like implementation of reaction Y 1 O Y 2 , via Holliday junction molecules, does not require buffer species, and in principle has a longer operational time when compared to the DNA strand-displacement implementation of Y 1 O Y 2 put forward in [14].

Discussion
In this paper, we have introduced a molecular controller (see section S2 in the electronic supplementary material and figure 1 for a general outline), called stochastic morpher, that can, under suitable conditions, gradually transform (morph) the PMF of a given black-box input network into a predefined shape. Two forms of the controller have been put forward: the lower-resolution controller, given by (S13)-(S14) in section S3 in the electronic supplementary material, that morphs input PMFs into linear combinations of Poisson distributions, and its higher-resolution counterpart, given by (S13)-(S15) in section S3 in the electronic supplementary material, that achieves linear combinations of Kronecker-delta distributions. The control can be accomplished explicitly on the target species, or implicitly on the residual species via the target species. We have showcased the capabilities of the stochastic morpher on particular biochemical networks in § §2-4. More broadly, we have established general properties of the stochastic morpher in section S4 in the electronic supplementary material using singular perturbation theory [61], and specialized the results to explicit and implicit control in sections S5 and S6 in the electronic supplementary material, respectively. More precisely, we have shown in theorems S5.1-S5.2 that, under the assumption that the output networks display a suitable form of stability, both lower-and higher-resolution stochastic morphers achieve the desired explicit control. In particular, the PMF of the explicitly controlled species converges exponentially fast in time to the desired target distribution, up to an error that decreases linearly in the asymptotic parameter ε. Importantly, the target distributions are independent of the input rate coefficients, and the convergence is robust to the initial conditions, implying long-time quasi-robustness of the stochastic morpher. In theorem S5.3, we have proved that the required stability condition, and hence the results from theorems S5.1-S5.2, hold for any first-order input network. Analogous properties have been established for a class of implicitly controlled residual networks in theorem S6.2. To the best of our knowledge, stochastic morpher is the first controller that can exert control beyond the mean and variance, being able to achieve any probability distribution; in particular, one can achieve desired multi-modal distributions (weak control) whose sample paths have controlled timing and mode-switching pattern (strong control). Furthermore, stochastic morpher can also accomplish implicit control of the residual species, i.e. control of biochemical species that do not directly interact with the controlling species; in contrast, the direct interaction is necessary in some other works, e.g. in the robust controller from [39].
The lower-resolution stochastic morpher incorporates simple production and degradation of the target species as its faster sub-network. Faster production and degradation processes are omnipresent in living cells, where the continuous turnover of RNA molecules and proteins naturally shields the cells against internally and externally induced fluctuations, and reduces the retro-active load placed on the molecular circuits [52]. In control theory context, rapid production and degradation processes are known as high-gain negativefeedback controllers [62]. A high-gain feedback controller has been previously put forward for destroying all but one stable equilibrium, in a given class of multi-stable gene-regulatory networks, thus achieving uni-stability [41,42]. In this paper, we have instead combined high-gain feedback control with noise-induced mixing to systematically morph probability distributions into any predefined shape, with a particular focus on multi-modality/multi-stability. On the other hand, the higher-resolution stochastic morpher shares some similarities with the work presented in [63], where Kronecker-delta distributions are also implemented using timescale separations and catalysis similar to [25,64]. However, while we have realized Kronecker-delta distributions with experimentally feasible bi-molecular networks, implementations from [63] involve higher than bi-molecular networks. In electronic supplementary material, section S4, using the results from [28], we have proved that the bi-molecular construction put forward in this paper reduces in an appropriate limit to the higher-molecular one from [63]. Furthermore, in this paper, we have focused on biochemical control, involving combining multiple biochemical networks together, some of which are black-box (unknown), in order to achieve suitable dynamical behaviours. Therefore, our analyses take into a consideration realistic effects biochemical networks experience in applications, such as retro-activity [52,65]; such effects have not been included in [63].
In §5, we have put forward a blueprint for an in vitro experimental realization of the bi-modal stochastic morpher (2.4) using a DNA-based biochemical network encapsulated inside a nano-scale vesicle, thus implementing a bi-phenotypic synthetic cell shown in figure 6. We have proposed to appropriately assemble the vesicles so that the reaction 2Y 1 → Y 1 from (2.4) can be eliminated. However, let us note that such a reaction can be achieved within the proposed set-up by replacing the Holliday junction molecules with suitable DNA dimers. In particular, we propose assembling sufficiently rigid dimers containing the Holliday junction and a corresponding complementary motif (deactivator) on the opposite end, see [66] for details on multimer assembly of DNA nanostructures. When multiple dimers are present, a sequence of associations between the dimers would take place, sequestering all but one active Holliday junction, thus implementing 2Y 1 → Y 1 . Let us also stress that, in this paper and §5 in particular, we have assumed that the molecular abundances are spatially homogeneous (well-mixed) inside reactors at any given time. This assumption holds if the stochastic morpher is implemented in sufficiently small compartments. However, for larger compartments, such as living cells, spatial heterogeneity of molecular abundances can play important roles [67][68][69][70][71]. To address such challenges, further theoretical and experimental considerations are required, which are beyond the scope of this paper.
Neglecting spatial heterogeneity in the molecular abundances, the lower-resolution stochastic morpher can also in principle be implemented in vivo via synthetic RNA-based networks encoded genetically inside cells; for concreteness, we focus on the bi-modal case (2.4). We propose introducing a synthetic plasmid, coding for a target micro-RNA species, into an Escherichia coli cell. In this context, the slowly interconverting unit copy-number controlling species are realized via different states of a suitable plasmid gene, and the switching between the two states is regulated by the slower binding and unbinding of suitable transcription factors giving rise to different degrees of promoter activity. The production of the target RNA is achieved via transcription from the plasmid gene, with the maximum transcription rate of the order of 10 −1 s −1 [72], and we propose exploiting the natural intracellular degradation of RNA molecules, occurring at a rate of around 3 × 10 −3 s −1 [73]. These estimations imply that the stochastic morpher R b,g can achieve multi-modal distributions with modes of up to approximately 30 RNA molecules per cell. To ensure quasi-robustness of the controller R b,g , we require that the desired input network R a fires at a slower timescale. For the interactions between the RNA species at a concentration of 30 molecules per cell, 0.1 reactions per second corresponds to a rate constant of the order of 10 6 M À1 s À1 , whereas the speed limit of nucleic-acid reactions inside cells is of the order of 10 8 M À1 s À1 [74]. If R a is a synthetic RNA-based network, the quasi-robust limit can be achieved by slowing down the input network via sequestration of the key domains in the secondary structures of the underlying species [17,54]. On the other hand, if R a is a native RNA-based network, the quasi-robust limit might be challenging to achieve. In this context, an approach for speeding up the stochastic morpher is by transcribing sufficiently long RNA chains that are broken down into multiple target RNA molecules post-transcription via suitable ribozymes [75].
Data accessibility. Data supporting this paper are contained within the paper, and as the electronic electronic supplementary material.