Exploiting the dynamic properties of covalent modification cycle for the design of synthetic analog biomolecular circuitry
 Mathias Foo^{1}Email authorView ORCID ID profile,
 Rucha Sawlekar^{1} and
 Declan G. Bates^{1}
DOI: 10.1186/s1303601600361
© The Author(s) 2016
Received: 3 June 2016
Accepted: 19 October 2016
Published: 14 November 2016
Abstract
Background
Cycles of covalent modification are ubiquitous motifs in cellular signalling. Although such signalling cycles are implemented via a highly concise set of chemical reactions, they have been shown to be capable of producing multiple distinct inputoutput mapping behaviours – ultrasensitive, hyperbolic, signaltransducing and thresholdhyperbolic.
Results
In this paper, we show how the set of chemical reactions underlying covalent modification cycles can be exploited for the design of synthetic analog biomolecular circuitry. We show that biomolecular circuits based on the dynamics of covalent modification cycles allow (a) the computation of nonlinear operators using far fewer chemical reactions than purely abstract designs based on chemical reaction network theory, and (b) the design of nonlinear feedback controllers with strong performance and robustness properties.
Conclusions
Our designs provide a more efficient route for translation of complex circuits and systems from chemical reactions to DNA strand displacementbased chemistry, thus facilitating their experimental implementation in future Synthetic Biology applications.
Keywords
Covalent modification cycle Chemical reaction networks Analog synthetic biomolecular circuits Linear and nonlinear operators Feedback control systems Synthetic biology applicationsBackground
The emerging field of Synthetic Biology [1, 2] has provided a number of recent examples of the successful construction of digital circuits in cells, including biomolecular transistors [3, 4] and logic gates [5–7]. A major open challenge associated with such digital circuitry is to ensure robust functionality of the designed circuit in the presence of interactions with the host cell [8, 9]. One possible solution to this problem is to develop analog circuitry, [10], since in general analog designs require far fewer devices to carry out a given computation at the moderate precision needed in cells, resulting in lower resource requirements and a reduced metabolic burden on the cell [11, 12]. There is also a growing realisation that in many cases the use of purely digital designs can fail to exploit the hybrid approach to information processing used by cells, which often combines digital computation with analog processing of signals with different levels of gradation [13]. The ability to process signals with different levels of gradation will be a key requirement for the development of more complex synthetic circuits with advanced monitoring and control capabilities.
Transduction of external perturbations via signalling cascades is a classical examples of analog signal processing in the cell, [14–16]. One of the most ubiquitous motifs seen in cell signalling cascades is the cycle of covalent modification [17, 18], examples of which include phosphorylation/dephosphorylation of enzymes [19], DNA methylation [20] and monoclonal antibodies [21]. The covalent modification cycle is implemented via a highly concise set of chemical reactions, which have been shown in [17, 22] to potentially exhibit highly sigmoidal inputoutput characteristics, generating socalled ultrasensitive or hyperbolic responses (see also [23, 24]). In [25], the authors systematically examine the steadystate and dynamic responses of covalent modification cycles to time varying perturbations and demonstrate the existence of two additional types of regimes, termed signaltransducing and thresholdhyperbolic, giving a total of four distinct mapping regimes (see Figure 2 of [25]).
Thus, by modifying only the reaction rates for the chemical reactions governing the covalent modification cycle, four highly distinct inputoutput mapping behaviours can be obtained. Here, we show how this flexibility can be exploited for the design of synthetic analog biomolecular circuitry. From an engineering design point of view, this versatile inputoutput mapping property is highly attractive, as different combinations of these inputoutput mappings can be used to design circuits that can perform many different types of computation, information processing, and control. From an experimental implementation point of view, the extremely concise set of chemical reactions used in covalent modification is also very advantageous, since it facilitates circuit designs with fewer reactions and components, hence easing their implementation using DNAbased chemistry.
Methods
Chemical reactions underlying the covalent modification cycle motif
Whereas signals in systems and control theory can have both positive and negative values, this is not the case for biomolecular concentration, as they can only take nonnegative values. This issue can be addressed following the framework suggested in [26]. Here, we present a summarised version of this framework, for full details the reader is referred to [26].
A signal x is represented as the difference in concentration of two chemical species, x ^{+} and x ^{−}. Thus, x ^{+} and x ^{−} are, respectively, the positive and negative components of x such that x=x ^{+}−x ^{−}. From an implementation point of view, e.g. using DNAbased chemistry as proposed in [27], x ^{+} and x ^{−} can each represent a DNA strand and the final concentration of x can be recovered using x=x ^{+}−x ^{−}. As an example, consider an initial DNA strand, x ^{+} with 10nM concentration present in the system. Then, another DNA strand, x ^{−} with 20nM concentration is added to the system and the resulting signal x becomes negative. In terms of the underlying chemical reactions, ensuring proper implementation of x=x ^{+}−x ^{−}, requires a fast annihilation reaction between x ^{+} and x ^{−}, i.e. \(x^{+} + x^{} \xrightarrow {\eta } \emptyset \) with η is the annihilation rate.
In this work, we adopt this formalism as it enables the realisation of negative signals. While there are alternative formalism available from standard chemical reaction network theory (see e.g. [28]) these cannot deal with negative signals or realise twosided subtraction operators, which are required, for example in feedback control design to compute the difference between two signals, where the resulting difference can result in either a positive or negative signal.
with total concentration x _{ total }=x _{ e }+x _{ C2} is constant and these ODE’s produce four distinct mapping regimes for the different choices of reaction rates, [25] (see Fig. 1(d)).
Implementation of chemical reactions via DNAbased chemistry
A number of recent studies have described how synthetic circuits composed of abstract chemical reactions may be readily implemented in DNAbased chemistry (see e.g. [29, 30, 34, 35]). As highlighted in [35], these chemical reactions can serve as a programming language for designing synthetic circuits based on DNAbased chemistry. Mathematically expressed components of circuits designed using DNA can be derived from biologically synthesised plasmids, in principle enabling the in vitro implementation of those circuits. A particular advantage of employing DNAbased chemistry lies in the ease of implementation, given that the design relies on the choice of relevant sequences following the wellknown WatsonCrick (i.e. adeninethymine and guaninecytosine) pairing. In [30], it is shown that unimolecular and bimolecular chemical reactions can be compiled into DNA strand displacement (DSD)based chemistry to achieve the desired behaviour of the considered biomolecular circuit. Here, we present a summarised version of the framework and refer interested readers to [30] for more details.
A simple bimolecular DSD reaction can be described by the following reaction,
where δ _{ b } and δ _{ ub } are the binding and unbinding rate of the DNA strand respectively. This reaction begins when Q, called the invader strand, binds with X in a complementary manner at the toehold domain of X. When this binding takes place, parts of the strand of X are displaced and this disengagement results in product Y and waste, R. The partially double stranded product, Y can then bind with other toehold domains of other DNA complexes. Usually, the toehold domain has short nucleotides (610 nucleotides) to ensure the reaction is fast and reliable. By varying the value of δ _{ ub } and δ _{ b }, one can control the rate of reaction. Specifically, the binding and unbinding rates can be changed by changing the length of the DNA strands. To further elaborate, the reaction rates from chemical reactions can be mapped into these binding and unbinding rates of the DNA reactions following the framework of [30]. Based on these binding and unbinding rates, the designer can determine the length of DNA strands required to achieve that reaction.
where G, O, T, L, H, and B are auxiliary species with appropriate initial concentrations C _{ max }, q=δ/C _{ max } is the partial strand displacement rate and q _{ max } is the maximum strand displacement rate.
The set of chemical reactions governing the covalent modification cycle is made up exclusively of unimolecular and bimolecular reactions, and thus, following the framework shown above, all of the circuits described in the remainder of this paper can be implemented via nucleic acids. Note, however, that the experimental challenges associated with building such circuitry increase with the number of reactions that are required for a given circuit. It is therefore imperative that circuit designs utilise as few reactions as possible, in order to ease experimental implementation and maximise scalability of the resulting synthetic systems.
Results and discussions
Computing nonlinear operators
Here, we present analog biomolecular circuit designs, based on the chemical reactions underlying the covalent modification cycle, that can compute three important nonlinear operators – the logarithm of arbitrary base, the signum function and the absolute value of a signal. We show that these designs achieve a dramatic reduction in circuit complexity when compared with designs based on purely abstract chemical reaction networks.
Computing a logarithm of arbitrary base
where l is the order of the series. The larger the order l is, the better the approximation, but the higher the complexity of the circuit. In this paper, we choose l=10 as this order allows us to compute the logarithm of numbers up to 10.
Multiplication: Consider the following multiplication: y=u _{1} u _{2}. At steadystate, the set of abstract chemical reactions that realise this operator is given by \(u_{1}^{\pm } + u_{2}^{\pm } \xrightarrow {\gamma _{M}} u_{1}^{\pm } + u_{2}^{\pm } + y^{+}\), \(u_{1}^{\pm } + u_{2}^{\mp } \xrightarrow {\gamma _{M}} u_{1}^{\pm } + u_{2}^{\mp } + y^{}\), \(y^{\pm } \xrightarrow {\gamma _{M}} \emptyset \) and \(y^{+} + y^{} \xrightarrow {\eta } \emptyset \), where γ _{ M } is the multiplication reaction rate. In total, 7 abstract chemical reactions are required to realise the multiplication operator, whose ODE is given by \(\frac {dy}{dt} = \gamma _{M}(u_{1}u_{2}  y)\).
where γ _{ P } is the power exponent reaction rate. In total, 7(n−1) abstract chemical reactions are required to realise the power exponent operator, and the corresponding ODE is given by \(\frac {dy}{dt} = \gamma _{P}\Big \{\Big (\prod _{l=1}^{n} u\Big)  y\Big \}\).
With l=10, we require 13 summation and subtraction operators, one multiplication operator, 10 power exponent operators with exponents 3,5,…,21 and 12 gain operators. This results in a total of \(13(7) + 1(7) + \sum _{n=3,5,\ldots,21} 7(n1) + 12(5) = 928\) abstract chemical reactions. To compute the logarithm of arbitrary base, as shown in Fig. 2(c), we require one more Φ that computes the second natural logarithm and one each for the subtraction, multiplication and gain operator. Thus, this circuit requires a total of 2(928)+1(7)+1(7)+1(5)=1875 abstract chemical reactions, which makes it completely intractable from an experimental point of view.
The huge number of abstract chemical reactions required to implement the circuit described above prompted us to seek alternative, more efficient, designs. We note that the characteristic of a natural logarithm resembles the hyperbolic regime of the covalent modification cycle (see Fig. 1(d)i) thus, making this regime potentially useful for computing the natural logarithm. Interestingly, this response is not governed by the order of the series approximation. Thus, as long as one can obtain the appropriate reaction rates for k _{1} to k _{4}, we can compute the natural logarithm. Moreover, this approach requires only 14 abstract chemical reactions. To compute the logarithm of arbitrary base using this approach, we replace the Φ block in Fig. 2(c) with the covalent modification cycle reactions in the hyperbolic regime, as shown in Fig. 2(b). This results in a total of 2(14) + 7 + 7 + 5 = 47 abstract chemical reactions, a 97 % reduction in circuit complexity.
Simulation results for computing log105 and log28 are shown in Fig. 2(Di) and (Dii) respectively. To implement the hyperbolic response, reaction rates are k _{1}=0.22 /M/s, k _{2}=0.43 /s, k _{3}=1.03 /M/s, k _{4}=35.10 /s, and x _{ e }=1 M.
For both approaches, the computed logarithms are close to the actual value, however the circuit based on the covalent modification cycle motif is significantly faster in settling to the correct steadystate value, even though it uses far fewer chemical reactions. An alternative approach for the biological computation of logarithms has been designed and implemented in [13]. This approach utilises transcriptional regulation, which requires a host cell, while our approach can be implemented in cellfree conditions (e.g. using DNA Strand Displacement (DSD) framework). Moreover, [13] considers the computation of the natural logarithm, while our approach enables the computation of logarithms of arbitrary base.
Computing the signum function
Now, the characteristics of a signum function closely resemble the ultrasensitive regime of the covalent modification cycle. Thus, as shown in Fig. 3(b), the corresponding chemical reactions can also be utilised to compute the signum function. While the ultrasensitive regime alone could be used to compute the signum function, we go a step further by appending a gain operator, K that can be used for scaling purposes. With the introduction of the gain operator, the total number of reactions required to compute the signum function is now only 14 + 1(5) = 19 reactions, a reduction in complexity of 77 %.
In the simulation results shown in Fig. 3(c), the input to our signum function is in the range of 10^{−6} while the output is in the range of 10^{−3}. We could have specified the signum function to produce outputs in the range of ±1 instead of ±1×10^{−3}, however, there could arise situations where such a response is not achievable due to the physical constraints on the system. In view of this, the gain component can be used to scale the output from ±10^{−3} to ±1. This ultrasensitive response is obtained with k _{1}=5,000 /M/s, k _{2}=5 /s, k _{3}=50 /M/s, k _{4}=0.05 /s and x _{ e }=0.1μ M.
Fig. 3(c) shows the simulation results of the signum function computed using both circuits. The input is changed every 10,000 seconds starting with 2.0×10^{−6}, 5.5×10^{−6}, −3.0×10^{−6} and −0.5×10^{−6}. The simulation result shows that while both circuits are able to compute the signum function accurately, the much simpler circuit based on covalent modification exhibits smaller transient overshoots in response to changes in the value of the input signal.
Computing the absolute value of a signal
We shall now illustrate how the remaining two regimes of the covalent modification cycle, i.e. signaltransducing and thresholdhyperbolic, can be utilised to compute the absolute value of a signal, for which the block diagram is shown in Fig. 4(b). The thresholdhyperbolic regime has a deadzone, (i.e. a nonresponsive region given any input signal) followed by a hyperboliclike region. To compute the absolute value, the deadzone range must not respond to input signals, u that are strictly negative and then respond to input signals that are nonnegative in a linear manner. To achieve this, note that in the thresholdhyperbolic regime, any hyperbolic response has an “almost linear” region when the input signal is small. By taking advantage of this property, we can ensure that our required thresholdhyperbolic regime has a linear instead of hyperbolic response after the deadzone region. On the other hand, the signaltransducing regime has a linear region followed by a plateau region. This makes this regime suitable for responding only to nonpositive input signals and not to strictly positive input signals. We also introduce two gain components, K _{1} and K _{2} for scaling purposes. By combining these two regimes with two gain components and one subtraction operation, 2(14) + 2(5) + 7 = 45 reactions are required to compute the absolute value, a reduction in circuit complexity of 21 %.
Since the thresholdhyperbolic response has a linear response with unity gradient, there is no requirement for the gain block, K _{1} or equivalently, K _{1}=1. To achieve this thresholdhyperbolic response, we set k _{1}=0.0027 /M/s, k _{2}=16,640 /s, k _{3}=0.043 /M/s, k _{4}=0.008 /s and x _{ e }=3.5 M. For the signaltransducing response, suppose that due to the limitations imposed by the system, a unity gradient of the linear response cannot be achieve, resulting in the gradient of the linear response to be 20. In this case, the gain component is set to K _{2}=1/20. To achieve this signaltransducing response, we set k _{1}=5 /M/s, k _{2}=100 /s, k _{3}=5 /M/s, k _{4}=630 /s and x _{ e }=1.8 M. Fig. 4(c) shows the simulation results for six different input signals, u ranging from +1 M to +6 M. At time 10,000 s, these input signals, u are switched to their negative counterpart ranging from 1 M to 6 M.
From Fig. 4(c), we see that the circuit using the configuration shown in Fig. 4(a) performs very well. For the circuit using the covalent modification cycle, the performance is also excellent, although when u=±1 M and ±6 M some deviations are observed. This is because the thresholdhyperbolic and signaltransducing responses achieved are not a perfect match to the ideal responses, as shown in Fig. 4(d).
Remarks on choosing the reaction rates of the covalent modification cycle
In all our illustrations above, the reaction rates (i.e. k _{1} to k _{4}) of the covalent modification cycle are obtained via numerical optimisation using the MATLAB ^{ Ⓡ } function ‘fminsearch’. The numerical optimisation aims to find the reaction rates (within biologically valid ranges) that minimise the difference between the desired mapping regime of the covalent modification cycle and the one obtained from Eq. (2).
Designing nonlinear controllers
Here, we illustrate how the chemical reactions underlying the covalent modification cycle can be used for the design of nonlinear biomolecular feedback controllers.
Controller descriptions
The modules involved in the biomolecular feedback controllers are as follows.
PI controller: The classical PI controller considered here is designed according to the methodology of [26]. The PI controller is made up of one integrator, one proportional gain and one summation operator. These three submodules require a total of 15 abstract chemical reactions to implement as follows:
Integrator: \(e^{\pm } \xrightarrow {K_{I}} e^{\pm } + n^{\pm }\) and \(n^{+} + n^{} \xrightarrow {\eta } \emptyset \), where K _{ I } is the integral gain of the PI controller and η is the annihilation rate.
Proportional gain: \(e^{\pm } \xrightarrow {\gamma _{K}K_{P}} e^{\pm } + m^{\pm }\), \(m^{\pm } \xrightarrow {\gamma _{K}} \emptyset \) and \(m^{+} + m^{} \xrightarrow {\eta } \emptyset \), where K _{ P } is the proportional gain of the PI controller, γ _{ K } is the gain reaction rate.
Summation junction: \(m^{\pm } \xrightarrow {\gamma _{Sm}} m^{\pm } + u^{\pm }\), \(n^{\pm } \xrightarrow {\gamma _{Sm}} n^{\pm } + u^{\pm }\), \(u^{\pm } \xrightarrow {\gamma _{Sm}} \emptyset \), and \(u^{+} + u^{} \xrightarrow {\eta } \emptyset \), where γ _{ Sm } is the summation reaction rate.
The tuning of this PI controller involves adjusting K _{ P }, K _{ I } and the reaction rates γ _{ K } and γ _{ Sm }.
CMC controller: The chemical reactions required to implement the CMC controller are given by Eq. (1) with e=x _{ in } and u=x _{ out }. We choose values for the CMC controller’s reaction rates that place it in its signaltransducing inputoutput mapping regime, which closely resembles the steadystate inputoutput mapping of a PI Controller (see Fig. 5(b)). Note that the CMC controller requires 14 reactions to implement, 1 fewer than the PI controller.
Closedloop system and ODE approximations
Comparative performance of the two controllers is evaluated for two biomolecular processes  the first a simple first order linear process and the second a more complex second order nonlinear process. The abstract chemical reactions for both processes are given by
Linear process: \(u^{\pm } \xrightarrow {k_{p1}} u^{\pm } + y^{\pm }\), \(y^{\pm } \xrightarrow {k_{p2}} \emptyset \) and \(y^{+} + y^{} \xrightarrow {\eta } \emptyset \), where k _{ p1} and k _{ p2} are the catalysis and degradation rates of the process.
Nonlinear process: \(u^{\pm } + p^{\pm } \xrightarrow {k_{r1}} q^{+}\), \(u^{\pm } + p^{\mp } \xrightarrow {k_{r1}} q^{}\), \(q^{\pm } \xrightarrow {k_{r2}} y^{\pm } + p^{\pm }\), \(y^{\pm } \xrightarrow {k_{r3}} \emptyset \) and \(y^{+} + y^{} \xrightarrow {\eta } \emptyset \), where p and q are intermediate species involved in the second order process reaction. k _{ r1}, k _{ r2} and k _{ r3} are respectively the binding, catalytic and degradation rates of the process. For the subtraction operator, the abstract chemical reactions are \({r}^{\pm } \xrightarrow {\gamma _{Sb}} {r}^{\pm } + {e}^{\pm }, {y}^{\pm } \xrightarrow {\gamma _{Sb}} {y}^{\pm } + {e}^{\pm }, {e}^{\pm } \xrightarrow {\gamma _{Sb}} \emptyset \) and \({e}^{+} + {e}^{} \xrightarrow {\eta } \emptyset \).
Using generalised mass action kinetics, the ODEs corresponding to the abstract chemical reactions employed in the modules of the feedback control system are given as:
According to the formalism of [26], the gain and summation operators used in the PI controller require multiple identical reaction rates to be used in their sets of abstract chemical reactions. However, implementing this requirement in an experimental setting is unlikely to be feasible, as experimental biologists are rarely able to specify the reaction rates of chemical reactions exactly. Additionally, in practice, as highlighted in [26], unregulated chemical devices or leaky expressions could potentially affect production and degradation rates and subsequently alter the behaviour of the designed component. To investigate these issues, we perform a formal robustness analysis of both controllers, focussed on the effect of uncertainties in the implemented reaction rates on the closedloop stability and performance properties of the feedback system.
Performance analysis of controllers with linear process
To analyse the performance and robustness of the closedloop responses achieved by the feedback controllers with the linear process, step response tests and Monte Carlo simulations are performed, respectively. For the Monte Carlo simulations, all the parameters are randomly drawn from a uniform distribution. The number of Monte Carlo simulations required to achieve various levels of estimation uncertainty with known probability are calculated using the wellknown Chernoff bound [38]. Following the guidelines provided in [39], an accuracy level of 0.05 and a confidence level of 99 % are chosen for the Monte Carlo simulation analysis, which requires a total number of 1060 simulations [38, 40]. To investigate the effect of different levels of uncertainty we vary the parameters within ranges of 20 %, 50 %, 100 % and 120 % around their nominal values. Mathematically, we have p(1+Δ P(x)), where p∈{γ _{ Sb1},γ _{ Sb2},γ _{ Sb3},γ _{ K1},γ _{ K2},γ _{ Sm1},γ _{ Sm2},γ _{ Sm3},K _{ I },K _{ P },k _{1},k _{2},k _{3}, k _{4},k _{ p1},k _{ p2}}, P(x) is the probability distribution and Δ∈{0.2,0.5,1.0,1.2}. Note that we split reaction rates γ _{ Sb }, γ _{ K } and γ _{ Sm } according to the number of chemical reactions in which they are involved.
In our simulations, a step change in the concentration of the reference species, r from 0 M to 1 M occurs at time 0 s and the purpose of the controller is to ensure that the process output reaches this new desired concentration. As quantitative measures of the control system performance, the step response characteristics, which comprise rise time, t _{ r }, settling time, t _{ s }, percentage of overshoot, M _{ OV } and steadystate error, e _{ ss } are used [41]. For good closedloop performance, it is desirable to achieve a small t _{ r }, t _{ s } and M _{ OV } as well as having e _{ ss }=0. As a benchmark for comparison, we first calculate the step response characteristics without parameter uncertainty. Hereafter, we refer to these as the set of results for the nominal system. The parameters for the nominal system in the required abstract chemical reactions are:
Linear process: k _{ p1}=0.1 /s, k _{ p2}=0.1 /s.
PI controller: γ _{ Sb1}, γ _{ Sb2}, γ _{ Sb3}=0.4 /s, γ _{ Sm1}, γ _{ Sm2}, γ _{ Sm3}=0.8 /s, γ _{ K1}, γ _{ K2}=0.0004 /s, K _{ P }=1 and K _{ I }=0.045.
CMC controller: k _{1}, k _{3}=0.00185 /M/s, k _{2}, k _{4}=0.5 /s. x _{ p }+u+x _{ C1}+x _{ C2}=27.5 M and x _{ e }+x _{ C2}=0.033 M.
Subtractor dynamics: γ _{ Sb1}, γ _{ Sb2}, γ _{ Sb3}=0.4 /s.
Step response characteristics and worstcase parameter ranges for PI and CMC controllers + linear process
PI controller  
Characteristics  Nominal  Δ=0.2  Δ=0.5  Δ=1.0  Δ=1.2 
t _{ r } (s)  29  44  75  157  173 
t _{ s } (s)  96  113  175  499  Unstable 
M _{ OV } (%)  9.14  22.83  52.25  114.17  Unstable 
e _{ ss } (M)  0.00  0.19  0.49  0.91  Unstable 
Parameters  Nominal  Δ=0.2  Δ=0.5  Δ=1.0  Δ=1.2 
γ _{ Sb1} (/s)  0.400  0.431–0.476  0.532–0.595  0.475–0.791  0.584–0.599 
γ _{ Sb2} (/s)  0.400  0.401–0.471  0.400–0.584  0.413–0.569  0.428–0.875 
γ _{ Sb3} (/s)  0.400  0.401–0.473  0.414–0.593  0.466–0.721  0.412–0.853 
K _{ I }  0.045  0.048–0.054  0.048–0.061  0.052–0.086  0.058–0.085 
K _{ P }  1.000  1.016–1.165  1.130–1.359  1.137–1.549  1.159–1.367 
γ _{ K1} (/s) [10 ^{−3}]  0.400  0.436–0.477  0.424–0.515  0.434–0.674  0.505–0.795 
γ _{ K2} (/s) [10 ^{−3}]  0.400  0.403–0.466  0.401–0.454  0.473–0.666  0.570–0.683 
γ _{ Sm1} (/s)  0.800  0.809–0.948  0.863–1.099  0.827–1.544  0.825–1.410 
γ _{ Sm2} (/s)  0.800  0.835–0.943  0.904–1.012  0.849–1.205  1.152–1.548 
γ _{ Sm3} (/s)  0.800  0.832–0.958  0.823–1.140  0.841–1.536  0.853–1.279 
k _{1} (/s)  0.100  0.101–0.116  0.106–0.142  0.111–0.174  0.127–0.208 
k _{2} (/s)  0.100  0.101–0.114  0.102–0.139  0.103–0.199  0.123–0.211 
CMC controller  
Characteristics  Nominal  Δ=0.2  Δ=0.5  Δ=1.0  Δ=1.2 
t _{ r } (s)  29  37  65  89  117 
t _{ s } (s)  97  116  155  202  353 
M _{ OV } (%)  10.12  25.3  44.63  60.55  75.00 
e _{ ss } (M)  0.00  0.18  0.46  0.92  1.12 
Parameters  Nominal  Δ=0.2  Δ=0.5  Δ=1.0  Δ=1.2 
γ _{ Sb1} (/s)  0.400  0.403–0.476  0.450–0.595  0.412–0.781  0.626–0.857 
γ _{ Sb2} (/s)  0.400  0.401–0.477  0.406–0.580  0.407–0.752  0.403–0.745 
γ _{ Sb3} (/s)  0.400  0.403–0.466  0.402–0.594  0.577–0.726  0.425–0.806 
k _{ b1} (/M/s) [10 ^{−2}]  0.185  0.186–0.221  0.202–0.274  0.199–0.342  0.224–0.365 
k _{ b2} (/s)  0.500  0.514–0.565  0.516–0.683  0.621–0.737  0.509–0.747 
k _{ b3} (/M/s) [10 ^{−2}]  0.185  0.187–0.213  0.187–0.245  0.242–0.354  0.230–0.318 
k _{ b4} (/s)  0.500  0.516–0.596  0.679–0.723  0.553–0.851  0.662–1.008 
k _{1} (/s)  0.100  0.100–0.119  0.104–0.148  0.130–0.199  0.101–0.174 
k _{2} (/s)  0.100  0.101–0.114  0.105–0.150  0.106–0.198  0.109–0.209 
The performance of the two nominal closedloop systems is rather similar, which reflects the fact that the CMC controller is designed to reproduce the steadystate inputoutput mapping of the original PI controller. Interestingly, however, we can clearly see a significantly improved robustness of the system when the CMC controller is used. With the PI controller, the closedloop system become unstable when Δ=1.2, while for the CMC controller, the closedloop system becomes unstable only when Δ=1.8, showing that the CMC controller is able to tolerate more than a 50 % larger variability in the values of the reaction rates in the underlying chemical reactions.
Performance analysis of controllers with nonlinear process
We proceed to analyse how the two controllers fare in controlling a more complex secondorder nonlinear process. The same step test and Monte Carlo simulations are carried out, with the parameters for the nominal system in the required abstract chemical reactions given as:
Nonlinear process: k _{ r1}=0.00005 /M/s, k _{ r2}=1.6 /s, k _{ r3}=0.0008 /s, with the total concentration constrained so that p+q=5.5 M.
PI controller: γ _{ Sb1}, γ _{ Sb2}, γ _{ Sb3}, γ _{ Sm1}, γ _{ Sm2}, γ _{ Sm3}, γ _{ K1}, γ _{ K2}=0.0004 /s, K _{ P }=0.65 and K _{ I }=0.3.
CMC controller: k _{1}==0.0000055 /M/s, k _{3}=0.000018 /M/s, k _{2}=12.50 /s, k _{4}=140 /s, x _{ p }+u+x _{ C1}+x _{ C2}=66 M and x _{ e }+x _{ C2}=0.00012 M.
Subtractor dynamics: γ _{ Sb1}, γ _{ Sb2}, γ _{ Sb3}=0.4 /s.
Step response characteristics and worstcase parameter ranges for PI and CMC controllers + nonlinear process
PI controller  
Characteristics  Nominal  Δ=0.2  Δ=0.5  Δ=1.0 
t _{ r } (s)  11,139  18,959  31,673  Unstable 
t _{ s } (s)  26,304  44,138  48,482  Unstable 
M _{ OV } (%)  2.42  24.88  44.44  Unstable 
e _{ ss } (M)  0.00  0.19  0.48  Unstable 
Parameters  Nominal  Δ=0.2  Δ=0.5  Δ=1.0 
γ _{ Sb1} (/s) [10 ^{−3}]  0.400  0.4450.480  0.4920.596  0.4720.730 
γ _{ Sb2} (/s) [10 ^{−3}]  0.400  0.4020.475  0.4000.453  0.5330.735 
γ _{ Sb3} (/s) [10 ^{−3}]  0.400  0.4270.467  0.4330.591  0.4020.473 
K _{ I } [10 ^{−3}]  0.300  0.3010.353  0.3100.420  0.4860.504 
K _{ P }  0.650  0.6730.771  0.7900.908  0.8061.282 
γ _{ K1} (/s) [10 ^{−3}]  0.400  0.4320.461  0.4460.573  0.4980.747 
γ _{ K2} (/s) [10 ^{−3}]  0.400  0.4010.422  0.4460.586  0.4040.740 
γ _{ Sm1} (/s) [10 ^{−3}]  0.400  0.4270.479  0.4240.587  0.6410.791 
γ _{ Sm2} (/s) [10 ^{−3}]  0.400  0.4230.462  0.4690.553  0.6150.730 
γ _{ Sm3} (/s) [10 ^{−3}]  0.400  0.4090.478  0.4180.539  0.4010.431 
k _{ r1} (/M/s) [10 ^{−4}]  0.500  0.5090.594  0.5360.734  0.7370.945 
k _{ r2} (/s)  1.600  1.6331.865  1.7492.232  1.8602.843 
k _{ r2} (/s) [10 ^{−3}]  0.800  0.8120.904  0.8191.102  0.8440.891 
CMC controller  
Characteristics  Nominal  Δ=0.2  Δ=0.5  Δ=1.0 
t _{ r } (s)  11,147  15,501  25,753  30,838 
t _{ s } (s)  28,848  28,324  42,494  49,196 
M _{ OV } (%)  2.84  13.12  26.25  56.37 
e _{ ss } (M)  0.00  0.19  0.46  0.98 
Parameters  Nominal  Δ=0.2  Δ=0.5  Δ=1.0 
γ _{ Sb1} (/s) [10 ^{−3}]  0.400  0.4260.479  0.5340.597  0.5420.798 
γ _{ Sb2} (/s) [10 ^{−3}]  0.400  0.4000.478  0.4040.511  0.4030.798 
γ _{ Sb3} (/s) [10 ^{−3}]  0.400  0.4060.457  0.4250.578  0.4030.633 
k _{ b1} (/M/s) [10 ^{−5}]  0.550  0.5670.644  0.5770.808  0.6191.075 
k _{ b2} (/s)  12.50  12.6414.75  16.5717.84  16.9917.95 
k _{ b3} (/M/s) [10 ^{−4}]  0.180  0.1820.207  0.2030.238  0.1960.274 
k _{ b4} (/s)  140.00  144.48163.09  165.38205.25  143.15278.30 
k _{ r1} (/M/s) [10 ^{−4}]  0.500  0.5030.593  0.5230.712  0.5380.807 
k _{ r2} (/s)  1.600  1.6351.893  1.8391.950  1.7892.861 
k _{ r2} (/s) [10 ^{−3}]  0.800  0.8030.943  0.8041.162  0.8081.525 
The performance of the two nominal closedloop systems are rather similar, which again reflects the fact that the CMC controller is designed to reproduce the steadystate inputoutput mapping of the original PI controller. The closedloop system with the CMC controller retains closedloop stability up until Δ=1.6, again demonstrating a significantly higher level of robustness than exhibited by the linear PI controller.
Flexible inputoutput mapping improves robustness
This intriguing observation leads us to ask how the gradient of this mapping of steadystate inputoutput signals is related to the robustness of the controller? Here, the gradient is associated with the straight line equation, y=mx+c, where m is the gradient and c is the intersection of the yaxis given that the mapping of the steadystate inputoutput signals are made up of a straight line.
As Eq. (13) is in scalar form, the overall process is stable if the real part of the eigenvalue of A (i.e. k _{ p1} K−k _{ p2}) is less than 0, hence, the following condition, \(K < \frac {k_{p2}}{k_{p1}}\) must hold. In other words, if the controller gain, K is less than the ratio of the process parameters k _{ p2} to k _{ p1}, we have a stable system. In our simulation, the process parameters of the nominal system are k _{ p1},k _{ p2}=0.1, thus, for the system to be stable, we require \(K < \frac {k_{p2}}{k_{p1}} = 1\). Looking at Fig. 8(b), a zoomedin version using M _{ OV } as illustration, the gradient of both the controllers’ inputoutput mapping are less than 1 (i.e. ≈0.34); the closedloop system is stable. Note that for the nominal system, both the controllers’ inputoutput mappings are very similar, as expected, since the CMC controller was designed to reproduce the PI controller’s steadystate inputoutput mapping.
We now consider the effect of increasing levels of variability in the values of the parameters in the chemical reactions implementing the feedback control system. For the PI controller, at Δ=1.2, the process parameters change from k _{ p1}=0.100→0.208 and k _{ p2}=0.100→0.124. Thus, the ratio \(\frac {k_{p2}}{k_{p1}}\) changes from 1 →0.596. Likewise, from Fig. 8(b), we see the gradient of the PI controller’s steadystate inputoutput mapping changes to 1.213 \(> \frac {k_{p2}}{k_{p1}} = 0.596\), which accounts for the observed unstable behaviour.
On the other hand, the change of gradient for the CMC controller is smaller compared to the PI controller. At Δ=1.2, the process parameters change from k _{ p1}=0.1→0.174 and k _{ p2}=0.1→0.109, leading the ratio \(\frac {k_{p2}}{k_{p1}} \) to change from 1 →0.628. However, the gradient of the CMC controller’s steadystate inputoutput mapping changes to 0.588 \(< \frac {k_{p2}}{k_{p1}} = 0.628\), thus preserving the stability of the system.
As the process is now nonlinear, the notion of eigenvalue no longer applies while the notion of stability for a nonlinear system is also more mathematically involved and beyond the scope of this paper. However, we can informally explain the difference in the robustness of both controllers by extending our arguments on the ‘gradient’ of the steadystate inputoutput mapping, as was done for the linear process. Figure 9(b) shows the nominal and worstcase deviation in the inputoutput mapping for both controllers at Δ=1.0.
From Fig. 9(c), we can see that despite both controllers having very similar mapping of inputsignal signals for the nominal system, when subjected to parameter uncertainty, the gradient of the PI controller’s steadystate inputoutput mapping becomes steeper and subsequently affects the stability of the system. On the other hand, not only does the CMC controller’s inputoutput mapping show a smaller change in response to uncertainty, it becomes more hyperbolic. The CMC controller’s innate ability to achieve hyperbolic behaviour seems to be able to prevent the adverse effect of parameter uncertainty, as it enables the gradient of its inputoutput mapping when subjected to parameter uncertainty to remain small. In this particular case, we show through an extensive analysis that the change in the mapping behaviour of the CMC controller from the linear signaltransducing regime to the hyperbolic regime actually acts to improve the robustness of the CMC controller.
Conclusions
In this paper, we have shown how the set of chemical reactions underlying the covalent modification cycle motif may be used to design a range of analog biomolecular circuits for computation, information processing and control. By exploiting the four distinct inputoutput mapping behaviours of the covalent modification cycle, we have designed biomolecular circuits to compute complex nonlinear operators and implement nonlinear feedback controllers. Our design approach results in a dramatic reduction in circuit complexity compared to the use of purely abstract reactions from standard chemical reaction network theory, and requires far fewer chemical reactions to be implemented experimentally. Our designs also demonstrated significantly greater robustness to variability in circuit parameters that will inevitably arise in experimental implementations of synthetic circuitry. Given the range of inputoutput mappings that can be produced by the set of chemical reactions underlying the covalent modification cycle, it is likely that they could be used to efficiently design many other types of operators and controllers. As the chemical reactions concerned are all represented either in unimolecular or bimolecular form, the resulting circuits can then be readily implemented using DNAbased chemistry either in vitro or in vivo for future Synthetic Biology applications.
Abbreviations
 DNA:

Deoxyribonucleic acid
 DSD:

DNA strand displacement
 ODE:

Ordinary differential equation
 PI:

proportionalintegral
 CMC:

Covalent modification cycle
Declarations
Acknowledgements
Not applicable.
Funding
The publication of this article is funded by the Biotechnology and Biological Sciences Research Council (BBSRC), UK via research grants BB/M017982/1 and the Engineering and Physical Sciences Research Council (EPSRC), UK.
Availability of data and materials
Data sharing not applicable to this article as no datasets were generated or analysed during the current study.
Authors’ contributions
MF, RS participated in data analysis. MF, DGB conceived and designed the study. All authors wrote the manuscript. All authors gave final approval for publication.
Competing interests
The authors declare that they have no competing interests.
Consent for publication
Not applicable.
Ethics approval and consent to participate
Not applicable.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver(http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
Authors’ Affiliations
References
 Purnick PEM, Weiss R. The second wave of synthetic biology: from modules to systems. Nat Rev Mol Cell Biol. 2009; 10:410–22.View ArticleGoogle Scholar
 Kahl LJ, Endy D. A survey of enabling technologies in synthetic biology. J Biol Eng. 2013; 7(1):1–19.View ArticleGoogle Scholar
 Liu Y, Chakrabartty S. Factor graphbased biomolecular circuit analysis for designing forward error correcting biosensors. IEEE Trans Biomed Circ Syst. 2009; 3:150–9.View ArticleGoogle Scholar
 Woo SS, Kim J, Sarpeshkar R. A cytomorphic chip for quantitative modeling of fundamental biomolecular circuits. IEEE Trans Biomed Circ Syst. 2015; 9:527–42.View ArticleGoogle Scholar
 Sharma P, Roy S. Alloptical biomolecular parallel logic gates with bacteriorhodopsin. IEEE Trans NanoBioscience. 2004; 3(2):129–36.View ArticleGoogle Scholar
 GoñiMoreno A, Amos M. A reconfigurable nand/nor genetic logic gate. BMC Syst Biol. 2012; 6(1):1–11.View ArticleGoogle Scholar
 Marchisio MA. In silico design and in vivo implementation of yeast gene boolean gates. J Biol Eng. 2014; 8(1):1–15.View ArticleGoogle Scholar
 Chen BS, Wu CH. A systematic design method for robust synthetic biology to satisfy design specifications. BMC Syst Biol. 2009; 3(1):1–18.View ArticleGoogle Scholar
 Weiße AY, Oyarzún DA, Danos V, Swain PS. Mechanistic links between cellular tradeoffs, gene expression, and growth. Proc Natl Acad Sci USA. 2015; 112(9):1038–47.View ArticleGoogle Scholar
 Teo JJY, Song SW, Sarpeshkar R. Synthetic biology: a unifying view and review using analog circuits. IEEE Trans Biomed Circ Syst. 2015; 9:453–74.View ArticleGoogle Scholar
 Sarpeshkar R. Analog versus digital: extrapolating from electronics to neurobiology. Neural Computation. 1998; 10(7):1601–38.View ArticleGoogle Scholar
 Sauro HM, Kim K. Synthetic biology: it’s an analog world. Nature. 2013; 497(7451):572–3.View ArticleGoogle Scholar
 Daniel R, Rubens JR, Sarpeshkar R, Lu TK. Synthetic analog computation in living cells. Nature. 2013; 497(7451):619–24.View ArticleGoogle Scholar
 Xiong L, Schumaker KS, Zhu JK. Cell signaling during cold, drought and salt stress. Plant Cell. 2002; 14:165–83.View ArticleGoogle Scholar
 Poli G, Leonarduzzi G, Biasi F, Chiarpotto E. Oxidative stress and cell signalling. Curr Med Chem. 2004; 11(9):1163–82.View ArticleGoogle Scholar
 Gomperts BD, Kramer IM, Tatham PER. Signal transduction. Massachusetts: Academic Press; 2004.Google Scholar
 Goldbeter A, Koshland Jr DE. An amplified sensitivity arising from covalent modification in biological systems. Proc Natl Acad Sci USA. 1981; 78(11):6840–4.MathSciNetView ArticleGoogle Scholar
 Ventura AC, Jiang P, Wassenhove LV, Del Vecchio D, Merajver SD, Ninfa AJ. Signaling properties of a covalent modification cycle are altered by a downstream target. Proc Natl Acad Sci USA. 2010; 107(22):10032–7.View ArticleGoogle Scholar
 Krebs EG, Beavo JA. Phosphorylationdephosphorylation of enzymes. Annu Rev Biochem. 1979; 48:923–59.View ArticleGoogle Scholar
 Miller CA, Sweatt JD. Covalent modification of DNA regulates memory formation. Neuron. 2007; 53:857–69.View ArticleGoogle Scholar
 Rodwell JD, Alvarez VL, Lee C, Lopes AD, Goers JWF, King HD, Powsner HJ, McKearn TJ. Sitespecific covalent modification of monoclonal antibodies: in vitro and in vivo evaluations. Proc Natl Acad Sci USA. 1983; 83:2632–6.View ArticleGoogle Scholar
 Goldbeter A, Koshland Jr DE. Ultrasensitivity in biochemical systems controlled by covalent modification. interplay between zeroorder and multistep effects. J Biol Chem. 1984; 259(23):14441–7.Google Scholar
 Ferrell Jr JE, Ha SH. Ultrasensitivity part I: Michaelian responses and zeroorder ultrasensitivity. Trends Biochem Sci. 2014; 39:496–503.View ArticleGoogle Scholar
 Heinrich R, Neel BG, Rapoport TA. Mathematical models of protein kinase signal transduction. Molecular Cell. 2002; 9:957–60.View ArticleGoogle Scholar
 GomezUribe C, Verghese GC, Mirny LA. Operating regimes of signaling cycles: statics, dynamics, and noise filter. PLoS Comput Biol. 2007; 3(12):246.MathSciNetView ArticleGoogle Scholar
 Oishi K, Klavins E. Biomolecular implementation of linear I/O systems. IET Syst Biol. 2011; 5:252–60.View ArticleGoogle Scholar
 Yordanov B, Kim J, Petersen RL, Shudy A, Kulkarni VV, Philips A. Computational design of nucleic acid feedback control circuits. ACS Synth Biol. 2014; 3(8):600–16.View ArticleGoogle Scholar
 Buisman HJ, ten Eikelder HMM, Hilbers PAJ, Liekens AML. Computing algebraic functions with biochemical reaction networks. Artif Life. 2009; 15:1–15.View ArticleGoogle Scholar
 Seelig G, Soloveichik D, Zhang DY, Winfree E. Enzymefree nucleic acid logic circuits. Science. 2006; 314(5805):1585–8.View ArticleGoogle Scholar
 Soloveichik D, Seelig G, Winfree E. DNA as a universal substrate for chemical kinetics. Proc Natl Acad Sci. USA. 2010; 107(12):5393–8.View ArticleGoogle Scholar
 Lakin MR, Youssef S, Polo F, Emmott S, Phillips A. Visual DSD: a design and analysis tool for DNA strand displacement systems. Bioinformatics. 2011; 27(22):3211–3.View ArticleGoogle Scholar
 Horn F. Necessary and sufficient conditions for complex balancing in chemical kinetics. Arch Ration Mech Anal. 1972; 49:172–86.MathSciNetView ArticleGoogle Scholar
 Klonowski W. Simplifying principles for chemical and enzyme reaction kinetics. Biophys Chem. 1983; 18:73–87.View ArticleGoogle Scholar
 Qian L, Winfree E. A simple DNA gate motif for synthesizing largescale circuits. J. Royal Soc. Interface. 2011. doi:10.1098/rsif.2010.0729.
 Chen YJ, Dalchau N, Srinivas N, Phillips A, Cardelli L, Soloveichik D, Seelig G. Programmable chemical controllers made from DNA. Nat Nanotechnol. 2013; 8(10):755–62.View ArticleGoogle Scholar
 Abramowitz M, Stegun IA. Handbook of mathematical functions. Appl Math Ser. 1966; 55:62.MathSciNetGoogle Scholar
 Foo M, Sawlekar R, Kim J, Bates DG, Stan GB, Kulkarni VV. Biomolecular implementation of nonlinear system theoretic operators. In: Proceedings of European Control Conference: 2016. p. 1824–31.
 Vidyasagar M. Statistical learning theory and randomised algorithms for control. IEEE Control Syst. 1998; 18(6):69–85.View ArticleGoogle Scholar
 Williams PS. A Monte Carlo dispersion analysis of the x33 simulation software. In: Proceedings of AIAA Conference on Guidance, Navigation and Control: 2001. p. 6–9.
 Menon PP, Postlethwaite I, Bennani S, Marcos A, Bates DG. Robustness analysis of a reusable launch vehicle flight control law. Control Eng Prac. 2009; 17(7):751–65.View ArticleGoogle Scholar
 Franklin GF, Powell JD, EmamiNaeini A. Feedback control of dynamic systems, 7th edn. New York: Pearson; 2015.MATHGoogle Scholar
 Friedland B. Control system design: an introduction to statespace methods. New York: McGrawHill; 2005.MATHGoogle Scholar