Signal differentiation with genetic networks 1

Proceedings of the 20th World Congress Proceedings of the 20th World The International Federation of Congress Automatic Control Proceedings of the the 20th World World Congress Proceedings of 20th The International of Congress Automatic Control Toulouse, France,Federation July 9-14, 2017 Available online at www.sciencedirect.com The International Federation of The International of Automatic Automatic Control Control Toulouse, France,Federation July 9-14, 2017 Toulouse, France, July 9-14, 2017 Toulouse, France, July 9-14, 2017

ScienceDirect

IFAC PapersOnLine 50-1 (2017) 10938–10943  Signal differentiation with genetic networks  Signal differentiation with genetic networks  Signal with networks Signal differentiation differentiation with genetic genetic networks ∗ ∗ ∗

Wolfgang Halter ∗ , Zoltan Wolfgang Halter ∗∗ ,, Zoltan Wolfgang Wolfgang Halter Halter , Zoltan Zoltan ∗ Institute for Systems Theory and ∗ ∗ Institute for Systems Theory and ∗ Institute for Theory Germany of corresponding Institute (e-mail for Systems Systems Theory and and Germany Germany (e-mail (e-mail of of corresponding corresponding Germany (e-mail of corresponding

A. Tuza ∗ , Frank Allg¨ ower ∗ A. Tuza o wer ∗∗ ∗ , Frank Allg¨ ∗ , Frank Allg¨ A. Tuza o A. Tuza , Frank Allg¨ ower wer Automatic Control, Stuttgart, 70569 Automatic Control, Stuttgart, 70569 Automatic Control, author: [email protected]). Automatic Control, Stuttgart, Stuttgart, 70569 70569 author: [email protected]). author: [email protected]). [email protected]). author:

Abstract: Biological counterparts of differential operators are not only likely to exist in natural Abstract: counterparts of differential are not only Therefore, likely to exist natural Abstract: Biological counterparts of operators are not likely in natural systems butBiological also desirable to realize as modulesoperators in synthetic we in present Abstract: Biological counterparts of differential differential operators are biology. not only only Therefore, likely to to exist exist in naturalaa systems but also desirable to realize as modules in synthetic biology. we present systems but also desirable to realize as modules in synthetic biology. Therefore, we present genetic regulatory network which approximates the differential operator. For the in silico design,aa systems but also desirable to realize as modules in synthetic biology. Therefore, we present genetic regulatory network which approximates the differential operator. For the in silico design, genetic regulatory network which approximates the differential operator. For the in silico design, we use a modeling framework which takes into account the limitations imposed by finite pools of genetic regulatory network which approximates the differential operator. For theby infinite silicopools design, we use a modeling framework which takes into account the limitations imposed of we use a modeling framework which takes into account the limitations imposed by finite pools of resources; we synthesize the differentiator module by combining genetic regulatory parts which we use a modeling framework which takes into account the limitations imposed by finite pools of resources; we synthesize the differentiator module by combining genetic regulatory parts which resources; we synthesize the differentiator module by combining genetic regulatory parts which realize basic input/output functions such as a gain, integrator and signal difference. The resulting resources; we synthesize the differentiator module by combining genetic regulatory parts which realize basic input/output functions such as a gain, integrator and signal difference. The resulting realize basic functions as integrator and difference. The resulting two-gene-network approximates thesuch transfer function of a differential with realize basic input/output input/output functions such as aa gain, gain, integrator and signal signaloperator difference. Theadditional resulting two-gene-network approximates the transfer function of aawell differential operator with additional two-gene-network approximates the transfer function of differential operator with additional low-pass filter. The functionality of this small network, as as its robustness towards changes two-gene-network approximates the transfer function of a differential operator with additional low-pass filter. The functionality of this small network, as well as its robustness towards changes low-pass filter. The functionality of this small network, as well as its robustness towards in the cellular environment, is investigated numerically. low-pass filter. environment, The functionality of this small network, as well as its robustness towards changes changes in the cellular is investigated numerically. in the cellular environment, is investigated numerically. in the cellular environment, is investigated numerically. © 2017, IFAC (International Federation of Automatic Control) Hosting by Elsevier Ltd. All rights reserved. Keywords: Bio Control, Differentiators, Integrators, Transfer functions, Robust performance Keywords: Bio Bio Control, Differentiators, Differentiators, Integrators, Transfer Transfer functions, Robust Robust performance Keywords: Keywords: Bio Control, Control, Differentiators, Integrators, Integrators, Transfer functions, functions, Robust performance performance 1. INTRODUCTION 1. INTRODUCTION INTRODUCTION y u 1. K 1. INTRODUCTION y u − K y u K y − The main goal of synthetic biology is to realize predefined u K − The main main in goal of synthetic synthetic biology is to to realize realize predefined − The goal of biology is predefined functions biological systems. Control engineers provide The main in goal of synthetic biology is to realize predefined functions biological systems. Control engineers provide functions in systems. Control engineers system theory based tools and methods rationalprovide design functions in biological biological systems. Control for engineers provide system theory based tools and methods for rational design Go (s) system based tools and methods for rational design of such theory functions. For example, based on biochemical reG system theory based tools and methods for rational design o (s) of such functions. For example, based on biochemical reG Goo (s) (s) of such functions. For example, based on biochemical reaction networks limited to first order reactions, Oishi and of such functions. For example, based on biochemical reaction networks networks limited to toafirst first order reactions, reactions, Oishi and action limited order Oishi and Klavins (2011) introduced framework to realize arbitrary action networks limited toafirst order reactions, and Fig. 1. Approximation of the inverse of Go (s). K needs to Klavins (2011)functions. introduced framework toframework realizeOishi arbitrary Klavins (2011) introduced aa framework realize input/output Therefore, theirto estab- Fig. 1. Approximation of the inverse of G needs to o (s). K Klavins (2011)functions. introduced framework toframework realize arbitrary arbitrary Fig. be 1. of chosen large enough. input/output Therefore, their estab1. Approximation Approximation of the the inverse inverse of of G Goo (s). (s). K K needs needs to to input/output functions. Therefore, their framework establishes a direct link between linear control circuit design and Fig. be chosen large enough. input/output functions. Therefore, their framework estabbe chosen chosen large large enough. enough. lishes a direct link between linear control circuit design and be theoretic mechanisms leading to this behavior is a crucial lishes a direct link between linear control circuit design and biochemical reaction networks. Building on this approach lishes a directreaction link between linearBuilding control circuit and theoretic mechanisms leading to this behavior is a crucial biochemical networks. on based thisdesign approach mechanisms leading to is step towards the biological realization of a differential opbiochemical Building on this Yordanov et reaction al. (2014)networks. implemented a DNA integral theoretic theoretic mechanisms leading to this this behavior behavior is a a crucial crucial biochemical reaction networks. Building on based this approach approach step towards the biological realization of aa differential opYordanov et al. (2014) implemented a DNA integral step towards the biological realization of differential operator. Despite the importance of the differential operator Yordanov et al. (2014) implemented a DNA based integral feedback controller which achieves perfect adaptation of an step towards the biological realization of a differential opYordanov et al. (2014) implemented a DNA based integral erator. Despite the importance of the differential operator feedback controller which achieves perfect adaptation of an erator. Despite the importance of the differential operator in a biological context, this problem has been investigated feedback controller which achieves perfect adaptation of an output signal to certain inputs. While Briat et al. (2016) erator. Despite the importance of the differential operator feedback controller which achieves perfect adaptation of an in aa biological context, this problem has been investigated output signal signal togeneric certainsolution inputs. to While Briat et al. (2016) (2016) in context, this problem has been investigated sporadically (Harris al., 2015; Sontag, output certain inputs. While al. offered a moreto the Briat same et model class, rather in a biological biological context, this et problem has Lang been and investigated output signal togeneric certainsolution inputs. to While Briat et al. (2016) rather sporadically (Harris et al., 2015; Lang and Sontag, offered a more the same model class, rather sporadically (Harris et al., 2015; Lang and Sontag, 2016). offered a more generic solution to the same model class, Ang et al. (2010) provide integral feedback control for rather sporadically (Harris et al., 2015; Lang and Sontag, offered aal. more generic solution to the same model class, 2016). Ang et (2010) provide integral feedback control for 2016). We thus adapt the approach of Oishi and Klavins (2011) in Ang et al. (2010) provide integral feedback control for genetic regulatory networks using transcriptional control 2016). Ang et al. (2010) provide integral feedback control for We thus adapt the approach of Oishi and Klavins (2011) in genetic of regulatory networks using using transcriptional control We thus adapt the approach of Oishi and Klavins (2011) order to derive a genetic motif which realizes a differential genetic regulatory networks transcriptional control instead DNA interactions. The interactions in such sysWe thus adapt the approach of Oishi and Klavins (2011) in in genetic regulatory networks using transcriptional control order to derive a genetic motif which realizes a differential instead of DNA interactions. The interactions in such sysorder to to derive derive a genetic genetic motif which PID realizes a differential differential operator and may be used to realize control or other instead DNA interactions. The such systems areof usually modeled using Hillinteractions kinetics — in a nonlinear order a motif which realizes a instead of DNA interactions. The interactions in such sysoperator and may be used to realize PID control or other tems are are usually usually modeled using Hill kinetics —applicability nonlinear gradient and used PID or dependent in a biological system. We tems using kinetics aa phenomenological model — thusHill limiting the— operator and may may be befunctions used to to realize realize PID control control or other other tems are usually modeled modeled using Hill kinetics a nonlinear nonlinear operator gradient dependent functions in aa biological system. We phenomenological model — — thus limiting the—applicability applicability gradient dependent functions in biological system. first review the system theoretic idea how to approximate phenomenological model thus limiting the of linear design methods. One way to approach this probgradient dependent functions in idea a biological system. We We phenomenological model — thus limiting the applicability first review the system theoretic how to approximate of linear design methods. One way to approach this probreview the system idea to approximate afirst differential operator intheoretic frequency andhow time domain. Subof linear design methods. to approach problem is presented by Harris One et al.way (2015) where theythis linearize first review the system theoretic idea how to approximate of linear design methods. One way to approach this proba differential operator in frequency and time domain. Sublem dynamics is presented presented by Harris Harris et al. (2015) (2015) where they linearize a differential differential operator in frequency frequency andoftime time domain. Subweoperator discuss the adaptations the domain. framework of lem is by al. linearize the around the et steady statewhere of thethey system and asequently, in and Sublem dynamics is presented by Harris et al. (2015) where they linearize sequently, we discuss the adaptations of the framework of the around the steady state of the system and sequently, we discuss the adaptations of the framework of and Klavins (2011) before we provide an implementhe dynamics the steady state of the system then switch to around the frequency domain representation forand the Oishi sequently, we discuss the adaptations of the framework of the dynamics around the steady state of the system and Oishi and Klavins (2011) before we provide an implementhen switch to the frequency domain representation for the Oishi and Klavins (2011) before we provide an implementation of an approximating differentiator with a nonlinear then switch to the frequency domain representation for the controller design. Further examples and a comprehensive Oishi and Klavins (2011) before we provide an implementhen switch to the frequency domain representation for the tation of an approximating differentiator with aa nonlinear controller design. Further examples and comprehensive tation of an controller Further examples aa review of design. the current challenges of and control in synthetic genetic tation ofnetwork. an approximating approximating differentiator differentiator with with a nonlinear nonlinear controller design. Further examples and a comprehensive comprehensive genetic network. review of the current challenges of control in synthetic genetic network. review of the current challenges of control in synthetic biology is given by Del Vecchio et al. (2016). genetic network. review of the current challenges of control in synthetic biology is given by Del Vecchio et al. (2016). 2. SYSTEM THEORETIC IDEA biology is by al. The biological of a et differential i.e. biology is given given realization by Del Del Vecchio Vecchio al. (2016). (2016).operator, 2. SYSTEM THEORETIC IDEA The biological biological realization ofsignals, a et differential operator, i.e. 2. The realization of a differential operator, i.e. determination of gradients of typically used in PID 2. SYSTEM SYSTEM THEORETIC THEORETIC IDEA IDEA The biologicalofrealization ofsignals, a differential operator, i.e. determination gradients of typically used in PID Since the exact realization of a signal is determination of gradients of signals, typically used in PID controllers to improve controller performance, remains elu- Since the exact realization of a signal differentiator determination of gradients of signals, typically used in PID differentiator is controllers to improve controller performance, remains eluSince the exact realization of a signal differentiator is not causal, only approximating designs are possible and controllers to controller performance, remains elusive. A potential application of such an operator in biologthe exact realization of adesigns signal are differentiator is controllers to improve improve controller performance, remains elu- Since not causal, only approximating possible and sive. A potential application of such an operator in biolognot causal, only approximating designs are possible and usually causality is achieved by adding a low-pass filter. sive. A potential application of such an operator in biological systems would be the detection ofanthe rate of in change of not causal, only approximating designs arelow-pass possiblefilter. and sive. A potential application of such operator biologusually causality is achieved by adding a systems would be be thecondition. detectionAof of the the rate rate of changefor of In usually causality is achieved achieved by adding adding low-pass filter. frequency domain, such a system can aabelow-pass expressed as ical systems would detection change of aical certain environmental example causality is by filter. ical systems would be the thecondition. detectionAofcommon the rate of of changefor of usually In frequency domain, such aa system be expressed as a certain environmental common example s can In frequency domain, such system can be expressed as athis certain environmental condition. A common example for is chemotaxis, where individual cells orient themselves In frequency domain, such a system can be expressed as s G (1) (s) = a certain environmental condition. A common example for d this is chemotaxis, where individual cells orient themselves s G (1) (s) = τ s + 1 s d (s) = this is chemotaxis, where individual cells orient themselves with respect to a concentration gradient of the food source G (1) this isrespect chemotaxis, where individual cells of orient themselves d (s) = τ s + 1 G (1) with to a concentration gradient the food source dτ a design τ s + 1 with the time constant parameter. An alternawith respect to a concentration gradient of the food source (Alon, 2007).toIta concentration is therefore widely accepted thatsource such τ s + 1 with respect gradient of the food with the time constant τ a design parameter. An alterna(Alon, 2007). It is therefore widely accepted that such with the constant design parameter. alternaapproach derive ττ(1)aa from systemsAn is depicted (Alon, 2007). therefore widely accepted such operators existIt in is nature and the identification of system with the time time to constant designbasic parameter. An alterna(Alon, 2007). therefore widely accepted that that such tive tive approach to derive (1) from basic systems is depicted operators existItin in isnature nature and the the identification of system system tive approach to derive (1) from basic systems is depicted (s) is apin Fig. 1, where the inverse of an arbitrary G operators exist and identification of o tive approach to derive (1) from basic systemsGis depicted operators exist in nature and the identification of system (s) is apin Fig. 1, where the inverse of an arbitrary o (s)  Supported by the research cluster BW2 (www.bwbiosyn.de) of the is in Fig. 1, where the inverse of an arbitrary G proximated. Therein, the transfer function from u to is o (s) is yyapapin Fig. 1, where the inverse of an arbitrary G  Supported by the research cluster BW2 (www.bwbiosyn.de) of the o u to proximated. Therein, the transfer function from is 2 (www.bwbiosyn.de)  proximated. Therein, the transfer function from u to y given as Ministry for Science, Research and Art Baden-W¨ u rttemberg. of the Supported by the research cluster BW  Supported by the research cluster BW2 (www.bwbiosyn.de) of the proximated. Therein, the transfer function from u to y is is given as Ministry for Science, Research and Art Baden-W¨ urttemberg. given Ministry u given as as Ministry for for Science, Science, Research Research and and Art Art Baden-W¨ Baden-W¨ urttemberg. rttemberg. Copyright © 2017, 2017 IFAC 11425 2405-8963 © IFAC (International Federation of Automatic Control) Hosting by Elsevier Ltd. All rights reserved. Copyright © 2017 IFAC 11425 Copyright © 2017 IFAC 11425 Peer review under responsibility of International Federation of Automatic Control. Copyright © 2017 IFAC 11425 10.1016/j.ifacol.2017.08.2463

Proceedings of the 20th IFAC World Congress Wolfgang Halter et al. / IFAC PapersOnLine 50-1 (2017) 10938–10943 Toulouse, France, July 9-14, 2017

DNA

f (t, λ, Φ)

mRNA

g(t, λ, η, Φ)

ν Ø

Rrnap

protein

λc

λ

δ

x1

10939

λc x2

λc x3

Rrnap + M

λc ...

λc ν

xn

→Ø M−

Ø Rrib

Fig. 2. Protein synthesis depending on transcription and translation initiation rates λ and η, decay rates ν and δ as well as environmental parameters Φ. K Y = G(s) = (2) U 1 + Go (s)K and for large K approximates the inverse of Go (s), i.e. 1 lim G(s) = . (3) K→∞ Go (s) If Go (s) is now chosen to be the integrator 1s , we obtain the overall transfer function Ks (4) G(s) = s+K which is equal to (1) with K = τ1 and therefore approximately realizes a differentiator. We note that the approximate behavior is only achieved if K is chosen large enough. From this interpretation and from Fig. 1, we can see that a genetic network implementation of a differentiator can be realized with three basic input/output functions: signal difference, gain and integrator. While the approach followed by Harris et al. (2015) is based on transforming the genetic network model into frequency domain to match the transfer functions of network and desired function, viz. (1), we rather follow the approach of Oishi and Klavins (2011) to synthesize the controller based on the basic input/output functions. A representation in frequency domain is therefore not necessary, it rather remains to define the genetic implementations of these basic functions. 3. GENETIC DIFFERENTIATOR We first introduce our modeling framework for gene regulatory networks and then our approach to realize the necessary input/output functions with genetic interactions before we synthesize the genetic differentiator. 3.1 Modeling protein production and GRNs In general, the interactions between several genes (DNA) through their gene products (proteins) create a network of interaction. The commonly used modeling frameworks for gene regulatory networks assume that the protein production rate is a direct function of other proteins. This simple picture sometimes fails to capture the complexity of the interactions as the dynamics of protein production depend on several environmental conditions (e.g. available amounts of enzymes). Further, they are subject to time delays caused by the multistage production mechanism which is depicted in Fig. 2. In principle, the DNA is first transcribed into several copies of mRNA with a time dependent rate f depending on the transcription initiation rate λ and several environmental parameters Φ ∈ RL . These parameters include among others the total amount of RNA polymerase (RNAP) and ribosomes, the

×G ηc ηc

ηc

η z1

z2

z3

Rrib + P

ηc ...

ηc zm

δ

→Ø P −

×M Fig. 3. Scheme of the model with xi and zi the densities of RNAP/ribosomes on DNA/mRNA templates. transcription and translation elongation rates and other host dependent variables. After transcription, the mRNA copies are processed again and the protein is produced at a time dependent rate g, again depending on λ and Φ and additionally the translation initiation rate η. Further, both mRNA and protein are subject to first order decay with rates ν and δ, respectively. It is important to notice that for both reading processes, namely transcription and translation, RNAP and ribosomes have to first bind to the templates (DNA and mRNA) with initiation rates λ and η before they move along the templates while producing mRNA and protein respectively. In Halter et al. (2017), we introduced a modeling framework partially based on the results of Reuveni et al. (2011) which captures the above process on a relative high level of detail while remaining computationally tractable. Fig. 3 depicts the scheme of this model, emphasizing the movement of both RNAP and ribosomes along the templates. Therein, Rrnap and Rrib stand for the amounts of currently available RNAP and ribosomes, xi and zi for the densities of RNAP and ribosomes at the i-th location on the DNA and mRNA templates and the variables G, M and P for the amounts of DNA, mRNA and protein, respectively. Finally, λc and ηc are the elongation rates of transcription and translation. For a detailed description of this modeling framework, the reader is referred to Halter et al. (2017). We note that the transcription and translation rates f and g depicted in Fig. 2 are clearly depending on environmental conditions and, given the chosen modeling framework, in particular found as f (t, λ, Φ) = Gλc · xn (t, λ, Φ) g(t, λ, η, Φ) = M (t, λ, Φ) · ηc · zm (t, λ, η, Φ) . For the purpose of modeling networks of genes, initiation and degradation rates may be chosen dependent on mRNA or protein levels of other genes in the network such that e.g. λj = λj (M, P) where the subscript j stands for the j-th gene in a network with q distinct genes and M = [M1 . . . Mq ] 



P = [P1 . . . Pq ] the vectors representing all gene products (mRNA and protein) of the network.

11426

Proceedings of the 20th IFAC World Congress 10940 Wolfgang Halter et al. / IFAC PapersOnLine 50-1 (2017) 10938–10943 Toulouse, France, July 9-14, 2017

3.2 Input/output functions As pointed out in Section 2 and in Fig. 1, it is essential to realize three basic input/output functions — signal difference, gain and integration — as molecular interactions to approximate a genetic signal differentiator. Additionally, as molecule numbers can only take positive values, we need a representation for both positive and negative signals. In the following, we proceed in a fashion very similar to the approach of Oishi and Klavins (2011). We will repeatedly assume that some part of the protein synthesis process operates in its linear regime, which is necessary to exactly realize the desired input/output functions. However, several types of nonlinearities do not severely interfere with the desired behavior as long as no saturation effects occur which will be shown by our robustness analysis at the end of the paper. Signal summation and difference For adding or subtracting two signals u1 and u2 , we assume that these signals are somehow represented by the molecular amount of a certain mRNA or protein which can influence the transcription or translation initiation rate of a gene. It is crucial that these interactions are mutually exclusive to realize linear combinations of the signals, e.g. λ(u1 + u2 ) = λ(u1 ) ± λ(u2 ) has to hold for a transcription initiation rate which realizes this function. This for instance excludes cooperative transcription factors which rather multiply their action such that the absence of either signal yields a zero output. Additionally, we assume that changes in the initiation rate directly (ideally linearly) cause changes in the mRNA and/or protein levels. This is not a trivial assumption and the operation point of the gene has to be chosen appropriately, ideally in the linear range. Gain For realizing the amplification of a signal u with a biological system, where the dynamics is mainly governed by production and degradation rates, the most important feature this system needs to incorporate is that it reaches its steady state in reasonable time. In other words, production and degradation terms need to be of similar order of magnitude, consequently the ratio of production and degradation determines the amplification factor. In context of genetic networks, these features might most easily be achieved with the process of transcription as the degradation rates of mRNAs are usually larger than the ones of proteins. According to Fig. 2 we get M˙ = f (t, λ(u), Φ) − νM

and choose ν large enough and f (t, λ(u), Φ) ideally in a linear range with respect to u, s.t. f (t, λ(u), Φ) ≈ αu. In that case, the steady state of M can be expressed as f (t, λ(u), Φ) α MSS = = u ν ν resembling the amplification or attenuation of u if ν is chosen larger or smaller α, respectively.

Integration If, in contrast to the realization of a gain, the degradation term is chosen very small and the production term still is approximately linear with respect to a signal u, the input/output behavior rather resembles the one of an integrator. As small degradation rates are typical

for proteins, which in general are more chemically stable compared to mRNA, a realization with the process of translation seems favorable. If therefore thinking of M , the amount of mRNA, as input, assuming it to be piecewise constant and therefore zm = const and further δ ≈ 0, then we find P˙ = g(t, λ, η, Φ) − δP P˙ ≈ M ηc zm . Thus, the amount of protein P is proportional to the integral of M .

Positive and negative signals In order to represent both positive and negative signals with absolute numbers of certain biochemical species, Oishi and Klavins (2011) suggest to split signals into positive and negative parts, such that u = u+ − u− . As the individual signals u+ and u− can be scaled arbitrarily, the representation of a certain u therefore is not unique anymore. For example, u+ = 10 and u− = 11 represent the same signal as u+ = 0 and u− = 1. To avoid this, additional degradation reactions are introduced, e.g. (5) u+ + u− → Ø, therefore resulting in a minimal representation of a signal u which is unique (cf. Oishi and Klavins (2011) for details). In case of genetic networks however, representing each input and output signal of each function by distinct biochemical species can be quite restrictive in terms of resource consumption. An alternative to this conservative approach can be achieved when the distinction between positive and negative is shifted from the signal to the function space. For example, in the case of the differentiator, it suffices to generate two different genetic parts which track positive and negative slopes of a common input signal, respectively. Instead of splitting the input signal into two separate species, only the influence of this input on the two parts differs: while signal u activates the part which tracks positive slopes, it inhibits the one tracking negative slopes as depicted in Fig 4. At the same time, the output still splits signals into two: one acting as an activator, the other as a repressor for downstream parts. Again, to arrive at a minimal representation of this output, the two species should somehow inactivate each other, as shown in (5), e.g. when two mRNAs or proteins form a complex which either is inactive or degraded at a significantly higher rate compared to the single mRNA/protein. 3.3 Implementation of a genetic differentiator We use the just presented ideas to realize the genetic differentiator, consisting of a gain, an integrator and a negative feedback (cf. Fig.1), by choosing the degradation and initiation rates of the protein synthesis model appropriately. The environmental parameters Φ which mainly depend on the choice of the biological host have been chosen similarly as in Halter et al. (2017), viz. such that they represent the parameters of E. coli. The overall layout of a genetic differentiator is depicted in Fig 4. First, we only focus on a gene G1 which shall capture positive slopes of an input u and we choose its parameters similar to an average E. coli gene with approximately 1064 nucleotides and mRNA and protein half-lifes of ∼ 1min

11427

Proceedings of the 20th IFAC World Congress Wolfgang Halter et al. / IFAC PapersOnLine 50-1 (2017) 10938–10943 Toulouse, France, July 9-14, 2017

10941

u G1 u

Ø

y

DNA

f (t, λ(u, P ), Φ)

mRNA

g(t, λ, η, Φ)

ν

G2

δ

Ø Fig. 4. Genetic differentiator: Genes G1 and G2 tracking positive and negative slopes of u. Proteins produced by G1 and G2 neutralize each other. Associated mRNAs activate/inhibit the production of a reporter y. and ∼ 20h, respectively (see, e.g. (Halter et al., 2017)). Further, the basal transcription and translation initiation rates are set to 0 and 1e−3 events per minute, respectively. As pointed out before, these large differences between the degradation rates of mRNA and protein lead to different input/output functions of the processes of transcription and translation. The mRNA production therefore can be seen as a gain, while the protein production rather functions as an integrator. Fig. 5 depicts how the negative feedback and the input are connected for the gene to track positive slopes of u. The transcription initiation rate is linearly dependent on both u and the protein P . These interactions are positive for the input and negative for the protein and further need to be mutually exclusive, e.g. λ1 (u, P1 ) = −βP1 + u with β = 1e−5 chosen appropriately and the index 1 indicating the association with G1 . This results in the desired part with transcription representing the gain, translation the integrator and the mutually exclusive dependence of the initiation rate on both input and the own protein as the signal difference depicted in Fig. 1. The output of this part is the amount of mRNA or a protein expressed by this mRNA with a high degradation rate. As this amount cannot represent negative values, the given part can only capture positive slopes of u and it remains to develop a part for negative slopes. To do so, only two adaptations have to be made: the input u also has to act as an inhibitor for another gene which also incorporates a basal transcription initiation rate α. This gives e.g. λ2 (u, P2 ) = α − βP2 − u where we choose α = 5e−3, which determines the constant protein level of this gene for no input and therefore the maximum range of slopes this part can track. The assembly of the parts tracking positive and negative slopes is depicted in Fig. 4. In order to realize the minimal representation, a link between G1 and G2 is added. This direct interaction between the two genetic parts can be implemented either as mRNAmRNA interaction, protein-protein interaction or as both. Besides the purpose of a minimal representation, a direct interaction on the protein level however has the additional effect of counteracting the overflow of the integrator part and therefore is crucial for the long-term functionality of the differentiator module. In particular, as soon as the slope of u switches its sign, the integrator of the inactive part is being drained and thus the capacity is restored, i.e. the protein-protein interaction acts as anti-windup. We

protein

Ø

Fig. 5. Gene G1 tracking positive slopes of u. Input u activates transcription while protein P suppresses it. Amount of mRNA indicates positive slopes of u. therefore only implement the reaction k

→Ø P1 + P2 − with k = 0.02 chosen fast enough and neglect the mRNAmRNA interaction as it did not improve the overall performance of the module. The output of the overall system is now given by the difference of the amount of the mRNA of the two genes and to observe this output directly, we added a reporter gene y whose transcription is initiated by λy (M1 , M2 ) = αy + βy (M1 − M2 ) with αy = 1e−3 and βy = 1e−4 and whose protein is degraded with a very high rate to act as a direct indicator of the difference M1 −M2 . The basal rate αy causes a basal protein level P¯y to be expressed whenever M1 = M2 , thus the output we are interested in is given by the relative changes of Py with respect to the baseline P¯y , e.g. y(t) = Py (t) − P¯y . (6) The realization of such a reporter is not the focus of this work and might be done in several different ways, e.g. by directly expressing a protein from M1 and an according protease from M2 . Any mechanism indicating the difference between these two mRNAs is suitable. 4. EXAMPLES AND ANALYSIS We simulated the differentiation module together with an additional background gene which collects all the nominal activities and acts as the basal transcriptional and translational burden of an average E. coli cell, details of this background gene approach can be found in Halter et al. (2017). For evaluating the functionality of the differentiator, a test signal u(t) is used which consists of several components: a constant component of 1e−4, several ramps with slopes of 5e−5 min−1 and 1e−4 min−1 up to the time point t = 400min and subsequently a sinusoidal signal with an amplitude of A = 1e−4 and an angular frequency of ω = 5e−2 rad min−1 up to t = 800min. In order to compare the calculated derivative of the input and the reported derivative, it is crucial to normalize the signals to a meaningful range. As the amplitude of the derivative of u in the sinusoidal range is given by Aω, we use this value as normalization factor for the calculated derivative and define du 1 u˙ norm = · . dt Aω As in (6), the reported signal Py is first shifted by the amount of basal protein level P¯y and subsequently nor-

11428

Proceedings of the 20th IFAC World Congress 10942 Wolfgang Halter et al. / IFAC PapersOnLine 50-1 (2017) 10938–10943 Toulouse, France, July 9-14, 2017

·10−3

2

0 u(t) u˙ norm (t) ynorm (t) 0

100

200 time in min

-20 300

·10−4

1 0

−1 −2 400

500

600 time in min

700

−1 800

Fig. 6. Input signal u(t) (left axis), normalized calculated derivative u˙ norm (t) (right axis) and normalized reported derivative ynorm (t) (right axis). malized with respect to the sinusoidal part of the signal, namely Py (t) − P¯y ynorm (t) = . max (Py (t) − P¯y ) t∈[400,800]

Fig. 6 shows the input signal in solid blue, the normalized calculated derivative in dashed red as well as the normalized simulated output signal in solid red. We note that we split the time axis for better visualization of the different magnitudes of the signal. As expected, a time delay between the real and the genetically approximated derivative can be observed, additionally over- and undershoots are present at the discontinuities of the ramp signals. Still, the presented genetic differentiator is capable of approximating the real derivative of an input signal to a satisfying extent. We quantify the match between the normalized calculated and simulated derivative by calculating the mean error N 1  e= |ynorm,i − u˙ norm,i | N i=1

where ynorm,i = ynorm (ti ) and u˙ norm,i = u˙ norm (ti ) with {t1 , . . . , tN } the sampling time points of the signals and N the total number of sampling instances. The value of e represents the average deviation of the normalized simulated output from the normalized calculated derivative and in the presented case takes a value of 1.264. 4.1 Sensitivity and Robustness

In Section 2, we pointed out that the desired differentiating behavior is only achieved for sufficiently large gain K, or in other words, the approximation quality increases with increasing K. Consequently, to tune the behavior of the

0.1 590 5663 26852

.5 2148 4243 6302

1 3221 3221 3221

1.5 3868 2590 2163

2 4300 2165 1628

Table 2. Mean error e for different parameters Factor λc λBG ν1 and ν2

400

1

0

Factor λc λBG ν1 and ν2

normalized amount

0

−5

signal strength

20

normalized amount

signal strength

5

Table 1. Gain K for different parameters

0.1 4.107 0.913 3.977

.5 1.559 0.985 1.462

1 1.264 1.264 1.264

1.5 1.171 1.586 1.656

2 1.126 1.925 2.220

genetic differentiator, one has to amplify the production of mRNA of the two genes while ensuring that the fundamental input/output functions are not altered in their nature, i.e. a gain still behaves as a gain and so on. The three main possibilities to do so are increasing the transcription elongation rate λc , providing more RNAP for quicker transcription or decreasing the degradation rates ν1 and ν2 . While the first possibility might be difficult to achieve in particular for one specific gene, the second possibility is strongly dependent on the basal transcriptional activity of the cell. This reveals that in turn the performance of this option is dependent on the overall transcriptional load of the cell. The third option, changing the degradation rates, is critical as low values for ν are desired for a large gain but too small values lead to an integrating behavior instead of an amplifying one. This suggests that the degradation rates are good candidates for future in silico optimizations as there has to exist an optimal value for which both the function is assured and the gain is maximized. The simulation of the genetic differentiator and a background gene together enables us to study all three options for the effect of increasing the gain. By increasing the transcription initiation rate of the background gene λBG , we alter the available RNAP resources which remain for the genetic differentiator (cf. Halter et al. (2017)). We therefore changed the parameters λc , λBG and ν1 and ν2 relatively to the nominal values we used in Section 3.3 and calculated the gain K as the fraction of M1 , the amount of mRNA of G1 , and the constant input value which is applied at the beginning of the signal, thus M1 (t1 ) K= . u(t1 ) Table 1 and 2 summarize the values of gain K and mean error e for the variation of the parameters. While the K values increase with λc and decrease with λBG and ν1 and ν2 , the average error follows our previously outlined argumentation that a larger gain in general yields better performance of the genetic differentiator. For the degradation rates however, the gain cannot be chosen arbitrarily large as the amplifying behavior is lost which is also indicated by higher values of the mean error. The next step is to assess the robustness of the genetic differentiator. The above analysis suggests that variations in the cellular load may influence the performance significantly. To study the robustness towards all environmental parameters Φ, we further perturb all components of Φ and compare the mean error for these perturbed systems. For that purpose, we use the measure of total parameter variation (T P V ) introduced by Barkai and Leibler (1997) and given by

11429

κ

Proceedings of the 20th IFAC World Congress Wolfgang Halter et al. / IFAC PapersOnLine 50-1 (2017) 10938–10943 Toulouse, France, July 9-14, 2017

1 0.8 0.6 0.4 0.2 0

0

0.5

1 1.5 log(T P V )

10943

The derived genetic signal differentiator also paves the way for the genetic realization of more complex control and optimization techniques. A PID controller for instance may be implemented on a genetic level in order to introduce a control mechanism for synthetically designed molecular processes inside the cell. Also, one might think of implementing gradient based optimization algorithms for an individual self-optimization of cells, avoiding the ensemble control problem of a population of cells.

2

REFERENCES Fig. 7. Ratio κ of configurations with sufficiently good performance over total parameter variation.   L  Φi log(T P V ) = | log | Φnom,i i=1

with L the number of different environmental parameters and Φi and Φnom,i the i-th perturbed and nominal parameters, respectively. A T P V value of 10 for instance is obtained when one parameter is increased or decreased 10fold, √ or when two parameters are increased or decreased 10-fold and so forth. First, we define all configurations where the mean error e < 1.9 as sufficiently good. This limit was chosen heuristically. We then define κ(q) as the ratio between sufficiently good and all configurations for which log(T P V ) ∈ [q − , q + ]. The results are shown in Fig. 7 where we analyzed 3000 different parameter perturbations. For total parameter variations in the range 0 < log(T P V ) < 0.2, resulting in a maximum parameter change of roughly 58%, we still observe that 98% of the configurations lead to a satisfying performance. Although this ratio drops in the next range of 0.2 < log(T P V ) < 0.4 down to 72%, this analysis shows that the presented genetic differentiator performs robustly with respect to a certain extent of parameter uncertainty. 5. CONCLUSION In general, we adapted a design approach by Oishi and Klavins (2011) which now may be used for the conception of genetic regulatory networks realizing any arbitrary transfer function. Therefore, the desired transfer function is decomposed into basic linear I/O functions, each realized by certain biological mechanisms of protein production. In particular, we derived a set of rules for the synthesis of the I/O functions signal difference, gain and integrator and further illustrated that a signal differentiator can be approximately realized by a two-gene genetic network which resembles the system theoretic structure of a high gain with integral feedback. The input signal may be the concentration of any arbitrary molecular species while the output is given by the difference of two mRNA levels. This difference then indicates the temporal derivative of the input signal. We further demonstrated with an in silico approach that the resulting network indeed is functional and that it is also robust with respect to a certain level of parameter uncertainties. A rigorous in vitro proof of concept is currently missing. Such an experimental validation of our results however would yield interesting insights and shall be conducted in future studies.

Alon, U. (2007). An Introduction to Systems Biology Design Principles of Biological Circuits. Chapman & Hall/CRC. Ang, J., Bagh, S., Ingalls, B.P., and McMillen, D.R. (2010). Considerations for using integral feedback control to construct a perfectly adapting synthetic gene network. Journal of Theoretical Biology, 266(4), 723– 738. doi:10.1016/j.jtbi.2010.07.034. Barkai, N. and Leibler, S. (1997). Robustness in simple biochemical networks. Nature, 387(6636), 913–7. doi: 10.1038/43199. Briat, C., Zechner, C., and Khammash, M. (2016). Design of a Synthetic Integral Feedback Circuit: Dynamic Analysis and DNA Implementation. ACS Synthetic Biology, acssynbio.6b00014. doi:10.1021/acssynbio.6b00014. Del Vecchio, D., Dy, A.J., and Qian, Y. (2016). Control theory meets synthetic biology. Journal of The Royal Society Interface, 13(120), 20160380. doi: 10.1098/rsif.2016.0380. Halter, W., Montenbruck, J.M., Tuza, Z.A., and Allg¨ ower, F. (2017). A resource dependent protein synthesis model for evaluating synthetic circuits. Journal of Theoretical Biology. doi:10.1016/j.jtbi.2017.03.004. Harris, A.W.K., Dolan, J.A., Kelly, C.L., Anderson, J., and Papachristodoulou, A. (2015). Designing Genetic Feedback Controllers. IEEE Transactions on Biomedical Circuits and Systems, 9(4), 475–484. doi: 10.1109/TBCAS.2015.2458435. Lang, M. and Sontag, E. (2016). Scale-invariant systems realize nonlinear differential operators. In 2016 American Control Conference (ACC), 1, 6676–6682. IEEE. doi:10.1109/ACC.2016.7526722. Oishi, K. and Klavins, E. (2011). Biomolecular implementation of linear I/O systems. IET systems biology, 5(4), 252–260. doi:10.1049/iet-syb.2010.0056. Reuveni, S., Meilijson, I., Kupiec, M., Ruppin, E., and Tuller, T. (2011). Genome-scale analysis of translation elongation with a ribosome flow model. PLoS Computational Biology, 7(9). doi:10.1371/journal.pcbi.1002127. Yordanov, B., Kim, J., Petersen, R.L., Shudy, A., Kulkarni, V.V., and Phillips, A. (2014). Computational design of nucleic acid feedback control circuits. ACS Synthetic Biology, 3(8), 600–616. doi:10.1021/sb400169s.

11430