PID Control of Biochemical Reaction Networks
Max Whitby, Luca Cardelli, Marta Kwiatkowska, Luca Laurenti, Mirco, Tribastone, Max Tschaikowski

TL;DR
This paper introduces a novel derivative component in biochemical reaction networks, enabling PID control in synthetic biological systems, demonstrated by regulating gene expression in a microRNA model.
Contribution
First implementation of a derivative component in CRNs, facilitating PID control for biochemical systems.
Findings
Successfully implemented a derivative component in CRNs.
Demonstrated PID control in gene expression regulation.
Showed potential for advanced control in synthetic biology.
Abstract
Principles of feedback control have been shown to naturally arise in biological systems and successfully applied to build synthetic circuits. In this work we consider Biochemical Reaction Networks (CRNs) as a paradigm for modelling biochemical systems and provide the first implementation of a derivative component in CRNs. That is, given an input signal represented by the concentration level of some species, we build a CRN that produces as output the concentration of two species whose difference is the derivative of the input signal. By relying on this component, we present a CRN implementation of a feedback control loop with Proportional-Integral-Derivative (PID) controller and apply the resulting control architecture to regulate the protein expression in a microRNA regulated gene expression model.
Peer Reviews
No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.
Videos
No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.
PID Control of Biochemical Reaction Networks
Max Whitby1, Luca Cardelli1, Marta Kwiatkowska1, Luca Laurenti1, Mirco Tribastone2, Max Tschaikowski3 1 Department of Computer Science, University of Oxford, UK2 IMT School for Advanced Studies, Lucca, Italy3 Department of Computer Engineering, TU Wien, Austria
Abstract
Principles of feedback control have been shown to naturally arise in biological systems and successfully applied to build synthetic circuits. In this work we consider Biochemical Reaction Networks (CRNs) as a paradigm for modelling biochemical systems and provide the first implementation of a derivative component in CRNs. That is, given an input signal represented by the concentration level of some species, we build a CRN that produces as output the concentration of two species whose difference is the derivative of the input signal. By relying on this component, we present a CRN implementation of a feedback control loop with Proportional-Integral-Derivative (PID) controller and apply the resulting control architecture to regulate the protein expression in a microRNA regulated gene expression model.
I INTRODUCTION
Biochemical Reaction Networks (CRNs) are a widely used formalism to describe biochemical systems [1]. More recently, they have also been employed as a formal programming language for synthetic circuits made of DNA [2, 3]. Due to the numerous potential applications, ranging from smart therapeutics to biosensors, the construction of CRNs that exhibit prescribed dynamics is a major goal of synthetic biology. However, achieving a desired behaviour by designing a CRN is difficult due to the complexity of such systems and limited knowledge of their dynamics [4, 5].
Negative feedback and Proportional-Integral-Derivative (PID) control are widely used in engineering to control the dynamics of a system due to their ability to achieve accurate set-point tracking and robustness to disturbances, even with only partial knowledge of the system. Because of these properties, such mechanisms have also been applied with success in the construction of synthetic bio-molecular systems [6, 7]. Moreover, molecular implementation of control systems has been shown to naturally occur in living organisms [8, 9, 10]. For example, integral control occurs in E.coli chemotaxis [11, 12], while CheY proteins regulate the bacteria’s tumbling frequency by implementing a derivative control [13]. As a consequence, in view of the potential applications, CRN designs that implement control mechanisms are sought for [12, 6]. CRNs implementing proportional and integral control have been proposed [14, 12]. However, a biochemical implementation of a full PID control is still missing due to the lack of a CRN implementing the derivative component.
In this work we first present a CRN implementation of a derivative component. That is, we provide a CRN such that, given an input signal represented by the concentration level of some species, the output is the concentration of two species whose difference gives the derivative of the input signal. We use this as a building block for a PID controller, and show how negative feedback with PID controller can be implemented in CRNs. We show the effectiveness of this architecture on a microRNA regulated gene expression example [15, 16], where we control the time evolution of a protein by acting on the expression of mRNA and microRNA.
In summary, we make the following contributions:
- •
We present a CRN that computes the derivative of an input signal and prove its asymptotic correctness.
- •
We extend the correctness results of CRN encodings [14] of proportional and integral signals to the case of nonlinear dynamics.
- •
We provide for an arbitrary CRN plant a CRN encoding of the PID feedback controller.
- •
We show the effectiveness of our control architecture on a microRNA regulated gene expression model.
II Biochemical Reaction Networks
In this section we provide some background about the deterministic mass-action semantics of a CRN based on the reaction-rate equations. Then we review the notion of dual rail encoding for a species. Finally we fix a graphical representation of CRNs that will be used throughout the paper.
II-A Deterministic Mass-action semantics
A CRN is a pair of finite sets, where is an ordered set of species, denotes its size, and is an ordered set of reactions. Species in interact according to the reactions in . A reaction is a triple , where is the reactant complex, is the product complex and is the coefficient associated with the rate of the reaction. Complexes and represent the stoichiometry of reactants and products. We denote the -th component of complex by ; the zero complex is denoted by . Given a CRN with species set , a reaction will be denoted by . The state change associated to is defined by . For example, the state change of the reaction above is .
We consider the deterministic interpretation of a CRN based on the well-known reaction-rate equations with mass-action kinetics. Given a CRN and an initial condition representing the initial concentration of each species, the time course of the concentrations can be described as the solution of an initial value problem with the following system of ODEs
[TABLE]
and initial condition . For a species we denote by the concentration of at time .
In this paper we synthesise PID controllers with mass-action kinetics for CRNs. We will also assume that the plant is represented by a mass-action CRN, although our results carry over to plants given in terms of smooth control systems.
II-B Dual Rail Encoding
The plant is a CRN, hence its output is given by non-negative solutions. However, the PID controller will involve quantities that are negative such as the error, i.e., difference between the set-point and the output, as well as its derivative. In order to handle this we will use the so-called dual rail encoding [14], by which a a signal is decomposed into a “positive” and “negative” species component whilst preserving a law of mass action kinetics such that each individual species concentrations cannot be negative. Specifically, for a signal we denote the two distinct component species by and , representing the positive and negative signals, respectively. Owning to the fact that the plant is described by a mass-action CRN, we wish to point out that its non-negative output can (and is) captured by single rail encoding.
II-C Graphical Representation of CRNs
A CRN can be represented as a labelled directed bipartite graph according to the usual Petri net representation with species and reaction nodes [17]. A reaction node is labelled with a rate coefficient. There is an edge from a species node to a reaction node if the species is in the reactant complex, with label equal to its multiplicity; similarly, an edge from a reaction node to a species node indicates the presence of a species in the product complex. For example, we have the following representation for the reaction :
Throughout the rest of the paper, for ease of presentation we do not draw labels on edges if the related complex multiplicity is 1. We also remove the black box representing reaction nodes to reduce clutter. Finally, we introduce a short-hand notation for recurring reaction patterns as shown in Figure 1, where each arc is either a pointed arrow () or a rounded arrow (), with the source represented by the flat edge and the target represented by the arrow head. A pointed arrow represents a non-catalytic reaction and a rounded arrow represents a catalytic reaction.
III CRN implementation of the PID controller
We introduce the CRN implementation of the proportional, integral, and derivative components of a PID controller. We describe them as blocks where the input species are (which will indicate the dual-rail error signal between the species representing the set-point and the plant output). The output of the PID controller is denoted by . The proportional and integral components have been already introduced for linear control systems [14]. Here we prove their correctness in the presence of non-linearity. In addition we detail the CRN implementation of the derivative block, which is a novel contribution to the best of our knowledge.
As shown in Figure 2, as in a classic feedback loop, the signals synthesised by the PID controller act on the CRN plant, whose output is measured, and sent back as input of the PID controller after comparison with the reference signal. The output of the plant is always given by a species . The objective of the control is to have follow the reference signal.
To this end, we let denote the mass-action CRN representing the plant and construct a CRN encoding of a PID feedback law as indicated in Figure 2. Given that output and control signals are vectors in general, the blocks of Figure 2 are multidimensional components in general. In particular, assuming that denotes the dimension of the output and control vector, the CRN encoding of the PID feedback law is given by interconnected
- •
subtraction blocks
- •
addition blocks
- •
proportional blocks
- •
integral blocks
- •
derivative blocks
- •
dual rail converter blocks .111Please note that imposing the presence of a , and block for every coordinate is without loss of generality because a block can be removed by setting its multiplier to zero.
With this, the overall CRN is given by
[TABLE]
where the feedback law CRN is defined by
[TABLE]
In what follows next, we address the correctness of each block type.
III-A Proportional, Addition, Subtraction and Dual Rail Converter Blocks
We begin by presenting the proportional block which computes an output signal that is proportional to the input signal.
Definition 1** (Proportional block [14]).**
For input species , output species , parameters , and the multiplier , the proportional block is a CRN composed by the following reactions
[TABLE]
Theorem 1**.**
On any bounded time interval, the ODE system of (2) converges to an ODE system satisfying and for all if in all proportional blocks.
Proof (Sketch).
Note that
[TABLE]
This motivates to interpret all and as fast variables in the sense of Tikhonov’s theorem [18, Section 8.2]. To see that the requirements of the theorem are satisfied, we note that the fast ODE system (i.e., the one consisting of fast variables) admits, for any fixed vector of slow variables (i.e., all variables that are not fast), exactly one possible equilibrium point. Moreover, it is an asymptotically stable equilibrium of the fast ODE system. We finish the proof by noting that smooth exogenous reference signals can be captured because Tikhonov’s theorem applies to non-autonomous smoooth ODE systems. ∎
Remark 1**.**
Theorem 1 extends the result of [14] to nonlinear control systems. The same holds true for the other blocks of this section.
Remark 2**.**
The proof of Theorem 1 reveals that reaction is not strictly needed to ensure correctness. However, this reaction precludes and from attaining excessively large values (recall that is large), thus reducing the impact of numerical errors.
The correctness of proportional, subtraction and converter blocks introduced next and depicted in Figure 2 is shown similarly to Theorem 1.
More specifically, the addition block is given by.
Definition 2** (Addition block [14]).**
For input species and , output species , and parameters the addition block is a CRN composed by the following reactions
[TABLE]
The subtraction block, instead, is defined as follows.
Definition 3** (Subtraction block [14]).**
For input species and , output species , and parameters the subtraction block is a CRN composed by the following reactions
[TABLE]
At last, the converter block is described by the following.
Definition 4** (Dual rail converter block).**
For input species , output species , and parameters the single to dual rail converter block is a CRN composed by the following reactions
[TABLE]
III-B Integral Block
The integral component computes a multiple of the integral of the input signal. This action is widely used in control systems due to its ability to collect past information about the error to be corrected. A CRN implementation of the integral component has been proposed in [14] and is reported in Definition 5.
Definition 5** (Integral block [14]).**
For input species , output species , some constant and multiplier , the integral block is given by the following CRN
[TABLE]
Proposition 1**.**
The integral blocks introduced in Definition 5 are correct. More formally, for all , it holds that
[TABLE]
Proof.
Straightforward via differentiation (the values at are chosen appropriately). ∎
Proposition 1 can be seen as a (straightforward) extension of the corresponding result in [14] to nonlinear CRN plants.
III-C Derivative Block
The derivative block computes a multiple of the derivative of the input signal. This component is used in control systems to predict the future error given its current trend, and thus to help dampen oscillations introduced by P and I components.
Building a derivative module by chemical reactions is challenging because on-the-fly differentiation can only be done by comparing a signal at two time points, inherently requiring an approximation dependent on the time difference. This is resolved by the circuit in Figure 2(d), which handles dual rail input and output. Intuitively, the inputs and are sampled at two time points, , and , , respectively, and a multiple of their difference is provided via , . In Figure 2(d), the two reactions and cause to track with a (slight) delay dependent on . Intuitively, this yields . The symmetric two reactions similarly cause to track , which ensures that . The three reactions , , and cause to track with delay dependent on . The symmetric three reactions similarly cause to track . Thus tracks with delay dependent on . This and the above discussion allow us then to conclude that .
The next theorem formalises the above considerations using Tikhonov’s theorem.
Definition 6**.**
For input species , auxiliary species , output species , parameters and the multiplier , the derivative block is a CRN composed by the following reactions
[TABLE]
The following theorem shows that the above CRN is such that, under certain scaling of the rates, produces a correct approximation of the derivative of .
Theorem 2**.**
The derivative blocks are asymptotically correct. In particular, the following holds.
The solution of (2) converges, on any bounded time interval, to an ODE system which satisfies , for all , if in all derivative blocks. 2. 2.
The solution of (2) converges, on any bounded time interval, to an ODE system which satisfies and for all , if in all derivative blocks.
Proof (Sketch).
To see , we first note that Definition 6 yields
[TABLE]
Since this implies
[TABLE]
this motivates us to add to the ODE system of (2) the additional ODE
[TABLE]
and to replace each instance of in the ODE system with . By processing the other derivative blocks in a similar fashion, we introduce new ODE variables . To see that the requirements of Tikhonov’s theorem are satisfied, we note that the fast ODE system (i.e., the one containing , and ) admits, for any fixed vector of slow variables, exactly one equilibrium point. Additionally, it is an asymptotically stable equilibrium of the fast ODE system. To see , instead, we apply Tikhonov’s theorem in the case where, in every derivative block, and are treated as fast variables (while all other variables are considered to be slow) and in all derivative blocks. ∎
IV PID control of gene expression
In this section we apply the PID feedback control architecture developed in this paper to a gene expression model. In particular, we consider a regulated gene expression model from [16], for which synthetic implementations have already been proposed in [19]. The model is composed by the following reactions, where for simplicity we fixed unitary kinetic parameters
[TABLE]
That is, we have that catalyses the production of the protein and is down-regulated by an annihilation reaction with the .
The objective of our control is to have the protein to follow a reference signal. Given and , the control signals synthesised by the controller, we assume that these can act on the plant by regulating the expression rate of and , respectively. This assumption is justified by the fact that these mechanisms can be implemented synthetically [19]. We consider the following reactions to model such actuation:
[TABLE]
In this model, a high concentration of will increase the production rate of and so of , whereas a high concentration of will decrease the amount of by producing with a higher rate.
In the actuation model considered above we have that the control signals act on two different species. This is not a requirement of our architecture. Another possible actuation is that annihilates directly. This can be modelled with the following reactions
[TABLE]
Finally, another possibility is that and acts directly on the target species . In this case, the actuation is
[TABLE]
In Figure 3 we consider the different actuation mechanisms described above and compare the performance of PI and PID controllers for two different reference signals: a constant signal and an oscillatory signal. For all the plots in the figure we considered the same parameters for PI and PID controllers (reported in the Appendix). It is easy to observe that, whereas a negative feedback with PI controller can already track both signals correctly, in the case of a PID controller the time evolution of the concentration of has reduced oscillations around the reference signals. This is due to the action of the derivative block. In fact, while it is well known that the derivative component in a PID does not necessarily reach zero error at steady state, it can help to reduce the transient error between the output and the reference signals and to dampen oscillations around the set points.
V Conclusion
In this work we considered feedback control with PID controllers and proposed a CRN implementation for this control architecture. This relies on a novel CRN, which computes the derivative of an input molecular signal. We applied our framework to control the protein expression in a regulated gene expression model and showed improved performance compared to a PI feedback control. An interesting aspect, which has not been considered in this paper, is to study the effect that the proposed control system has on noise [20]. This is left as future work.
-A CRN for PID control of gene expression
We present our full Chemical Reaction Network PID feedback loop with gene expression plant and two reference signals introduced in Section IV. We include the three actuation mechanisms given with the plant. For each block we also give the initial conditions used to produce the simulations seen in Figure 3.
PID Controller
First we report the reactions and parameters of the CRN PID controller of which the proof of correctness and details of operation are outlined in Section 3. For all figures we used the same parameters
Proportional Block:
[TABLE]
Initial Conditions: (rates) (species) ,
Integral Block:
[TABLE]
Initial Conditions: (rates) , (species) ,
Derivative Block:
[TABLE]
Initial Conditions: (rates) , , (species) ,,,
-B Summation and Subtraction Blocks
We provide the CRNs for the two summation blocks which are highlighted in the green block in Figure 2 and noted upon in section 3. The first adds the output of the P and I blocks together. The second the PI and D blocks together which produces the output species of the controller (given as in the CRNs below).
Addition Block :
[TABLE]
Initial Conditions for summation block: (rates) ,, (species) ,
Addition Block :
[TABLE]
Initial Paramterisation for summation block: (rates) , (species) ,
Subtraction Block:
The subtraction block is used to compute the error of the output of the plant with the reference signal. It is detailed further in Section 3:
[TABLE]
Initial Conditions for difference block summation block: (rates)
-C Reference Signals
Next we introduce the intitial conditions and reactions for both the constant and sine wave reference signals outlined further in Section IV. These act as a reference which we are trying to control our plant to track.
Constant:
The constant signal can be given simply by stating a non-decaying species with a molecular count equal to the constant signal however it can also be given by the following CRN which is stated in Figure 3.
[TABLE]
Initial Conditions: (rates) (species) ,
Sine wave: The sine wave has a slow reaction rate to allow for the PID controller to properly track the signal.
[TABLE]
Initial Conditions: (rates) (species) , , ,
Plant and Actuators
We introduce the gene expression plant used within our model. We also include the three actuations methods from the controller seen in Figure 3 and discussed in section 4 used to interface with the plant.
Actuator 1:
[TABLE]
Initial Conditions: (rates) , , (species) ,
Actuator 2:
[TABLE]
Initial Conditions: (rates) , (species) ,
Actuator 3:
[TABLE]
Initial Conditions: (rates) , (species) ,
**Plant: **
[TABLE]
Initial Conditions: (rates) , (species) , ,
-D Single to Dual Rail Block
We introduce the block which takes the output of the plant and transforms into into a dual rail signal, see blue block Figure 2 and discussed in section 3.
[TABLE]
Initial Conditions: (rates) , (species) ,
The reference list from the paper itself. Each links out to its DOI / PubMed record.
- 1[1] Péter Érdi and János Tóth. Mathematical models of chemical reactions: theory and applications of deterministic and stochastic models . Manchester University Press, 1989.
- 2[2] David Soloveichik, Georg Seelig, and Erik Winfree. DNA as a universal substrate for chemical kinetics. Proceedings of the National Academy of Sciences , 107(12):5393–5398, 2010.
- 3[3] Yuan-Jyue Chen, Neil Dalchau, Niranjan Srinivas, Andrew Phillips, Luca Cardelli, David Soloveichik, and Georg Seelig. Programmable chemical controllers made from DNA. Nature nanotechnology , 8(10):755, 2013.
- 4[4] Luca Cardelli, Milan Češka, Martin Fränzle, Marta Kwiatkowska, Luca Laurenti, Nicola Paoletti, and Max Whitby. Syntax-guided optimal synthesis for chemical reaction networks. In International Conference on Computer Aided Verification , pages 375–395. Springer, 2017.
- 5[5] Tai-Yin Chiu, Hui-Ju K Chiang, Ruei-Yang Huang, Jie-Hong R Jiang, and François Fages. Synthesizing configurable biochemical implementation of linear systems from their transfer function specifications. Plo S one , 10(9):e 0137442, 2015.
- 6[6] Domitilla Del Vecchio, Aaron J Dy, and Yili Qian. Control theory meets synthetic biology. Journal of The Royal Society Interface , 13(120):20160380, 2016.
- 7[7] Ciarán L Kelly, Andreas W K Harris, Harrison Steel, Edward J Hancock, John T Heap, and Antonis Papachristodoulou. Synthetic negative feedback circuits using engineered small rnas. Nucleic acids research , 46(18):9875–9889, 2018.
- 8[8] Fuzhong Zhang, James M Carothers, and Jay D Keasling. Design of a dynamic sensor-regulator system for production of chemicals and fuels derived from fatty acids. Nature biotechnology , 30(4):354, 2012.
