Design of an embedded inverse-feedforward biomolecular tracking controller for enzymatic reaction processes.

Feedback control is widely used in chemical engineering to improve the performance and robustness of chemical processes. Feedback controllers require a 'subtractor' that is able to compute the error between the process output and the reference signal. In the case of embedded biomolecular control circuits, subtractors designed using standard chemical reaction network theory can only realise one-sided subtraction, rendering standard controller design approaches inadequate. Here, we show how a biomolecular controller that allows tracking of required changes in the outputs of enzymatic reaction processes can be designed and implemented within the framework of chemical reaction network theory. The controller architecture employs an inversion-based feedforward controller that compensates for the limitations of the one-sided subtractor that generates the error signals for a feedback controller. The proposed approach requires significantly fewer chemical reactions to implement than alternative designs, and should have wide applicability throughout the fields of synthetic biology and biological engineering.

Introduction

A major challenge in synthetic biology is to develop practically implementable design methods for the synthesis of feedback controllers that achieve reference tracking, i.e. force the output of a biomolecular process of interest to track desired changes in its concentration over time (Hsiao et al., 2015). The design of feedback controllers to control biochemical processes has received significant attention in the literature (see. e.g. Henson, 2003, Baldea et al., 2013), and the construction of synthetic control circuits has become a major focus of research in the new field of synthetic biology. Ideally, such circuits should be made up of well-defined modules consisting only of molecular reactions, in order to allow the realisation of embedded biomolecular control systems (Cosentino et al., 2016). A promising approach to facilitating the design of such circuits is provided by nucleic acid-based chemistry, wherein the design of biomolecular circuits can be done using abstract chemical reaction network (CRN) theory (e.g. Soloveichik et al., 2010), and then translated to deoxyribonucleic acid (DNA) using strand displacement reactions for implementations (Chen et al., 2013). A CRN is a collection of chemical reactions written in the formwhere γ is the reaction rate, the left-hand-side (LHS) of the reaction consists of reactants and the right-hand-side (RHS) of the reaction consists of products. Most of the chemical reactions considered in this paper are either unimolecular (i.e. one reactant on the LHS of (1)) or bimolecular (i.e. two reactants on the LHS of (1)). According to standard CRN theory (see e.g. Feinberg, 1986, Feinberg, 1988) a CRN with n species and m reactions can be represented by an ordinary differential equation (ODE) following generalised mass-action kinetic in the form ofwhere is the species concentration, is a function describing the reaction rates of the CRN, is the stoichiometric matrix that describe the dynamics of the species concentrations following their associated reaction rates, is the non-negative real number set, is the real number set, and n and m are positive integers.

In any reference tracking feedback system, it is imperative that an appropriate error signal can be computed such that the designed controller can take relevant control action to drive the process output towards the intended state. While such a requirement is trivial to satisfy in standard control theory, it is not in the context of CRN theory. This is because a two-sided subtractor (Fig. 1), which is an operator that is able to compute the difference between two input signals regardless of their relative magnitude, is yet to be realised using standard CRN's. For accurate reference tracking, the error e should be able to take both positive (r > y) and negative (r < y) values. Thus, the aforementioned constraint is a serious impediment to the design of functional biomolecular feedback control systems that will inevitably lead to poor quality reference tracking and potentially even instability.

Fig. 1

Subtraction operator.

To the best of the authors' knowledge, almost all previous designs for biomolecular subtractors using CRNs have resulted in only one-sided subtraction. We note that there is a literature on the design of half-subtractors or full-subtractors using digital logic gates realised using CRNs (see e.g. Xu et al., 2013, Lin et al., 2015), however, as our focus is on the design of analog biomolecular circuitry, we exclude this work from our discussion.

The subtraction operator used in our paper is based on the design presented in Buisman et al. (2009), which can be realised using a set of four chemical reactions (see Page 6 of Buisman et al. (2009)). The authors analysed the Jacobian matrix of the ODE associated with the subtraction operator and found that when the resulting subtraction is zero, the fixed point does not exist. Additionally, when the resulting subtraction is negative, the overall system diverges, as the fixed point is unstable. In view of this, the subtraction outputs a positive value when the magnitude of the first component is greater than the second component and zero when the condition is reversed. We further illustrate this point in Section 2.2 of our paper. Some other relevant results on biomolecular subtraction can be found in Salehi et al. (2016) and Song et al. (2016). The design considered in the former used the subtraction operator to realise a biomolecular computation of a Berstein polynomial and their subtraction is equivalent to the design of Buisman et al. (2009). The latter paper proposed frameworks to build operators using DNA strand displacement, and explicitly mentioned that their subtraction operator is one-sided. In Harris et al. (2015), the authors designed a feedback controller for gene expression regulation, which requires the use of a subtraction, but do not propose a detailed biomolecular implementation of the subtraction operator.

An alternative approach to the design of biomolecular subtraction operators can be found in Cosentino et al., 2013, Cosentino et al., 2016 and Bilotta et al., 2015, Bilotta et al., 2016. In this work, the authors developed subtractors that are used to compute the difference between two molecular fluxes, rather than two molecular concentrations. By satisfying certain conditions, and assuming all the reactions involve unitary stoichiometric coefficients with input fluxes constant, the output flux is shown to converge to the difference between the two input fluxes in an asymptotic manner. While the chemical reactions required to realise the subtraction operator in this framework are slightly different to the ones proposed in Buisman et al. (2009), the final ODE representation is exactly the same, and thus also yields a one-sided subtraction.

The only available partial solution to the problem mentioned above is to adopt the design framework proposed in Oishi and Klavins (2011). In this framework, each signal in the biomolecular circuit is implemented as the difference in the concentration of two chemical species. In this way, a two-sided subtraction operator can then be realised. As we show in the following section, however, this approach at least doubles the total number of chemical reactions required to implement the entire feedback circuit. This increase in the number of chemical reactions is highly undesirable as it presents a major challenge for wet lab implementation, and strongly limits the scalability of the design. Moreover, large numbers of chemical reactions potentially increases the probability of unwanted crosstalk interactions. For instance, a circuit whose implementation requires n molecular species will increase the potential bimolecular crosstalk interactions by n2. This has prompted researchers to look into ways to reduce crosstalk, such as requiring a certain number of mismatches for any two distinct recognition domains (see e.g. Qian and Winfree, 2011). Nonetheless, obtaining large numbers of well-behaved sequences with long domains is extremely difficult to achieve in practice.

In this paper, using the available framework for realising one-sided biomolecular subtraction using CRNs (see e.g. Buisman et al., 2009), we propose a design strategy that uses a model-inversion feedforward controller (see e.g. Devasia, 2002, Franklin and Powell, 2014) to circumvent the limitations of the feedback controller when using a one-sided subtraction operator. In addition, our controller design strategy also aims to utilise the minimal number of chemical reactions, to allow for a more scalable and feasible wet lab implementation.

Materials and methods

One-sided subtraction operator

To the best of authors’ knowledge, all current existing designs for biomolecular subtraction operators that utilise standard CRN theory can only implement one-sided subtraction. A comprehensive list of mathematical operators that can be implemented using CRN's, which includes the one-sided subtraction and its detailed analyses can be found in Buisman et al. (2009). Following the design of Buisman et al. (2009), the subtraction operator can be realised using the following abstract chemical reactions:where x and x are the two inputs, x is the resulting output, x is the intermediate states and γ is the reaction rate. Note that this one-sided subtraction operator is realised using four abstract chemical reactions. Using generalised mass-action kinetics (see e.g. Feinberg, 1986), these abstract chemical reactions can be represented by ODE's, where the corresponding ODE's for (2) are given byAt steady state, x = xx, leading to x = x − x. By analysing the Jacobian matrix of the ODE's relating to this subtraction operator, it has been shown in Buisman et al. (2009) that when the subtraction of two components results in either zero or a negative value, the fixed point does not exist or the system converges to an unstable fixed point, respectively. Hence, x = 0 when x < x and x = x − x > 0 when x ≥ x, making the subtraction one-sided. Other example of one-sided subtraction operators can be found in Cosentino et al. (2016) – in these cases the operators are used to compute the difference of molecular fluxes, rather than concentrations.

Two-sided subtraction operator

As mentioned in Section 1, the only known framework to date that offers a partial solution to the realisation of a two-sided subtraction operator can be found in Oishi and Klavins (2011). In that proposed framework, a signal, u is represented as a difference between two chemical species i.e. u : = u+ − u−. This representation results in the chemical species having positive and negative components. To illustrate how such representation can achieve a two-sided subtraction, we consider first how the summation operator can be realised using abstract chemical reactions. The abstract chemical reactions for the summation operator are given bywhere η is a reaction rate such that η ≫ γ. The corresponding ODEs are given bywhere at steady state (i.e. dx/dt = 0), x = x + x.

Now, for the subtraction operator, its abstract chemical reactions are given byNotice the difference of the superscripts + and − in the abstract chemical reaction compared to (4). The corresponding ODEs are given bywhere at steady state, x = x − x. Because the signal is represented as the difference of positive and negative components, x can be properly computed regardless of the relative magnitude of x and x. Note that both the summation and subtraction operators in (4), (6) are realised using nine abstract chemical reactions.

Now, note that the summation of two concentrations that is equivalent to (5) with positive signals could also have been realised by employing the following three abstract chemical reactions (Buisman et al., 2009): , and . This realisation does not require the signal to be represented using positive/negative components. Surprisingly, however, an equivalent way of representing the subtraction operator as in (7) seems not to exist, as there are no associated abstract chemical reactions to realise it. We shall further demonstrate this point below.

Consider the following two reactions: and . Their corresponding ODEs are dy/dt =+ γx and dy/dt =− γy respectively. Their final ODE expression can then be obtained by summing these two equations together, i.e. dy/dt = γ(x − y). Now, for the subtraction operator, the corresponding ODE is given by dy/dt = γ(x − x − y). We have already shown how dy/dt = γ(x − y) can be obtained. Therefore, we simply need another abstract chemical reaction that can achieve dy/dt =− γx. With the sign on the RHS of the ODE being negative, one would expect to write y on the LHS of the abstract chemical reaction. In addition, we require the multiplication of x with γ, which means x has to be on the LHS of the abstract chemical reaction as well. A natural first attempt would then be to write . However, a sum of reactants in an abstract chemical reaction leads to multiplication in the corresponding ODE, i.e. dy/dt =− γxy. If we are to move y to the RHS of the abstract chemical reaction, i.e. , its corresponding ODE would then be dy/dt =+ γx. Thus, there is no way to realise dy/dt =− γx using standard abstract chemical reactions. This is the reason why the positive/negative components formalism introduced by Oishi and Klavins (2011) is needed to realise a two-sided subtraction operator.

While the formalism proposed in Oishi and Klavins (2011) is able to realise two-sided subtraction, it also inevitably leads to a large increase in the number of abstract chemical reactions required, since it must be used to realise all operators and components in the circuit, even though in many cases there are simpler alternatives available (e.g. the summation operator). From the point of view of experimental implementation, it is highly desirable to keep the number of abstract chemical reactions required as small as possible, thus motivating our attempt to formulate a design strategy that can cope with the limitations of a one-sided subtraction operator and thus allow us to avoid using positive/negative components.

Enzymatic reaction processes

In this paper, we focus on the problem of controlling enzymatic reaction processes. Such processes are ubiquitous in cell biology, with some notable examples including protein phosphorylation by kinases, metabolic synthesis pathways and anaerobic fermentation of glucose to ethanol in yeast (see e.g. Segel, 1975, Galazzo and Bailey, 1990, Heinrich et al., 2002, Faeder et al., 2003, Hatakeyama et al., 2003, Charusanti et al., 2004, Gershon and Shaked, 2008).

The main function of enzymes is to act as the catalysts of biochemical systems (see e.g. Taylor et al., 2008). Disruption to the regulation of enzymatic reaction processes can therefore have significant adverse effect on biological function, and the ability to accurately control the dynamics of enzymatic reaction processes at their optimal levels is crucial to a wide range of natural and engineered biochemical systems.

The chemical reactions describing the enzymatic reaction process in its simplest form (see e.g. Segel, 1975, Gershon and Shaked, 2008) are given bywhere k, k and k are the process binding, catalytic and degradation rates respectively. x, x, x and x usually represent the substrate, enzyme, enzyme-substrate and product respectively. In a typical enzymatic reaction, the enzyme combines with its specific substrate to form the enzyme-substrate complex. After the reaction, this enzyme–substrate complex breaks up resulting in the associated product and the enzyme itself. This enzyme is unchanged and once separated from the complex is free to interact again with more substrate.

The corresponding ODE's for (8) are then given byHere, we assume that the total concentration of x and x represented as x = x + x is constant.

Theory/calculation

Derivation of inverse-feedforward controller

The basic aim of the inverse-feedforward controller is to invert the relevant dynamics of the process at steady state in order to provide accurate steady-state tracking of reference signals. From (9), at steady-state, (i.e. setting both dx/dt and dx/dt to zero), we have

After some algebraic manipulation, we getwhere α = kk, β = kkx and δ = kk. For reference tracking, the steady state value of the output, x should reach the reference signal, r. By substituting x = r and rewriting x = u into (11), we obtain the required control input signal such that the process is able to track the reference signal at steady state, which is given by

The control signal given in (12) essentially amounts to an open loop control strategy, in which the relevant dynamics of the process are inverted in order to achieve perfect steady-state tracking of reference signals. In practice, of course, the inversion will not be exact, due to model uncertainty, and there is also no way to control the transient dynamics of the closed loop system using this approach, which could lead to the presence of steady state errors. To address these limitations, the feedforward controller is used together with a classical proportional-integral (PI) feedback controller, as shown in Fig. 2. The main purpose of the feedback controller is to correct for errors introduced by model mismatch and to control the transient behaviour of the process. Thus, the final control signal acting on the process is given by the sum of the outputs of the feedforward and feedback controllers. Note that Subtraction I and Subtraction II shown in Fig. 2 are both one-sided subtraction operators.

Fig. 2

Block diagram configuration of the proposed tracking controller.

Abstract chemical reaction representation of the tracking controller

Here, we show how the tracking controller may be realised using abstract chemical reactions. As shown in Fig. 2, the feedforward controller requires two gain operators, one division operator and one (one-sided) subtraction operator. Following the design procedure for all those operators given in Buisman et al. (2009), the associated abstract chemical reactions for the feedforward controller part are given as follows:

[Gain, α:]where is α gain reaction rate.

[Gain, δ:]where is the δ gain reaction rate.

[Subtraction I:]where is Subtraction I reaction rate and x is the intermediate species involved in Subtraction I.

[Division:]where γ is the division reaction rate.

For the PI feedback controller, the associated abstract chemical reactions are given as follows:

[Integral gain:]where K is the integral gain.

[Proportional gain:]where K is the proportional gain and is the gain reaction rates.

[Summation I:]where is Summation I reaction rate.

The error resulting from model mismatch and transient effects is computed by Subtraction II, where the corresponding abstract chemical reactions are given by

[Subtraction II:]where is Subtraction II reaction rate and x is the intermediate species involved in Subtraction II.

The overall control signal to be applied to the process is the summation of the control signals from both the feedforward and feedback controllers, where the abstract chemical reactions are given by

[Summation II:]where is Summation II reaction rate.

With the control signal, x5 as the input to the process, the abstract chemical reactions for the process to be controlled are given by (8) with x = x5 and x = x6. Thus, the number of chemical reactions needed to realise the complete control circuit is 26. On the other hand, if the design framework that utilises two-sided subtraction as proposed in Oishi and Klavins (2011) is to be used to design just a standard PI feedback controller, 36 chemical reactions would be required in the circuit, an increase in complexity of 28%.

Ordinary differential equation representation of the tracking controller

Using generalised mass-action kinetics, the corresponding ODE's for all the abstract chemical reactions described in Section 3.2 are given by

[Feedforward controller:]

[Feedback controller:]

[Subtraction II:]

[Summation II:]

[Process:]with x = x + x is constant.

Results and discussion

Limitation of feedback control with one-sided subtraction

In this section, we first illustrate the effect of using a one-sided subtraction operator when a standard tracking feedback controller is used. Fig. 3(A) shows the configuration of a standard feedback system with a PI controller using a one-sided subtraction operator and the PI controller is tuned using the standard Ziegler–Nichol or Internal Model Control (IMC) methods (see e.g. Ogata, 2010, Rivera et al., 1986). From Fig. 3(B), at time 0 to 40,000s, the feedback controller attempts to track the reference value but as overshoot occurs, we have the situation that the reference value is smaller than the output value resulting in no control action before the undershoot where the situation is reversed. This leads to the oscillatory behaviour observed within that time span. At time 40,000s to 80,000s, the reference value is stepped down and we have the situation of the reference value being smaller than the output value. This means the error given by the one-sided subtraction is always zero and the PI controller exerts no control action, resulting in the large steady state error. The performance does not improve despite repeated attempts to tune the PI controller gain – one of the best tunings is shown in Fig. 3(C), where the PI controller is still unable to achieve proper reference tracking.

Fig. 3

(A) Block diagram configuration of a standard closed-loop feedback control system using a PI controller with a one-sided subtraction operator. (B) System response with PI feedback controller and one-sided subtraction operator with lower control gain, (C) with higher control gain.

Inverse-feedforward controller with one-sided subtraction

In Fig. 4, we show the results of repeating the above simulation with our inverse-feedforward controller architecture from Fig. 2. All the relevant parameters used in the simulation are summarised in Table 1. All units are assumed to be defined appropriately. The performance of the controller is good as it is able to track the reference signal properly. The contribution of the two controllers is as expected, where most of the control action is given by the feedforward controller while the PI controller is operative when dealing with transients. The parameters of the PI controller are obtained using standard tuning methods e.g., Ziegler–Nichol/IMC method (see e.g. Ogata, 2010, Rivera et al., 1986) followed by further fine tuning based on closed-loop simulations.

Fig. 4

System responses with the tracking controller. (A) Output and reference signals. (B) Control signals from inverse-feedforward (FF) and PI controllers.

Table 1

Parameters used in the closed-loop feedback control system.

Parameters	Values
Process
kr1	0.005
kr2	1.6
kr3	0.0008
xT	5.5
Inverse-feedforward controller
γ_G_α	1.0
γ_G_δ	1.0
γ_G_Sb_I	1.0
γ_G_D	1.0
PI controller
γ_G_K_P	1.0
KP	0.02
KI	2.5 × 10−8
γ_Sm_I	0.0004
Summation II and Subtraction II
γ_Sm_II	1.0
γ_Sb_II	3.0

Robustness and sensitivity analysis

Here, the robustness of the controller when implemented in closed-loop is investigated using Monte Carlo simulations, where all the parameters in (22), (23), (24), (25), (26) are randomly drawn from a uniform distribution and the above simulations are repeated. The number of Monte Carlo simulations required to obtain various levels of estimation uncertainty with known probability are determined following Chernoff bound (Vidyasagar, 1998). Following the guidelines given in Williams (2001), a total number of 1060 simulations are required to accomplish an accuracy level of 0.05 with confidence level of 99% (Vidyasagar, 1998, Menon et al., 2009). All the parameters are varied within ranges of 10% around their nominal values. Mathematically, this is written as θ(1 + 0.1Δ), where , K, K, k, k, k}, where i ∈ {1, 2}, j ∈ {1, 2, 3, 4, 5} and k ∈ {1, 2, 3}. Δ is a random number from the uniform distribution in [−1,1]. Note that all the associated reaction rates are split according to the number of chemical reactions in which they are involved.

The simulation results are shown in Fig. 5. The region shaded in grey is the output envelope encompassing all responses from 1060 Monte Carlo simulations for randomly perturbed parameters in the range of ±10% from the nominal values. The controller displays a good level of robust performance with no stability issues as a result of varying parameters. The simulation was repeated with perturbations up to ±100% (i.e. θ(1 + Δ) with Δ = 1) and no instability was observed.

Fig. 5

Robustness analysis of the tracking controller.

Despite the robustness analysis showing good performance with no stability issues, one notable result warranting further investigation is the presence of steady state error. It is known that any controller with integral action should eliminate steady state errors even in the presence of model uncertainty and this is not observed in our simulation. This is due to the effect of uncertainty on the biomolecular subtraction and summation operators. From (22), (23), (24), (25), (26), to achieve exact summation and subtraction, their associated reaction rates (γ and γ) are required to be identical. This can be seen by considering the following summation, i.e., dy/dt = γ1x1 + γ2x2 − γ3y. The summation is exact when γ1 = γ2 = γ3, where at steady state, y = x1 + x2. When the parameters are perturbed by uncertainty, however, this results in the summation reaction rates no longer being identical (γ1 ≠ γ2 ≠ γ3). Then, at steady state, y = (γ1/γ3)x1 + (γ2/γ3)x2, and the summation is no longer exact. Thus, these non-proper summation and subtraction operators are likely to result in the controller being unable to compute the correct control signal in response to the error, leading to the observed steady state error. To confirm this, we performed parameter sensitivity analysis, whereby each of the parameters shown in Table 1 are multiplied by a factor that ranges from 0.5 to 2.0 with increments of 0.1. Sensitivity is quantified by computing the relative steady state error at times 40,000s (step-up) and 80,000s (step-down) or mathematically,where is the output response subjected to parameter sensitivity. Table 2 shows the maximum relative steady state errors resulting from tracking the step-up and step-down of the reference signal.

Table 2

Parameter sensitivity analysis of the tracking controller. The maximum percentage relative steady state error has two rows for each parameter, where the upper and bottom rows denote Relative e and Relative e respectively.

Parameters	Rel. ess (%)	Parameters	Rel. ess (%)	Parameters	Rel. ess (%)
Subtraction II		Inverse-feedforward controller		Summation II
γ_Sb_II,1	29.58	γ_G_α,1	50.00	γ_Sm_II,1	50.00
	31.97		50.00		50.00
γ_Sb_II,2	20.95	γ_G_α,2	50.00	γ_Sm_II,2	0.25
	23.08		50.00		0.99
γ_Sb_II,3	0.25	γ_G_δ,1	0.25	γ_Sm_II,3	50.00
	0.99		0.00		50.00
γ_Sb_II,4	20.95	γ_G_δ,2	0.25	Process
	23.08		0.00	kr1	0.25
γ_Sb_II,5	20.95	γ_Sb_I,1	50.00		0.99
	23.08		50.00	kr2	0.00
PI controller		γ_Sb_I,2	0.25		0.00
KI	0.25		0.00	kr3	0.25
	0.99	γ_Sb_I,3	50.00		0.99
KP	0.00		50.00
	0.00	γ_Sb_I,4	0.25
γ_Sm_I,1	0.25		0.00
	0.99	γ_Sb_I,5	0.25
γ_Sm_I,2	0.00		0.00
	0.00	γ_G_D,1	50.00
γ_Sm_I,3	0.25		50.00
	0.99	γ_G_D,2	50.00
γ_G_K_P,1	0.00		50.00
	0.00
γ_G_K_P,2	0.00
	0.00

The results of the sensitivity analysis confirm that the steady state error is largely attributable to the Summation II and Subtraction II modules. These two operators are responsible for computing the total control signal to the process and the error to the PI controller respectively. The parameters associated with the feedforward controller are also highly sensitive. This is expected given that the feedforward controller operates by inverting the relevant dynamics of the process at steady state. Thus, any changes to the parameters within the feedforward controller could result in the computation of the incorrect control signal to negate the process dynamic, which subsequently leads to large steady state error. Care is thus required when specifying and implementing the reaction rates related to the controller and the summation and subtraction operators.

Simplified feedforward controller

The realisation of the inverse-feedforward controller involves a total of 26 abstract chemical reactions. While this produces a reduction of circuit complexity by 28% compared to using the framework proposed in Oishi and Klavins (2011), further reductions to the number of abstract chemical reactions would be highly desirable. This can be achieved using an alternative way to realise this feedforward controller, as follows. From (12), which we rewrite here,we note that r appears in both the numerator and denominator. This fractional representation can often be approximated by a polynomial to obtain a simplified representation. Taking the Taylor series expansion of u at r* = 0, we have the following:Neglecting the contribution of higher order terms and substituting r * =0 to (28), we getAs a remark, for r* ≠ 0, we can always perform a change of variables to ensure that the equilibrium is zero.

With this approximation, we have reduced the inverse-feedforward controller to requiring two gain operators, one summation operator and one polynomial operator. Substituting the relevant values of α, β and δ, which are related to the process parameters (i.e. k, k and k), we obtain u ≈ 0.029r + 0.000002r2. Given that the coefficient of the second term is very small, our final approximation of the inverse-feedforward controller is then given by u ≈ 0.029r. Hence, we have further reduced the complexity of the inverse-feedforward controller to requiring only one gain operator. With that, the block diagram configuration of the resulting simplified feedforward controller is shown in Fig. 6.

Fig. 6

Block diagram configuration with simplified feedforward controller.

Using the same variables as in the full feedforward controller case, the associated abstract chemical reaction for the simplified feedforward controller is given by

[Simplified feedforward controller:]where γ is the gain reaction rate and Γ = α/β. The associated ODE is given by

We repeat the reference tracking simulations as in the case of a full feedforward controller and the results are shown in Fig. 7. The results show that the performance of the simplified feedforward controller is similar to the case of a full feedforward controller. However, implementation of the simplified feedforward controller requires only 18 chemical reactions, a reduction of 31% and 50% from the full feedforward controller and framework proposed in Oishi and Klavins (2011) respectively.

Fig. 7

System responses to the tracking controller using simplified feedforward controller. Top: output, control and reference signals. Bottom: control signals from simplified feedforward (FF) and PI controllers.

We repeat the robustness analysis with the parameters , , where i ∈ {1, 2}, j ∈ {1, 2, 3} and k ∈ {1, 2, 3, 4, 5}. The result is shown in Fig. 8. Again, we observe similar performance compared to the full feedforward controller. In fact the envelope of possible responses for the simplified feedforward controller is smaller compared to the full feedforward controller, since the number of uncertain parameters is reduced. The sensitivity analysis is also repeated and the results are shown in Table 3.

Fig. 8

Robustness analysis of the tracking controller using the simplified feedforward controller.

Table 3

Parameter sensitivity analysis of the tracking controller using the simplified feedforward controller. The maximum percentage of relative steady state error has two rows, where the upper and bottom rows denote Relative e and Relative e respectively.

Parameters	Rel. ess (%)	Parameters	Rel. ess (%)
Subtraction II		Simplified feedforward controller
γ_Sb_II,1	29.58	γG,1	50.12
	31.97		50.25
γ_Sb_II,2	21.10	γG,2	50.12
	23.08		50.25
γ_Sb_II,3	0.50	Summation II
	0.99	γ_Sm_II,1	50.12
γ_Sb_II,4	21.10		50.25
	23.08	γ_Sm_II,2	0.50
γ_Sb_II,5	21.10		0.99
	23.08	γ_Sm_II,1	50.12
PI controller			50.25
KI	0.50	Process
	0.99	kr1	0.74
KP	0.50		0.99
	0.99	kr2	0.74
γ_Sm_I,1	0.50		0.99
	0.99	kr3	0.74
γ_Sm_I,2	0.50		0.99
	0.99
γ_Sm_I,3	0.50
	0.99
γ_G_K_P,1	0.50
	0.99
γ_G_K_P,2	0.50
	0.99

As in the case of the full feedforward controller, the most sensitive parameters are associated with the summation and subtraction operators as well as the gain of the simplified feedforward controller.

Design considerations and limitations for different biochemical processes

The dynamics of the enzymatic reaction process depend on k, k and k. From the experimental biology literature, it is clear that these three parameters can vary over several orders of magnitude - depending on the particular biochemical process in question values of k, k and k in the literature range from 10−9 to 106, 10−2 to 102 and 10−5 to 10−3 respectively (see e.g. Faeder et al., 2003, Hatakeyama et al., 2003, Charusanti et al., 2004). Such large variations in the process parameters can clearly significantly impact the design of the tracking controller and need to be taken into account in the design process. Here, we investigate the effect of each of these process parameters on the performance of the tracking controller using both the full and simplified feedforward controllers. We retain all the original parameters apart from the process parameters, which we will vary in our analysis. For each process parameter, we use the minimum and maximum values given above and consider their effect independently and jointly. Table 4 summarises the findings and the simulation results are shown in Fig. 9, Fig. 10.

Table 4

Effect of process parameters on reference tracking capability and remarks on the control design guidelines. For the reference tracking capability, Y denotes Yes and N denotes No.

Reference tracking capability (Y/N)	kr1	kr2	kr3	Design remarks
Inverse-feedforward controller
N	1 × 10−9	1 × 10−2	1 × 10−5
N	1 × 10−9	1 × 102	1 × 10−5
N	1 × 106	1 × 10−2	1 × 10−5
N	1 × 106	1 × 102	1 × 10−5
Y	1 × 10−9	1 × 10−2	1 × 10−3	Increase reaction rates γ
Y	1 × 10−9	1 × 102	1 × 10−3	Increase reaction rates γ
Y	1 × 106	1 × 10−2	1 × 10−3	Decrease KP and KI
Y	1 × 106	1 × 102	1 × 10−3	Decrease KP and KI
Simplified feedforward controller
N	1 × 10−9	1 × 10−2	1 × 10−5
N	1 × 10−9	1 × 102	1 × 10−5
N	1 × 106	1 × 10−2	1 × 10−5
N	1 × 106	1 × 102	1 × 10−5
Y	1 × 10−9	1 × 10−2	1 × 10−3	Increase KP and KI
Y	1 × 10−9	1 × 102	1 × 10−3	No change to existing parameters
Y	1 × 106	1 × 10−2	1 × 10−3	Decrease KP and KI
Y	1 × 106	1 × 102	1 × 10−3	Decrease KP and KI

Fig. 9

Effect of varying process parameters on reference tracking using inverse-feedforward controller. The notation ‘+’ and ‘−’ denotes respectively the maximum and minimum values of the process parameter. Red line: output response. Black line: reference (set-point). Magenta line: control signal from inverse-feedforward controller. Cyan line: control signal from PI controller. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of the article.)

Fig. 10

Effect of varying process parameters on reference tracking using simplified feedforward controller. The notation ‘+’ and ‘−’ denotes respectively the maximum and minimum values of the process parameter. Red line: output response. Black line: reference (set-point). Magenta line: control signal from simplified feedforward controller. Cyan line: control signal from PI controller. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of the article.)

In Fig. 9, Fig. 10, the signs (+) and (−) represent respectively the maximum and minimum values of k to k of the process. For example, the notation (+,−,−) means that we are analysing the scenario where k takes the maximum value of its range, i.e. 1 × 106, and k and k take the minimum value of their respective range, i.e. 1 × 10−2 and 1 × 10−5. Likewise, the notation (+,+,−) represents the analysis of the scenario where k and k take their maximum value of their respective range and k is at the minimum value of its range. In each of the subfigures ((A)–(H)), the top row represents the reference tracking capability of the output and the bottom row represents the control action given by both the feedforward and feedback control.

The obtained results show that the tracking controller using both full and simplified feedforward controller has difficulty in tracking reference inputs properly if the degradation term of the process k is very small (∼10−5). As shown in both Fig. 9, Fig. 10(A)–(D), the control action of PI (cyan line) stays at zero most of the time indicating that the subtraction is not working properly due to r < y. Now, with a small value of k, y will degrade slowly, meaning r < y for a longer time, hence the poor tracking (red line). On the other hand, if k is large (∼10−3), the time for which r < y is short as y degrades much faster, thus leading to accurate tracking as shown in both Fig. 9, Fig. 10(E)–(H).

For very small values of k, to achieve good reference tracking requires the reaction rates in the inverse-feedforward controller to be increased. For implementation using DNA-based chemistry, the ability to increase reaction rates could potentially be constrained by the physical binding property of the DNA, thus this point has to be taken into account when adjusting those reaction rates. For a large value of k, good reference tracking can be achieved by decreasing the PI controller gains K and K.

On the other hand, when the simplified feedforward controller is used, a good reference tracking can be achieved by just adjusting the PI controller gain. Finally, in the case when the value of k is small and the values of k and k are large, no adjustment is required to the controller parameters. In Fig. 9, Fig. 10(E) and (F), we note that despite the tracking controller being able to track the reference signal properly, potentially large control signals are required, whose feasibility needs to be checked during the design process with experimentalists.

Retroactivity

Here, we investigate how retroactivity affects the overall performance of the tracking controllers. In previous sections, we have assumed perfect modularity of the different elements in the closed-loop feedback control system shown in Fig. 2, Fig. 6. In other words, the dynamic responses are not affected by the interconnection of the components. While this assumption is widely made in the analysis of chemical reaction based systems, recent work (see e.g. Del Vechhio et al., 2008) has shown that such an assumption does not hold for many biomolecular feedback systems. As shown in Del Vechhio et al. (2008), the occurrence of different modules sharing the same molecular species is common, and this sharing of species can affect the overall dynamics of the processes upon their interconnection. To quantify the way that interconnection of two modules alters their dynamics with respect to their behaviour in isolation, the concept of retroactivity has been introduced (Del Vecchio and Jayanthi, 2008, Jayanthi et al., 2013). For the process considered here, it should be noted that no retroactivity effects are present due to the interconnection of modules involving unimolecular reactions. For example, for the closed-loop feedback control systems in Fig. 2, Fig. 6, retroactivity does not affect the interconnection of the module Subtraction II with the PI controller.

On the other hand, an interconnection of two modules, where one module comprises unimolecular reactions while the other module comprises bimolecular reactions, will feature a unidirectional retroactivity. Again, using an example in the context of the closed-loop feedback control system shown in both Fig. 2, Fig. 6, retroactivity affects the interconnection of the module Summation II with the process (see (8), (21)). To take into account the effect of retroactivity, the ODE representation of Summation II must consider the chemical reactions downstream of the process (see (25)). Thus the ODE representation of the Summation II module including the effect of retroactivity is given as follows:

We repeat the simulations in Section 4.2 using both the full and simplified feedforward controller taking into account the effects of retroactivity, and the results are shown in Fig. 11, Fig. 12. The results show that retroactivity does not significantly affect the performance of the overall closed-loop system other than to introduce a very small steady state error. The reason for this steady state error can be explained as the following. The retroactivity affects the input to the process, x5. At the same time, x5 is also the resulting control signal from both the feedforward and PI controllers. As both the controllers do not have information regarding this error resulting from retroactivity, these controllers are unable to provide the correct control signal to react to the effect of retroactivity, thus leading to the observed steady state error. Nevertheless, our results show that the performance of the tracking controller is not significantly affected by retroactivity, allowing us to undertake designs based on the assumption of modularity (as is the case with standard control applications).

Fig. 11

System responses with retroactivity using inverse-feedforward controller. Top: output and reference signals. Bottom: control signals from inverse-feedforward (FF) and PI controllers.

Fig. 12

System responses to retroactivity using simplified feedforward controller. Top: output, control and reference signals. Bottom: control signals from simplified feedforward (FF) and PI controllers.

Implementation of CRN via DNA strand displacement chemistry

Finally, we briefly provide an overview of the way in which the type of designs presented in this paper may be implemented experimentally. Recently, there have been a number of studies that described how synthetic circuits composed of abstract chemical reactions may be readily implemented using DNA-based chemistry (see e.g. Seelig et al., 2006, Soloveichik et al., 2010, Qian and Winfree, 2011, Chen et al., 2013). In Chen et al. (2013), the authors highlighted that chemical reactions can serve as a programming language for the design of DNA-based chemistry synthetic circuits. Therefore, circuit components synthesised using DNA that can be expressed mathematically can be derived from biologically synthesised plasmids. Consequently, in principle this enables the in vitro implementation of those circuits. The advantage of utilising DNA-based chemistry in the design lies in the ease of implementation, as the design depends on the choice of relevant sequences following the standard Watson-Crick pairing (i.e. adenine-thymine and guanine-cytosine or A-T and G-C). As a result, design frameworks and software tools that allow synthetic biomolecular circuitry described using CRN's to be readily implemented using DNA strand displacement (DSD) chemistry have been developed in recent years (see e.g. Seelig et al., 2006, Soloveichik et al., 2010, Qian and Winfree, 2011, Chen et al., 2013).

In particular, it has been shown in Soloveichik et al. (2010) that unimolecular (one reactant at the LHS of the chemical reaction) and bimolecular (two reactants at the LHS of the chemical reaction) chemical reactions can be compiled in DSD chemistry to accomplish the intended behaviour of their designed biomolecular circuit. All designs presented in this work are all made up of either unimolecular or bimolecular chemical reaction, making this framework suitable for the designed controller to be implemented experimentally. Here, we briefly describe their proposed framework and refer interested readers to Soloveichik et al. (2010) for details.

Consider the following bimolecular DSD reaction,where k is the binding reaction rate and k is the unbinding reaction rate of the DNA strand. The reaction commences when the invader strand Q binds in a standard Watson–Crick complementary manner to the toe-hold domain of strand X. As the binding occurs, portions of the strand of X are displaced whereby this separation results in the product Y and waste W. This partially double stranded product Y can then bind with to the toe-hold domain of other DNA strands and complexes for subsequent reactions. The rate of the overall reaction can be controlled by varying the binding and unbinding rates, k and k. Quite often, different DNA strands do not interact directly with one another, and hence auxiliary species with sufficiently large amount are required to facilitate their interaction. With the inclusion of auxiliary species, the DSD implementation for the respective unimolecular (X → Y + Z) and bimolecular (X + Y → Z) reactions are given by,andwhere G, O, T, L, H and B are the auxiliary species with appropriate initial concentrations C. The partial strand displacement rate is given by q = k/C and q is the maximum strand displacement rate.

Since the abstract chemical reactions describing the tracking controller presented here are made up exclusively of unimolecular and bimolecular reactions, the framework introduced in Soloveichik et al. (2010) enables our design to be implemented via DNA chemistry. Nevertheless, since the introduction of auxiliary species further increases the number of reaction needed for DSD implementations, it is imperative that circuit designs utilise as few reactions as possible. Hence, our results showing the feasibility of a radically simplified inverse-feedforward tracking controller offer strong potential for future wet lab implementations of embedded biomolecular control systems.

Conclusions

In this paper, we have shown for the first time how a controller architecture for implementing reference tracking, based on the use of an inverse-feedforward controller, can be adapted to the specific context of embedded biomolecular feedback systems. Our proposed approach circumvents a fundamental limitation of CRN-based systems from the point of view of tracking control, i.e. the inability of standard CRNs to implement a two-sided subtraction operator. The only currently available solution to this problem proposed by Oishi and Klavins (2011) requires the adoption of a non-standard modelling and design framework that represents signals as the difference in the concentration of two chemical species, an approach that significantly increases the complexity of the design process and doubles the total number of chemical reactions needed to implement a given circuit.

The complexity of the control problem considered here is significantly higher than in previous studies – almost all previous studies consider either a static or a very simple first order process that avoids the possibility of the reference value being larger than the output value (Cosentino et al., 2016) or employing control strategies that can reduce the duration of the reference value being smaller than the output value (Briat et al., 2016) or consider only a very simple static process with no dynamics (Yordanov et al., 2014). Using the existing standard realisation of one-sided subtraction, our approach of using inverse feedforward combined with feedback control produces highly accurate and robust reference tracking. By exploiting the biochemical structure of the feedforward controller, we are able to reduce its complexity using a Taylor series approximation. The resulting simplified feedforward control not only achieves similar performance to the full feedforward controller, but also utilises far fewer chemical reactions. The reduction in the number of chemical reactions is important, as it will significantly facilitate the experimental implementation of the proposed design in DNA-based chemistry either in vitro or in vivo.

Conflict of interest

The authors declare that they have no competing interests.