 Research
 Open Access
 Published:
Exploiting the dynamic properties of covalent modification cycle for the design of synthetic analog biomolecular circuitry
Journal of Biological Engineering volume 10, Article number: 15 (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.
Background
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.
For conciseness, we use the notation x ^{±}→x ^{±}+y ^{∓} to represent a pair of reactions, x ^{+}→x ^{+}+y ^{−} and x ^{−}→x ^{−}+y ^{+}. In [26] and [27] it was shown how sets of abstract chemical reactions implemented using this formalism can be used to design biomolecular circuits that compute a number of linear operators, such as scalar gains, summation/subtraction and integration. Circuits designed using this framework can then be implemented experimentally using DNAbased chemistry, as discussed in [29–31]. Using the mathematical formalism of [26], we can represent a covalent modification cycle of an enzymesubstrate pair by the following set of 14 abstract chemical reactions:
where all the x represent biomolecular species and k _{ j } (j=1,⋯,4) are the reaction rates. The covalent modification cycle described by Eq. (1) is illustrated in Fig. 1(a), and operates in the following manner. \(x_{p}^{\pm }\) represents the inactive component, which is associated with \(x_{in}^{\pm }\), (i.e., the enzyme kinase), forming an intermediate species, \(x_{C1}^{\pm }\) with reaction rate k _{1}. The intermediate species then produces the active component, \(x_{out}^{\pm }\) and the remaining unused \(x_{in}^{\pm }\) with reaction rate k _{2}. This reaction from inactive state to active state is called the forward reaction. Likewise, the active component, \(x_{out}^{\pm }\) is associated with x _{ e } (i.e. the enzyme phosphatase) with reaction rate k _{3} and produces another intermediate species, \(x_{C2}^{\pm }\). This intermediate species produces inactive component, \(x_{p}^{\pm }\) and the remaining unused x _{ e } with reaction rate k _{4}. This reaction from active state to inactive state is called the backward reaction. Note that as x _{ e } is externally introduced, it is not split into positive and negative components.
Using mass action kinetics, (see e.g. [32, 33]), the nonlinear ordinary differential equations (ODE’s) for the covalent modification cycle are given by:
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.
Given that different DNA strands do not interact directly with each other, their interaction is normally mediated by auxiliary species that are usually present in large amounts. The framework of approximating abstract chemical reactions to DSD presented in [30] considers the DNA implementation for unimolecular and bimolecular reactions, where these two types of reactions can be represented by Eqs. (4) and (??) respectively.
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
Consider the operation c=log_{ b } a, i.e. computing the logarithm of a to the base b. This logarithm can be computed through the change of logarithm base, i.e. \(c = \frac {\ln a}{\ln b}\), where ln denotes the natural logarithm. In other words, c can be realised as a ratio of lna and lnb. Several numerical methods exist to compute the natural logarithm. The most commonly used method is to use Taylor series, but this method accurately computes the logarithm of numbers only within the range 0<x<2. A more efficient method to compute the natural logarithm for x≥2 is based on the area hyperbolic tangent series approximation [36]. Thus, using such approximation, the natural logarithm can be computed as follows:
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.
The block diagram of a circuit that can compute the natural logarithm using the area hyperbolic tangent series approximation of order l=10 is shown in Fig. 2(a). This circuit uses a combination of several linear and nonlinear operators; summation, subtraction, gain, multiplication and power exponent, each of which may be implemented using a number of abstract chemical reactions. For details of the chemical reactions used in computing the linear summation, subtraction and gain operators, see [26], (the numbers of reactions needed are 7, 7, and 5 respectively). For the nonlinear multiplication and power exponent operators, their abstract chemical reactions (modified from [28] and [37]) are given as follows.
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)\).
Power exponent: Consider the following exponent: y=u ^{n}, where n is a positive integer. At steadystate, the set of abstract chemical reactions that realise this operator is given by
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
A general signum function outputs +1 for any positive real valued input and 1 for any negative real valued input. Mathematically, this is represented as X=sgn(X)X. The signum function can be designed using the circuit shown in Fig. 3(a). Details of the three abstract chemical reactions used to realise the integrator used in this circuit are provided in [26]. In this circuit, the input signal is first squared before taking its square root. The square root operator can be implemented using the wellknown NewtonRaphson method, which requires two subtraction operators, two multiplication operators, one power exponent operator of order 3, one integrator operator and one gain operator, resulting in a total of 2(7) + 2(7) + 1(14) + 1(3) + 1(5) = 50 reactions. After taking the square root, the reciprocal of this signal is calculated, and this is then multiplied by the input signal, which results in either +1 or 1 depending on the sign of the input signal. With a power exponent operator of order 2, one square root operator, one subtraction operator, one gain operator and two multiplication operators, the total number of reactions required to compute the signum function is 1(7) + 1(50) + 1(7) + 1(5) + 2(7) = 83.
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
The block diagram of a circuit that can compute the absolute value of a signal is shown in Fig. 4(a). As shown, the input signal is first squared before taking its square root. We have introduced the NewtonRaphson method previously for computing the square root and thus it can be seen that this circuit requires a total of 7 + 50 = 57 reactions.
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
Figure 5(a) shows the block diagram of a biomolecular feedback control circuit. In industrial control systems, the most commonly used controller is the linear proportionalintegral (PI) controller, and this type of controller has been successfully implemented for biomolecular systems using DNAbased chemistry in previous studies [26, 27]. Interestingly, the signaltransducing regime of the covalent modification cycle resembles the steadystate inputoutput mapping of the PI controller (Fig. 5(b)). Here, we compare the performance and robustness properties of a nonlinear ‘covalent modification cycle’ (CMC) controller, designed to operate in its signaltransducing regime, with those of a classical PI controller.
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:
Subtraction operator:
PI controller:
CMC controller:
Linear process:
Nonlinear process:
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.
The step response characteristics for both the nominal systems are tabulated in Table 1. For each of the analysed uncertainty sets, the worstcase values returned by Monte Carlo simulation for each of the step response characteristics and its associated parameter set are shown. Note that a range of parameters is given here as the parameter set associated with each worstcase characteristic is different. For example, the parameters yielding the worst t _{ r } may not yield the worst t _{ s }, M _{ OV } and e _{ ss } and vice versa. For illustration, the step responses depicting the nominal and worstcase responses for each step response characteristics for Δ∈{0.2,0.5,1.0,1.2} are shown in Fig. 6 for both PI and CMC controllers.
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.
The step response characteristics for both the nominal systems are tabulated in Table 2. As previously, the step responses depicting the nominal and worstcase responses for Δ∈{0.2,0.5,1.0} are shown in Fig. 7 for the PI and CMC controllers respectively. Note that we do not consider the case for Δ=1.2 as the closedloop system becomes unstable for Δ=1.0, when the PI controller + nonlinear process is used.
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
The results thus far have shown consistently better robustness from the CMC controller compared to the PI controller. To explain this, we analyse the mapping of steadystate inputoutput signals of these two controllers. Fig. 8(a) shows the mapping of steadystate inputoutput signals of both controllers as they were implemented when controlling the linear process. The mapping of inputoutput signals for the nominal system and the maximum deviation from this response when Δ=1.2 are shown in black solid line and magenta dashdotted line respectively. We observe a significantly greater change to the gradient of the PI controller’s inputoutput mappings compared to the CMC controller.
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.
Given the process to be controlled is a linear process, its ODE representation (with x:=y) is given by
Equation 12 is in the standard statespace representation (i.e. \(\frac {dx}{dt} = Ax + Bu, y = Cx + Du\)) with A=−k _{ p2} and B=k _{ p1}, C=1 and D=0. In linear control theory design using a statespace approach, [42], a standard control law can be written as u=Kx where K is the controller gain. This linear control law can be viewed as a mapping of input, x to output, u with K being the gradient. Substituting u=Kx into Eq. (12), we have
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.
What makes the CMC more robust (in terms of gradient change) to parameter uncertainty? Our simulation results using the nonlinear process shed some light on this matter. The steadystate mapping of inputoutput signals simulated 1060 times at Δ=1.0 for both the controllers when controlling the nonlinear process is shown in Fig. 9(a). For the nominal system both the controllers’ inputoutput mapping retains a the linear behaviour. While the PI controller’s steadystate inputoutput mapping stays linear for all 1060 uncertainty combinations, the CMC controller’s inputoutput mapping displays a ‘hyperbolic’ behaviour for some parameter combinations. Recall that this ‘hyperbolic’ behaviour is one of the inputoutput signal mappings reported in [25] (see also Fig. 5(b)). Thus, our simulation results seem to indicate that parameter uncertainty has the capacity to change the operating regime of the CMC controller from signaltransducing to hyperbolic. Thus, the question we are interested in is whether this change in the mapping regime accounts for the better robustness of the CMC controller.
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
References
 1
Purnick PEM, Weiss R. The second wave of synthetic biology: from modules to systems. Nat Rev Mol Cell Biol. 2009; 10:410–22.
 2
Kahl LJ, Endy D. A survey of enabling technologies in synthetic biology. J Biol Eng. 2013; 7(1):1–19.
 3
Liu Y, Chakrabartty S. Factor graphbased biomolecular circuit analysis for designing forward error correcting biosensors. IEEE Trans Biomed Circ Syst. 2009; 3:150–9.
 4
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.
 5
Sharma P, Roy S. Alloptical biomolecular parallel logic gates with bacteriorhodopsin. IEEE Trans NanoBioscience. 2004; 3(2):129–36.
 6
GoñiMoreno A, Amos M. A reconfigurable nand/nor genetic logic gate. BMC Syst Biol. 2012; 6(1):1–11.
 7
Marchisio MA. In silico design and in vivo implementation of yeast gene boolean gates. J Biol Eng. 2014; 8(1):1–15.
 8
Chen BS, Wu CH. A systematic design method for robust synthetic biology to satisfy design specifications. BMC Syst Biol. 2009; 3(1):1–18.
 9
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.
 10
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.
 11
Sarpeshkar R. Analog versus digital: extrapolating from electronics to neurobiology. Neural Computation. 1998; 10(7):1601–38.
 12
Sauro HM, Kim K. Synthetic biology: it’s an analog world. Nature. 2013; 497(7451):572–3.
 13
Daniel R, Rubens JR, Sarpeshkar R, Lu TK. Synthetic analog computation in living cells. Nature. 2013; 497(7451):619–24.
 14
Xiong L, Schumaker KS, Zhu JK. Cell signaling during cold, drought and salt stress. Plant Cell. 2002; 14:165–83.
 15
Poli G, Leonarduzzi G, Biasi F, Chiarpotto E. Oxidative stress and cell signalling. Curr Med Chem. 2004; 11(9):1163–82.
 16
Gomperts BD, Kramer IM, Tatham PER. Signal transduction. Massachusetts: Academic Press; 2004.
 17
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.
 18
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.
 19
Krebs EG, Beavo JA. Phosphorylationdephosphorylation of enzymes. Annu Rev Biochem. 1979; 48:923–59.
 20
Miller CA, Sweatt JD. Covalent modification of DNA regulates memory formation. Neuron. 2007; 53:857–69.
 21
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.
 22
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.
 23
Ferrell Jr JE, Ha SH. Ultrasensitivity part I: Michaelian responses and zeroorder ultrasensitivity. Trends Biochem Sci. 2014; 39:496–503.
 24
Heinrich R, Neel BG, Rapoport TA. Mathematical models of protein kinase signal transduction. Molecular Cell. 2002; 9:957–60.
 25
GomezUribe C, Verghese GC, Mirny LA. Operating regimes of signaling cycles: statics, dynamics, and noise filter. PLoS Comput Biol. 2007; 3(12):246.
 26
Oishi K, Klavins E. Biomolecular implementation of linear I/O systems. IET Syst Biol. 2011; 5:252–60.
 27
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.
 28
Buisman HJ, ten Eikelder HMM, Hilbers PAJ, Liekens AML. Computing algebraic functions with biochemical reaction networks. Artif Life. 2009; 15:1–15.
 29
Seelig G, Soloveichik D, Zhang DY, Winfree E. Enzymefree nucleic acid logic circuits. Science. 2006; 314(5805):1585–8.
 30
Soloveichik D, Seelig G, Winfree E. DNA as a universal substrate for chemical kinetics. Proc Natl Acad Sci. USA. 2010; 107(12):5393–8.
 31
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.
 32
Horn F. Necessary and sufficient conditions for complex balancing in chemical kinetics. Arch Ration Mech Anal. 1972; 49:172–86.
 33
Klonowski W. Simplifying principles for chemical and enzyme reaction kinetics. Biophys Chem. 1983; 18:73–87.
 34
Qian L, Winfree E. A simple DNA gate motif for synthesizing largescale circuits. J. Royal Soc. Interface. 2011. doi:10.1098/rsif.2010.0729.
 35
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.
 36
Abramowitz M, Stegun IA. Handbook of mathematical functions. Appl Math Ser. 1966; 55:62.
 37
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.
 38
Vidyasagar M. Statistical learning theory and randomised algorithms for control. IEEE Control Syst. 1998; 18(6):69–85.
 39
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.
 40
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.
 41
Franklin GF, Powell JD, EmamiNaeini A. Feedback control of dynamic systems, 7th edn. New York: Pearson; 2015.
 42
Friedland B. Control system design: an introduction to statespace methods. New York: McGrawHill; 2005.
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.
Author information
Affiliations
Corresponding author
Rights and permissions
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.
About this article
Cite this article
Foo, M., Sawlekar, R. & Bates, D.G. Exploiting the dynamic properties of covalent modification cycle for the design of synthetic analog biomolecular circuitry. J Biol Eng 10, 15 (2016). https://doi.org/10.1186/s1303601600361
Received:
Accepted:
Published:
DOI: https://doi.org/10.1186/s1303601600361
Keywords
 Covalent modification cycle
 Chemical reaction networks
 Analog synthetic biomolecular circuits
 Linear and nonlinear operators
 Feedback control systems
 Synthetic biology applications