Feedback control theory facilitates the development of self-regulating systems with desired performance which are predictable and insensitive to disturbances. Feedback regulatory topologies are found in many natural systems and have been of key importance in the design of reliable synthetic bio-devices operating in complex biological environments. Here, we study control schemes for biomolecular processes with two outputs of interest, expanding previously described concepts based on single-output systems. Regulation of such processes may unlock new design possibilities but can be challenging due to coupling interactions; also potential disturbances applied on one of the outputs may affect both. We therefore propose architectures for robustly manipulating the ratio/product and linear combinations of the outputs as well as each of the outputs independently. To demonstrate their characteristics, we apply these architectures to a simple process of two mutually activated biomolecular species. We also highlight the potential for experimental implementation by exploring synthetic realizations both in vivo and in vitro. This work presents an important step forward in building bio-devices capable of sophisticated functions.

1. Introduction

For more than two decades, we have witnessed significant advances in the highly interdisciplinary field of synthetic biology whose goal is to harness engineering approaches in order to realize genetic networks that produce user-defined cell behaviour. These advances have the potential to transform several aspects of our life by providing efficient solutions to many global challenges related to food security, healthcare, energy and the environment [1–6]. A fundamental characteristic of living systems is the presence of multi-scale feedback mechanisms facilitating their functioning and survival [7,8]. Feedback control enables a self-regulating system to adjust its current and future actions by sensing the state of its outputs, thus maintaining an acceptable response even in the face of unintended and unknown changes. This can be the answer to a number of major challenges [9–11] that prevent the successful implementation of synthetic genetic circuits and keep innovative endeavours in the field trapped at the laboratory stage. Control theory offers a rich toolkit of powerful techniques to design and manipulate biological systems and enable the reliable function of next-generation synthetic biology applications [12–16].

Engineering synthetic gene circuits aims at constructing modular biomolecular devices which are able to operate in a controllable and predictable way in constantly changing environments with a high level of metabolic burden and interactions (cross-talk) with endogenous signalling systems. It is therefore a requirement for them to be resilient to context-dependent effects and adapt to external environmental perturbations. Several control approaches inspired by both natural and technological systems have recently been proposed allowing for effective and robust regulation of biological networks in vivo and/or in vitro [17–23]. Despite some conceptual differences, all of these studies focus on biomolecular systems with one output of interest, such as the expression of a single protein.

Building advanced bio-devices capable of performing more sophisticated computations and tasks requires the design of genetic circuits where multiple inputs are applied and multiple outputs are measured. In control engineering, these types of systems are also known as multi-input multi-output or MIMO systems [24]. This may be the key for achieving control of the whole cell, which can be regarded as a very complex MIMO bio-device itself. Regulation of processes comprising multiple interacting variables of interest can be challenging since there may be interactions between inputs and outputs. Thus, a change in any input may affect all outputs. At the same time, when attempting to apply feedback control by ‘closing the loop’, a quandary arises as to which input should be connected with which output (input–output pairing problem). Addressing such problems therefore requires alternative, suitably adjusted regulation schemes which take into account the presence of mutual internal interactions in the network to be controlled (open-loop system).

The research area of MIMO control bio-systems has up until now remained relatively unexplored. There have been only a few studies towards this direction, associated with cybergenetic approaches where a computer is a necessary part of the control feedback loop [25,26]. By contrast, substantial progress has been made in a closely related area, namely MIMO logic bio-circuits which are able to realize Boolean functions [27,28] while ‘multi-layer/level’ control concepts for one-output processes [29,30] and resource allocation in gene expression [31] have also been proposed.

In this paper, we investigate regulation strategies for biomolecular networks with two outputs of interest which can correspond, for example, to the concentration of two different proteins inside the cell, assuming the presence of mutual interactions. Both the open-loop and the closed-loop system (open-loop system within a feedback control configuration) are represented by chemical reaction networks (CRNs) obeying the law of mass action [8]. Consequently, the entire regulation process takes place in the biological context of interest without the use of computer-aided methods. Our designs take advantage of the antithetic integral motif which was first introduced in [32] and whose properties and performance trade-offs have been extensively studied in various single-output biomolecular systems [33–52]. The antithetic integral motif is able to achieve robust steady-state tracking, which is equivalent to the biological principle of robust perfect adaptation (RPA) [17,53,54], via integral feedback control. A core element of this motif is an (ideally) irreversible sequestration reaction between two species representing a comparison operation at the molecular level. The memory function, necessary for any integral controller, is performed by a memory variable accumulating, through (mathematical) integration, the error between an output and a set-point of interest over time. In the general case, this memory variable is ‘hidden’ in the sense that it corresponds to a non-physical quantity defined as a (mathematical) combination of the (physical) controller species. The efficacy of this biomolecular mechanism has also been demonstrated experimentally in living cells, at both the cell population and the single-cell level, and in cell-free environments using either external (in silico) or embedded single-output control schemes [33,39,40,43,55–61]. Furthermore, in recent years, considerable attention has been given to topologies combining the antithetic integral controller with proportional and derivative control action or biomolecular buffering [39,40,43,60,62–66]. Such efforts seek to resolve commonly encountered issues associated with the standalone antithetic integral controller, such as instability, poor transient dynamics including overshoots, and long-lasting oscillations or increased variance.

One of the main objectives of this work is to show how this molecular sequestration mechanism can be used to regulate biomolecular processes with more than one output, expanding existing theoretical single-output approaches. Thus, we introduce novel strategies of biomolecular interconnections which are able to efficiently control multi-output biological systems in several ways and discuss important challenges and phenomena arising in such contexts. Focusing primarily on two-output biological systems, we present regulatory designs exploiting ‘multi-loop’ concepts based on two independent feedback loops as well as concepts where the control action is carried out jointly considering both outputs simultaneously. Our designs are scalable and, with appropriate modifications, can handle biological systems with an arbitrary number of outputs.

Specifically, we present regulatory architectures, which we refer to as regulators, capable of achieving one of the following control objectives: robustly driving (a) the ratio/product of the outputs, (b) a linear combination of the outputs and (c) each of the outputs to a desired value (set-point). At steady state, the architectures of (a) and (b) result in two coupled outputs which can still affect each other, albeit in a specific way dictated by the respective control approach. On the other hand, the architectures for (c) achieve steady-state decoupling, thus making the two outputs independent of each other. Our control schemes can be used for regulation of any arbitrary open-loop process provided that the resulting closed-loop system has a finite, positive steady state and the closed-loop system converges to that steady state as time goes to infinity (closed-loop (asymptotic) stability). Thus, the present analysis focuses exclusively on such scenarios. Furthermore, we mathematically and computationally demonstrate their special characteristics by applying these schemes to a simple, monomolecular, biological process of two mutually activating species. Finally, to highlight their biological relevance and motivate further experimental investigation, we explore potential implementations of our designs.

2. Control schemes with steady-state coupling

In figure 1a, we show a general biomolecular process with two outputs of interest for which we first present two bio-controllers aiming to regulate the ratio and an arbitrary linear combination of the outputs, respectively. The different types of biomolecular reactions as well as their graphical representations used in this work are presented in figure 1b.

2.1. Regulating the ratio of outputs

Figure 1c illustrates a motif which we call Ratio-Regulator (R-Regulator) and consists of the following reactions:

\begin{matrix} Y_{1} \overset{k_{1}}{⟶} Y_{1} + Z_{1}, Y_{2} \overset{k_{2}}{⟶} Y_{2} + Z_{2}, \\ Y_{2} + Z_{2} \overset{k_{3}}{⟶} Z_{2}, Z_{1} + Z_{2} \overset{η}{⟶} Ø . \end{matrix}}

2.1

This controller consists of two species, Z₁ and Z₂, which annihilate each other. The production of Z₁, Z₂ is catalysed by the target species Y₁, Y₂, respectively, while Y₂ is also inhibited by Z₂.

The dynamics of the R-Regulator are described by the following system of ordinary differential equations (ODEs):

{\dot{Z}}_{1} = k_{1} Y_{1} - η Z_{1} Z_{2}

2.2a

and

{\dot{Z}}_{2} = k_{2} Y_{2} - η Z_{1} Z_{2} .

2.2b

Equations (2.2a)–(2.2b) give rise to a non-physical memory variable which enables integration, i.e.

{\dot{Z}}_{1} - {\dot{Z}}_{2} = k_{1} Y_{1} - k_{2} Y_{2}

(Z_{1} - Z_{2}) (t) = k_{1} \int_{0}^{t} (Y_{1} (τ) - \frac{k_{2}}{k_{1}} Y_{2} (τ)) d τ .

2.3

As a result, assuming closed-loop stability (

{\dot{Z}}_{1}, {\dot{Z}}_{2}, {\dot{Y}}_{1}

{\dot{Y}}_{2} \to 0

as t →∞), we get

\frac{Y_{1}^{*}}{Y_{2}^{*}} = \frac{k_{2}}{k_{1}},

2.4

where the * notation indicates the steady-state concentration of a species. As can be seen, the integrand in equation (2.3) corresponds to an error quantity which converges to zero over time, thus guaranteeing that the output ratio (

Y_{1}^{*} / Y_{2}^{*}

) will converge to the set-point (k₂/k₁). It is important to note that the aforementioned stability depends on the structure of the open-loop process, which is unknown here, as well as the set of the reaction rates/parameter values we select for the closed-loop system.

As revealed by equation (2.4), the R-Regulator is characterized by a dynamic set-point tracking property regarding species Y₁ and Y₂. This property becomes more apparent if we examine the resulting closed-loop architecture from a different viewpoint. Imagine, for instance, that Y₁ represents an input species through which a ‘reference signal’ is applied while Y₂ represents an output (target) species. Then, $Y_{2}^{*}$ is able to track the changes of the set-point ( $k_{1} Y_{1}^{*} / k_{2}$ ) (and vice versa).

A modified version of the above control scheme can be obtained by replacing $Y_{1} \overset{k_{1}}{⟶} Y_{1} + Z_{1}, Y_{2} \overset{k_{2}}{⟶} Y_{2} + Z_{2}$ with $Ø \overset{k_{1}}{⟶} Z_{1}, Y_{1} + Y_{2} \overset{k_{2}}{⟶} Z_{2}$ in CRN (2.1). As a result, the memory variable becomes ${\dot{Z}}_{1} - {\dot{Z}}_{2} = k_{1} - k_{2} Y_{1} Y_{2}$ leading to $Y_{1}^{*} Y_{2}^{*} = k_{1} / k_{2}$ . This modified R-Regulator is able to regulate the product of two outputs, assuming both outputs represent species concentrations. Equivalently, this can be seen as regulation of the ratio of two outputs where one of them represents a species concentration and the other one the reciprocal of a species concentration (see also §S12 of the electronic supplementary material for further demonstration).

2.2. Regulating a linear combination of the outputs

In figure 1d, a second motif, which we call Linear Combination-Regulator (LC-Regulator), is depicted. The only difference from the R-Regulator is that species Z₁, Z₂ are also produced through two independent processes with constant rates θ₁, θ₂, respectively. More specifically, the corresponding reaction network is

\begin{matrix} Ø \overset{θ_{1}}{⟶} Z_{1}, Ø \overset{θ_{2}}{⟶} Z_{2}, Y_{1} \overset{k_{1}}{⟶} Y_{1} + Z_{1}, Y_{2} \overset{k_{2}}{⟶} Y_{2} + Z_{2}, \\ Y_{2} + Z_{2} \overset{k_{3}}{⟶} Z_{2}, Z_{1} + Z_{2} \overset{η}{⟶} Ø . \end{matrix}}

2.5

The dynamics of LC-Regulator is given by the set of ODEs

{\dot{Z}}_{1} = θ_{1} + k_{1} Y_{1} - η Z_{1} Z_{2}

2.6a

and

{\dot{Z}}_{2} = θ_{2} + k_{2} Y_{2} - η Z_{1} Z_{2} .

2.6b

Similar to before, in order to see the memory function involved, we subtract equations (2.6a)–(2.6b) and integrate to get

(Z_{1} - Z_{2}) (t) = \int_{0}^{t} ((k_{1} Y_{1} (τ) - k_{2} Y_{2} (τ)) - (θ_{2} - θ_{1})) d τ .

Under the assumption of closed-loop stability (

{\dot{Z}}_{1}

{\dot{Z}}_{2}

{\dot{Y}}_{1}

{\dot{Y}}_{2} \to 0

as t → ∞), we have at steady state:

k_{1} Y_{1}^{*} - k_{2} Y_{2}^{*} = θ_{2} - θ_{1} .

2.7

An interesting feature of LC-Regulator is that equation (2.7) can be adjusted as desired by modifying the production reactions regarding Z₁, Z₂. A more general formulation of this control scheme providing a full characterization of the possible (steady-state) output combinations is discussed in §S2 of the of the electronic supplementary material.

Finally, in figure 1e, we show an alternative version of the controllers presented above. Specifically, the inhibitory reaction $Y_{1} + Z_{1} \overset{k_{4}}{⟶} Z_{1}$ has been added to the R- or LC-Regulator. Note that this additional reaction does not change the dynamics of the controllers—equations (2.2a) and (2.2b) and (2.6a) and (2.6b) still hold for R-Regulator and LC-Regulator, respectively. Despite the increase in complexity, the additional reaction strengthens the regulatory ability of the controllers in the sense that control action is now applied on both target species. This could, for example, be useful to make closed-loop stability more robust. These slightly modified motifs are further discussed from a stability viewpoint in §6 and §S9 of the electronic supplementary material.

3. Control schemes with steady-state decoupling

We now present three alternative bio-controllers, which we call Decoupling-Regulator (D-Regulator) I, II and III, capable of achieving independent control of each output in the arbitrary biomolecular process (figure 1a). In particular, D-Regulators are able to drive each output species to a desired steady-state concentration unaffected by the behaviour of the other species.

D-Regulators I, II follow a decentralized approach exploiting a ‘multi-loop’ control strategy. More analytically, each of them uses two single-input single-output (SISO) integral controllers which can be constructed separately. This might be advantageous in certain applications in the sense that already-existing, successful SISO implementation techniques can be used. However, in the general case, the two SISO controllers cannot be analysed or tuned independently due to the existence of coupling interactions in the network to be controlled. D-Regulators I, II and, by extension, their resulting closed-loop architectures are MIMO systems and should be studied as such in order for a desirable overall behaviour to be achieved—for instance, in terms of closed-loop stability or dynamic performance of both output responses. Furthermore, in a later section, we investigate a ‘pairing problem’ between actuator and sensor species using a simple example based on one of the above regulators. Problems of such nature are very common in multi-loop contexts and can be difficult to address, especially for complex, strongly coupled networks. Moreover, in D-Regulator I, II the individual SISO controllers constitute alternative realizations of the antithetic integral motif [32]. We choose to focus on these specific versions because of their essential structural differences which can play a crucial role for a circuit designer implementation-wise. At the same time, other well-studied realizations of the antithetic integral motif in the literature appear to be more complex and use the aforementioned versions as a structural basis. A characteristic example is the (SISO) rein controller presented in [36], which is implemented as part of a D-Regulator discussed in §S3 of the electronic supplementary material.

On the other hand, D-Regulator III follows a centralized approach where some parts of the architecture jointly contribute to the realization of integral control on both output species. This control strategy can result in a structurally simpler topology with fewer controller species. Nevertheless, building such a topology might require more sophisticated biomolecular components.

3.1. D-Regulator I

The set of reactions describing D-Regulator I (figure 2a) is

\begin{matrix} Y_{1} \overset{k_{1}}{⟶} Y_{1} + Z_{1}, Y_{2} \overset{k_{2}}{⟶} Y_{2} + Z_{2}, Y_{1} + Z_{1} \overset{k_{3}}{⟶} Z_{1}, \\ Y_{2} + Z_{2} \overset{k_{4}}{⟶} Z_{2}, Ø \overset{θ_{1}}{⟶} Z_{3}, Ø \overset{θ_{2}}{⟶} Z_{4}, \\ Z_{1} + Z_{3} \overset{η_{1}}{⟶} Ø, Z_{2} + Z_{4} \overset{η_{2}}{⟶} Ø . \end{matrix}}

3.1

This design comprises four controller species. The target species Y₁, Y₂ catalyse the formation of two of them, Z₁, Z₂, which, in turn, inhibit the former. In addition, Z₃, Z₄, which are produced independently at a constant rate, participate in annihilation reactions with Z₁ and Z₂, respectively.

Figure 2. Control architectures with steady-state decoupling. Schematic of a general closed-loop architecture using (a) D-Regulator I (CRN (3.1)), (b) D-Regulator II (CRN (3.6)) and (c) D-Regulator III (CRN (3.7)).

The dynamics of D-Regulator I can be modelled using the following set of ODEs:

{\dot{Z}}_{1} = k_{1} Y_{1} - η_{1} Z_{1} Z_{3},

3.2a

{\dot{Z}}_{2} = k_{2} Y_{2} - η_{2} Z_{2} Z_{4},

3.2b

{\dot{Z}}_{3} = θ_{1} - η_{1} Z_{1} Z_{3}

3.2c

and {\dot{Z}}_{4} = θ_{2} - η_{2} Z_{2} Z_{4} .

3.2d

By contrast to the regulation strategies presented in the preceding section, D-Regulator I includes two memory variables which carry out integral action independently. Indeed, combining equations (3.2a) and (3.2c) results in

(Z_{3} - Z_{1}) (t) = k_{1} \int_{0}^{t} (\frac{θ_{1}}{k_{1}} - Y_{1} (τ)) d τ

3.3

while combining equations (3.2b) and (3.2d) gives

(Z_{4} - Z_{2}) (t) = k_{2} \int_{0}^{t} (\frac{θ_{2}}{k_{2}} - Y_{2} (τ)) d τ .

3.4

Consequently, the steady-state output concentrations under the assumption of closed-loop stability ( ${\dot{Z}}_{1}$ , ${\dot{Z}}_{2}$ , ${\dot{Z}}_{3}$ , ${\dot{Z}}_{4}$ , ${\dot{Y}}_{1}$ , ${\dot{Y}}_{2}$ → 0 as t → ∞) are

Y_{1}^{*} = \frac{θ_{1}}{k_{1}} and Y_{2}^{*} = \frac{θ_{2}}{k_{2}} .

3.5

3.2. D-Regulator II

By using four controller species as before, we construct D-Regulator II (figure 2b) consisting of the following reactions:

\begin{matrix} Y_{1} \overset{k_{1}}{⟶} Y_{1} + Z_{1}, Y_{2} \overset{k_{2}}{⟶} Y_{2} + Z_{2}, Ø \overset{θ_{1}}{⟶} Z_{3}, \\ Ø \overset{θ_{2}}{⟶} Z_{4}, Z_{3} \overset{k_{3}}{⟶} Z_{3} + Y_{1}, Z_{4} \overset{k_{4}}{⟶} Z_{4} + Y_{2}, \\ Z_{1} + Z_{3} \overset{η_{1}}{⟶} Ø, Z_{2} + Z_{4} \overset{η_{2}}{⟶} Ø . \end{matrix}}

3.6

In this case, species Z₃, Z₄ catalyse the formation of the target species Y₁, Y₂, respectively, and Z₃, Z₄ are produced at a constant rate. Furthermore, species Z₁, Z₂ are catalytically produced by Y₁, Y₂, respectively, while the pairs Z₁–Z₃ and Z₂–Z₄ participate in an annihilation reaction.

Note that the species of D-Regulator II are described by the same ODE model as D-Regulator I (equations (3.2a)–(3.2d)). Thus, the memory variables involved (equations (3.3) and (3.4)) as well as the steady-state output behaviour (equation (3.5)) are identical in these two motifs (provided that closed-loop stability is guaranteed). Nonetheless, in general, regulating the same open-loop process via the aforementioned controllers results in different output behaviour until an equilibrium is reached or, in other words, the transient responses differ. This is because of the different topological characteristics of the two motifs which cannot be captured by focusing only on the controller dynamics: considering closed-loop dynamics is required, which is addressed in a later section.

3.3. D-Regulator III

The last bio-controller presented in this study is D-Regulator III (figure 2c) whose structure is composed of the following reactions:

\begin{matrix} Y_{1} \overset{k_{1}}{⟶} Y_{1} + Z_{1}, Y_{2} \overset{k_{2}}{⟶} Y_{2} + Z_{2}, Ø \overset{θ_{1}}{⟶} Z_{3}, Z_{3} \overset{k_{3}}{⟶} Z_{3} + Y_{1}, \\ Y_{2} + Z_{2} \overset{k_{4}}{⟶} Z_{2}, Z_{1} + Z_{3} \overset{η_{1}}{⟶} C, Z_{2} + C \overset{η_{2}}{⟶} Ø . \end{matrix}}

3.7

Here, there are three controller species. Z₁, Z₃ interact with the target species Y₁ as well as with each other in the same way as in D-Regulator II. The complex C, which is formed by the binding of Z₁, Z₃, and the third controller species, Z₂, can annihilate each other. Finally, the target species Y₂ catalyses the production of Z₂ which, in turn, inhibits Y₂ analogous to D-Regulator I.

The dynamics of D-Regulator III can be described by the following set of ODEs:

{\dot{Z}}_{1} = k_{1} Y_{1} - η_{1} Z_{1} Z_{3},

3.8a

{\dot{Z}}_{2} = k_{2} Y_{2} - η_{2} Z_{2} C,

3.8b

{\dot{Z}}_{3} = θ_{1} - η_{1} Z_{1} Z_{3}

3.8c

and \dot{C} = η_{1} Z_{1} Z_{3} - η_{2} Z_{2} C .

3.8d

Similar to the other D-Regulators, the memory function responsible for the regulation of the output Y₁ is carried out by the (non-physical) quantity Z₃ − Z₁ (equation (3.3)). However, the memory variable related to the output Y₂ is realized in a different way from before. More specifically, combining equations (3.8b)–(3.8d) yields

{\dot{Z}}_{3} + \dot{C} - {\dot{Z}}_{2} = θ_{1} - k_{2} Y_{2}

(Z_{3} + C - Z_{2}) (t) = k_{2} \int_{0}^{t} (\frac{θ_{1}}{k_{2}} - Y_{2} (τ)) d τ .

Therefore, assuming closed-loop stability, i.e. ${\dot{Z}}_{1}$ , ${\dot{Z}}_{2}$ , ${\dot{Z}}_{3}$ , $\dot{C}$ , ${\dot{Y}}_{1}$ , ${\dot{Y}}_{2}$ → 0 as t → ∞, the steady-state output behaviour is

Y_{1}^{*} = \frac{θ_{1}}{k_{1}} and Y_{2}^{*} = \frac{θ_{1}}{k_{2}} .

3.9

4. Specifying the biological network to be controlled

We now turn our focus to a specific two-output open-loop network which will henceforward take the place of the abstract ‘cloud’ process in the preceding sections. This will allow us to implement in silico the proposed control motifs and demonstrate the properties discussed above (see §5). In addition, we will explore potential experimental realizations of the resulting closed-loop networks (see §7).

Figure 3a illustrates a simple biological network comprised of two general birth–death processes involving two target species, Y₁, Y₂. These species are coupled in the sense that each of them is able to catalyse the formation of the other. Such motifs of positive feedback action are ubiquitous in biological systems [68–70]. In particular, we have the reactions

\begin{matrix} Ø \overset{b_{1}}{⟶} Y_{1}, Ø \overset{b_{2}}{⟶} Y_{2}, Y_{1} \overset{d_{1}}{⟶} Ø, Y_{2} \overset{d_{2}}{⟶} Ø, \\ Y_{1} \overset{α_{2}}{⟶} Y_{1} + Y_{2}, Y_{2} \overset{α_{1}}{⟶} Y_{1} + Y_{2} \end{matrix}}

4.1

which can be modelled as

{\dot{Y}}_{1} = b_{1} - d_{1} Y_{1} + α_{1} Y_{2}

4.2a

and

{\dot{Y}}_{2} = b_{2} - d_{2} Y_{2} + α_{2} Y_{1} .

4.2b

Figure 3. Specifying the open-loop biomolecular network. (a) A simple biological process with two mutually activating output species Y₁, Y₂, described by CRN (4.1). (b) Simulated response of the topology in (a) using the ODE model (4.2) with the following parameters: b₁ = 2 nM min⁻¹, b₂ = 1 nM min⁻¹, d₁ = d₂ = 1 min⁻¹, α₁ = 0.1 min⁻¹, α₂ = 0.4 min⁻¹. At time t = 50 min, a disturbance on Y₁ is introduced which affects both output species. More specifically, the value of parameter b₁ changes from 2 to 4.

For any d₁d₂ > α₁α₂, ODE system (4.2a) and (4.2b) has the following unique positive steady state:

Y_{1}^{*} = \frac{α_{1} b_{2} + b_{1} d_{2}}{d_{1} d_{2} - α_{1} α_{2}} and Y_{2}^{*} = \frac{α_{2} b_{1} + b_{2} d_{1}}{d_{1} d_{2} - α_{1} α_{2}}

4.3

which is (globally) exponentially stable (see §S4 of the electronic supplementary material).

Note that for this system, a change in any of the reaction rates of network (4.1) due to, for instance, undesired disturbances, will affect the behaviour of both species Y₁ and Y₂ (figure 3b).

5. Implementing the proposed regulation strategies

We now demonstrate the efficiency of the bio-controllers introduced in §§2 and 3 by regulating the open-loop network (4.1) presented in §4 (see also §8 for regulation of a more complex network). A detailed analysis of the steady-state behaviour of the resulting closed-loop processes can be found in §S5 of the electronic supplementary material.

We show in figure 4 that R-Regulator and LC-Regulator are capable of driving the ratio and a desired linear combination of the output species to the set-point of our choice in the presence of constant disturbances, respectively. Similarly, we illustrate in figure 5 the ability of D-Regulators to robustly steer each of the output species towards a desired value independently, thus cancelling the steady-state coupling.

In the topology shown in figure 5b, there are two actuation reactions realized though Z₃ and Z₄. Due to the existence of coupling interactions in the network that we aim to control, it is evident that these actuator species act on both Y₁ and Y₂ simultaneously. Consequently, one could argue that an alternative way of closing the loop would be through a different species pairing (figure 6). In particular, an annihilation (comparison) reaction between Z₁, Z₄ and Z₂, Z₃ could be used instead (Z₁, Z₂ can be considered as sensor species measuring the outputs Y₁, Y₂, respectively). However, it can be demonstrated (see §S6 of the electronic supplementary material) that this control strategy is not feasible since there is no realistic parameter set that can ensure closed-loop stability.

Figure 6. A different feedback configuration regarding the topology shown in figure 5b. It is based on D-Regulator II with a different actuator–sensor species pairing which leads to instability.

Finally, in the electronic supplementary material, using three closed-loop architectures (one for each regulator type – R, LC and D), we demonstrate through simulations the robust steady-state tracking property of the systems by perturbing several model parameters (see §S7 of the electronic supplementary material). At the same time, we computationally investigate the effect of controller species degradation on their performance and how the latter can be mitigated via appropriate parameter tuning (see §S8 of the electronic supplementary material).

6. Closed-loop stability

As already emphasized, assuming the existence of a finite, positive equilibrium, the proposed regulation strategies require asymptotic closed-loop stability, at least around that equilibrium (locally). A commonly used approach to assess local stability of a nonlinear system is through (Jacobian) linearization. Specifically, we can study the resulting Jacobian matrix [8]. If its eigenvalues have strictly negative real parts, i.e. the matrix is Hurwitz, then the aforementioned equilibrium is locally asymptotically stable. Necessary and sufficient conditions for that can be determined via the Routh–Hurwitz criterion (see the sections of the electronic supplementary material associated with §§4 and 5).

Instead of analysing the system as a whole, we can alternatively examine it as an interconnection of two (or more) subsystems [71,72]. It is often possible to assess the overall stability by studying those subsystems separately. This could be beneficial when only an input–output property of the system to be controlled is known. To demonstrate this, we consider the R-Regulator and LC-Regulator with two inhibitory reactions controlling a general ‘cloud’ network in a negative feedback configuration, as shown in figure 1e. Focusing on the behaviour around an equilibrium of interest, we can show in both cases that if $k_{2} k_{4} Z_{1}^{*} = k_{1} k_{3} Z_{2}^{*}$ , then the ‘controller block’ corresponds to a positive real (PR) system. It is also known [71] that the negative feedback interconnection of a PR block and a weakly strictly PR (WSPR) one yields an overall asymptotically stable system. Consequently, for every WSPR ‘cloud block’, asymptotic closed-loop stability can be guaranteed. Further details including definitions of PR and WSPR concepts as well as proofs can be found in §S9 of the electronic supplementary material.

7. Experimental realization

To highlight the feasibility of experimentally realizing the proposed control schemes, this section describes potential in vivo and in vitro implementations of the open-loop and closed-loop circuits introduced earlier. We first focus on implementations using biological parts that have been characterized in Escherichia coli and further discuss a molecular programming approach.

Following the description in §4, the biological network to be controlled can be realized as shown in figure 7. In this implementation, Y₁ and Y₂ are heterologous sigma factors [73], which are fused to fluorescent proteins (GFP and mCherry) to facilitate tracking of the output. While genes encoding fusion proteins are shown for simplicity, bicistronic constructs could also be used and may be preferred in practice to avoid impairment of sigma factor activity by fusion to a fluorescent protein. Through a suitable choice of promoters, Y₁ mediates the expression of Y₂ and vice versa. Low levels of Y₁ and Y₂ are continuously produced from constitutive promoters, such as promoters from the BioBrick collection [74]. In all following figures, the biological parts underlying these interactions are not explicitly shown.

7.1. R-Regulator and LC-Regulator

For the proposed implementation of the R-Regulator (figure 8),Y₂ mediates expression of the hepatitis C virus protease NS3 fused to maltose-binding protein (MBP) (Z₂). Y₁ facilitates expression of a MBP-single-chain antibody (scFv) fusion (Z₁) that specifically binds to and thus inhibits NS3 protease. Inhibition of NS3 protease activity through coexpression with single-chain antibodies in the cytoplasm of E. coli has been demonstrated previously [75]. Adding a recognition sequence to Y₂ will further allow for its degradation by NS3. Importantly, this will require identification of sites in the Y₂ protein that allow for integration of the NS3 recognition sequence without compromising the catalytic activity of Y₂. An additional requirement for the LC-Regulator would be constitutive expression of malE-scFv and malE-scNS3 as indicated in the dashed boxes in figure 8. It is important to note that binding between the biomolecular species realizing the annihilation reaction should ideally be irreversible, which would likely require targeted engineering of a suitable antibody [76] or exploration of alternative protease–protease inhibitor pairs with exceptionally strong binding.

7.2. D-Regulators

Similar to R- and LC-Regulator, the implementation for D-Regulator I makes use of the interaction between NS3 protease and a suitable single-chain antibody (figure 9a). However, the antibody is solely expressed from a constitutive promoter in this case. As a second protease–protease inhibitor pair, we suggest the E. coli Lon protease and the phage T4 protease inhibitor PinA as discussed in our previous work [67]. For this purpose, a suitable degradation tag should be added to Y₁ and to avoid leaky integration due to endogenous Lon protease, a Lon-deficient E. coli strain, such as BL21(DE3) [77] should be used. Note that the latter protease–protease inhibitor pair can also be used for realizing the R-Regulator and LC-Regulator.

To realize the two annihilation reactions in D-Regulator II (figure 9b), we propose the use of sigma factors and anti-sigma factors as described previously [33,78]. Specifically, Z₃ could be the sigma factor SigW, which is constitutively expressed and mediates expression of SigF (Y₁). SigF mediates expression on the anti-sigma factor RsiW (Z₁), which binds to SigW. Analogous reactions are realized using SigM (Y₂), SigB (Z₄) and RsbW (Z₂).

The design for D-Regulator III may be more difficult to implement experimentally due to the requirement of a two-stage complex formation by three biomolecules (Z₁, Z₂ and Z₃) in addition to the requirement of Z₃ catalysing the production of Y₁ and Z₂ inhibiting Y₂. While it may be possible to achieve the desired behaviour of biomolecules using protein fusions and/or protein engineering, an alternative method to implement this design (as well as all the others) would be via molecular programming as discussed in the following section. In §S10 of the electronic supplementary material, we further discuss some challenges and limitations of such in vivo implementations accompanied by simulations based on more realistic (non-ideal) conditions.

7.3. Molecular programming implementation

In molecular programming, an abstract reaction network is realized by designing a concrete chemical reaction network using engineered molecules, so that the latter network emulates the kinetics of the former. At the edges of the abstract network, appropriate chemical transducers must be introduced to interface the abstract network with the environment. While such transducers are specific to each application, the core network is generic, and DNA (natural or synthetic) is commonly used to construct it. These systems are typically tested in vitro in controlled environments, with the eventual aim of embedding them in living cells, in synthetic cells [79], or in other deployable physical media. We refer to [21, §IV] for details of concrete synthetic DNA schemes in the context of biochemical regulation, and for literature overview. Suffice to say that all the reactions used in this paper can be systematically compiled into networks of synthetic molecules that well approximate the required mass action kinetics [80]. In particular, §S11 of the electronic supplementary material details the DNA strand-displacement realization of a bimolecular reaction $A + B ⇌ C + D$ . A collection of such reactions (and their unimolecular special cases) can then realize the chemical reaction networks used in this paper. Tools are available to simulate strand displacement systems, e.g. to evaluate their fidelity to the corresponding chemical reaction networks [21,22,80].

8. Discussion

In this paper, we address the challenge of regulating biomolecular processes with two outputs of interest which are, in the general case, co-dependent due to coupling interactions. This co-dependence means that disturbances applied to one of the outputs will also affect the other—each of the output species may be part of a separate, independent network and, by extension, be subject to different perturbations. Thus, we propose control schemes for efficient and robust manipulation of such processes adopting concepts based on both output steady-state coupling and decoupling. The proposed regulators describe biomolecular configurations with appropriate feedback interconnections which, under some assumptions, result in closed-loop systems where different types of output regulation can be achieved.

In particular, we present a variety of bio-controllers for regulating the ratio (R-Regulators) and linear combinations of the outputs (LC-Regulators) as well as each of the outputs individually (D-Regulators). At the core of their functioning lies a ‘hidden’ integral feedback action realized in suitable ways in order to meet the control objectives for each case. Integral control is one of the most widely used strategies in traditional control engineering since it guarantees zero control error and constant disturbance rejection at steady state. This is based on the fact that with this type of control, the existence of a positive/negative error, regardless of its magnitude, always generates an increasing/decreasing control signal. Essential structural components of these designs are production–inhibition loops [67] and/or annihilation reactions [32]. Moreover, to get a more practical insight, we consider a two-output biomolecular network with positive feedback coupling interactions. Treating the network as an open-loop system, we use our control designs to successfully manipulate its outputs under constant parameter perturbations and non-ideal conditions. At the same time, we discuss an alternative way of closing the loop in D-Regulator-II via a different controller species ‘pairing’. Although it may seem reasonable, we show that this feedback configuration leads to an unstable closed-loop system.

Assuming a biologically meaningful equilibrium, the proposed designs can be used to regulate arbitrary biological processes provided that the closed-loop topologies are asymptotically stable. We therefore anticipate that they will be useful for building complex pathways that robustly respond to environmental perturbations in synthetic biology applications. To this end, we extensively discuss ways of achieving local closed-loop asymptotic stability while, for R- and LC-Regulator, we also present specific sufficient conditions based on the concept of positive realness. Furthermore, we describe possible experimental implementations of all regulators using either biomolecular species in E. coli or molecular programming.

The regulation strategies presented in this work can be easily adapted to more complex networks to be controlled, than the one introduced in §4. In §S12 of the electronic supplementary material, we demonstrate the scalability of our control schemes to networks to be controlled with both monomolecular and bimolecular reactions, a high number of (strongly coupled) species and more than two outputs of interest. We also show that, in general, our control schemes do not require for actuator species to act directly on output species, as happens with the architectures discussed in the main text. In networks to be controlled with high number of accessible species, this can offer significant design flexibility, as different variations of our control schemes might be feasible. Our regulation strategies can be implemented through different biomolecular interconnections provided that the latter result in a stable closed-loop system with suitable memory variables and, by extension, in a desired steady-state output behaviour.

Biological networks are inherently stochastic due to the probabilistic nature of biomolecular interactions [8,81–84]. In the present study, we use deterministic mathematical analysis and simulations which offer a good approximation of the CRN dynamics when the biomolecular counts are sufficiently high. Thus, an interesting future endeavour would be to investigate the behaviour of our topologies within a stochastic mathematical framework examining, for instance, both the stationary mean and variance [64,84–88]. Another interesting extension of our work would be to study the non-local behaviour of our topologies. For example, the region of attraction for an equilibrium point of interest can be estimated via Lyapunov functions [89]. Additionally, treating those topologies as interconnections of suitably selected subsystems, dissipativity theory approaches based on storage functions can be used to assess the corresponding (local or global) stability [71,72].

Ethics

This work did not require ethical approval from a human subject or animal welfare committee.

Data accessibility

The programming codes supporting this work can be found at: https://github.com/emgalox/MIMO-bio-controllers.

The data are provided in electronic electronic supplementary material [90].

Authors' contributions

E.A.: conceptualization, formal analysis, methodology, software, writing—original draft, writing—review and editing; C.C.M.S.: conceptualization, methodology, writing—original draft, writing—review and editing; L.C.: conceptualization, methodology, supervision, writing—original draft, writing—review and editing; A.P.: conceptualization, methodology, supervision, writing—original draft, writing—review and editing.

All authors gave final approval for publication and agreed to be held accountable for the work performed therein.

Conflict of interest declaration

We declare we have no competing interests.

Funding

This work was supported by funding from the Engineering and Physical Sciences Research Council (EPSRC) (grant nos. EP/M002454/1 and EP/L016494/1) and the Biotechnology and Biological Sciences Research Council (BBSRC) (grant no. BB/M011224/1). C.C.M.S. was supported by the Clarendon Fund (Oxford University Press) and the Keble College De Breyne Scholarship. L.C. is supported by a Royal Society Research Professorship.

Footnotes

^†Present address: SecondCircle, Copenhagen, Denmark.

Electronic supplementary material is available online at https://doi.org/10.6084/m9.figshare.c.6767532.

Published by the Royal Society under the terms of the Creative Commons Attribution License http://creativecommons.org/licenses/by/4.0/, which permits unrestricted use, provided the original author and source are credited.

VIEW FULL TEXT DOWNLOAD PDF

This Issue

Journal of The Royal Society Interface cover image

August 2023

Volume 20Issue 205

Article Information

DOI:https://doi.org/10.1098/rsif.2023.0174
Published by:Royal Society
Online ISSN:1742-5662

History:

Manuscript received26/03/2023
Manuscript accepted06/07/2023
Published online02/08/2023

License:

Citations and impact

Keywords

Subjects

Your password must have 8 characters or more and contain 3 of the following:

Regulation strategies for two-output biomolecular networks

Abstract