—The validity of a complex reaction pathway pro- posed to treat Inﬂammatory Bowel Disease (IBD) was veriﬁed by a comprehensive time and frequency domain analysis. The model was taken to the frequency domain to study the effect and the signiﬁcance of the negative feedback loop introduced by the reaction pathways. It could be shown that such proposed probiotics have very interesting potentials that could be used extensively in near future. Index Terms —Synthetic Biology, Control theory, Negative feedback in Biological systems I. I NTRODUCTION Inﬂammatory Bowel Disease (IBD) is characterised by chronic inﬂammation of the intestine. The condition is asso- ciated with an imbalance in immune cell populations, notably Th17 and Treg. Existing immunosuppressive therapies, when successful, often elicit systemic side effects and require fre- quent re-administration. Therefore, the proposed solution is a probiotic strain that restores the Th17/Treg cell balance via secretion of IL-10 in response to Nitric oxide in the intestinal lumen. Overshoot is prevented by an adenine riboswitch- sRNA construct which responds to extracellular adenosine, an indicator of the Treg cell population. The two independent BBSRC, Wellcome Trust, IDT, Eppendorf, BioLab, Snap gene and Oxford university departments of biochemistry, Engineering and chemistryreaction pathways - responding separately to IL-10 deﬁciency and abundance - and are summarised in Figure 1. To avoid the complexity of multivariate control systems, the pathways have been modelled independently. Physiologically, an elevated serum concentration of NO can be associated with autoimmunity, whereas an elevated serum adenosine con- centration can signify immunoinsufﬁciency. The NO pathway will be referred to as the ‘negative pathway’, representing a negative concentration difference in IL-10. Likewise, the adenosine pathway has been termed the ‘positive pathway’. Integration of separate stimuli in a dual feedback loop enables a more dynamic, robust response to the immune state of the body. Thus the engineered probiotic bacteria would be capable of maintaining an equilibrium between the Th17 and Treg cell populations by virtue of the negative feedback system. A. Negative Pathway The negative pathway can be modelled in a canonical manner: nitric oxide triggers oxidation of the SoxR iron- sulphur cluster, allowing binding to the promoter of IL-10 in the circuit - pSoxS. This results in transcriptional activation of the system within approximately 0:01seconds [1]. Hence, the system has been modelled by a simple Hill activation function for IL-10 gene transcription. The SoxR/pSoxS system has been studied extensively [2] but previous modelling has been limited Fig. 1. Simpliﬁed schematic diagram of the negative and positive pathways to that of the SoxR/pSoxS binding interaction by the 2009 Stanford iGEM team [3]. B. Positive Pathway The positive pathway mechanism facilitates translational in- hibition of IL-10 to prevent overshoot. The anti-inﬂammatory molecule, adenosine, is produced by Treg cells in response to inﬂammation, with elevated levels indicating Treg overactivity and an immunoinsufﬁcient state. Adenosine is subsequently broken down into ribose and adenine, via an extracytoplasmic, nucleoside hydrolase, by the engineered cell. Following uptake into the cell, adenine binds to a riboswitch, resulting in an increased rate of RNA transcription. A self-splicing ribozyme separates the riboswitch from the remaining transcript - an sRNA molecule which is complementary to the 5’ region of the IL-10 mRNA. Consequently, translation of IL-10 mRNA is prevented by the binding of sRNA. It should be noted that NO is an essential input to the positive pathway and hence, combined model for both reaction pathways could be found below in Figure‘2. Fig. 2. Schematic representation of the pathways [4] II. D YNAMIC MODEL The dynamic model of both of the pathways were analysed using the Simbiology toolbox in MATLAB as well as ODE15s Solver and some of the modelling and parameters are based on the paper ”Frequency domain analysis of small non-coding RNAs” [4]. The dynamic time domain analysis of the system was conducted using the combined modelin Figure 2 with different initial conditions corresponding to negative and positive pathways separately. It is worth mentioning that data ﬁtting was used to determine some of the parameters. A list of the parameters and data is available at the end of paper. The differential equations describing the systems are as follows: d[sRNA ] dt= 1 1[sRNA ]k1[IL-10 mRNA ][sRNA ];(1) d[IL-10 mRNA ] dt= 2 2[IL-10 mRNA ]k1[IL-10 mRNA ][sRNA ]; (2) d[IL-10 - i ] dt=k2[IL-10 mRNA ] 3[IL-10 - i ]k3[IL-10 - i ]; (3) d[IL-10 - s ] dt=k3[IL-10 - i ] 4[IL-10 - s ]; (4) where iare the Hill functions describing the rate of tran- scription with inputs of Adenosine and NO for 1and2 respectively as represented by uiin (5). i= [ui]n Kn+ [ui]n; (5) In which is the maximal transcription factor, Kthe disso- ciation coefﬁcient and nthe Hill coefﬁcient. k1is the rate of sRNA binding to IL-10 mRNA which heavily depends on the length of the sRNA, making it easy to exploit this relationship for better control over the fate of the system. k2 is the translation rate of IL-10 mRNA and k3is the secretion and diffusion rate of IL-10 [5]. irepresents the dilution + degradation rate. It is important to notice that Hill functions introduce non- linearity into the system as well as sRNA binding. Therefore, Hill functions have been linearised for transfer function deriva- tion in the frequency domain, and are noted by as follows:  i= iniKni i[ui]ni1 (Kni i+uni i)2; (6) All the assumptions made for (1) to (4) are listed below: Adenosine substituted with adenine as the hydrolase reaction is believed to be much faster than the body response Concentration of adenine and NO kept constant for dy- namic analysis due to intra- and extracellular abundance The initial conditions used for time domain analysis correspond to maximal insufﬁciency based on [6], [7] and [8]. Stochastic response ignored due to the large number of E. Coli used. A list of all parameters and concentrations used in ﬁgures, charts and graphs can be found at the end of the section. A. Negative pathway dynamics Negative pathway dynamics are fully described by the ODEs stated previously. It should be noted that based on the second assumption, the level of NO would not change with time and henced[NO] dt= 0. An initial concentration of 19:88M was used for elevated level of NO corresponding to IL-10 deﬁciency in patients with IBD [6] as well as the nominal concentration of 18M for adenine [8]. The evolution of the system to NO stimuli is illustrated in Figure 3. 0 50 100 150 20005101520NO response dynamics Adenine sRNA IL10 - intracellular IL10 mRNA IL10 - secreted NO Fig. 3. NO response dynamics with nominal level of Adenine As it can be seen from Figure 3, [IL-10 ]reaches a concen- tration of 14:16M after roughly 180 seconds. It is important to note that the initial concentration of IL-10 is set to 0 for simplicity, meaning that in reality, the model can correct deﬁciencies of IL-10 to a margin of 14:16M. A control mechanism is essential in order for the probiotic bacteria to be able to maintain a healthy level of IL-10. To determine how system would respond and whether or not it would be able to stabilise the IL-10 concentration, it is necessary to construct a large-scale version of the model that incorporates the body as well as the system embedded in a negative feedback loop [9]. The control modelling is discussed in more detail later in the paper. Lastly, the model behaviour is as expected due to the constant production rate, as opposed to the degradation rate, and hence a bounded steady state value is expected. B. Positive pathway dynamics The complexity of the positive pathway is visible from Figure 2. The complexity of the system was maintained for the dynamic analysis based on ODEs stated previously. It should be noted that based on the second assumption, the concentration of adenine would not change with time and henced[Adenine ] dt= 0. An initial concentration of 100M was used for elevated level of adenine corresponding to IL-10 abundance in patients with IBD [8] as well as the nominal concentration of 13:24M for NO [6]. The evolution of the system to adenine stimuli is illustrated in Figure 4. It can be observed from Figure 4, that [IL-10 ]is inhibited signiﬁcantly and it could be assumed that a high concentration of adenine could inhibit translation of IL-10 completely. It should be noted that the length of sRNA and the promoter strength in both negative and positive pathways( k1; 1; 2) play a signiﬁcant role and hence they have been optimised for 0 50 100 150 200020406080100Adenine response dynamics Adenine sRNA IL10 - intracellular IL10 mRNA IL10 - secreted NOFig. 4. Adenine response dynamics with nominal level of NO maximum inhibition as well as maximum translation of IL- 10 in both positive and negative pathways respectively. The sensitivity analysis and sRNA models are discussed below. III. S ENSITIVITY ANALYSIS AND S RNA MODELLING Sensitivity analysis was done separately for each reaction pathway in order to ﬁnd the most signiﬁcant steps and pa- rameters, thus enabling these parameters to be exploited in order to improve the responsiveness of the system to smaller perturbations in concentrations. The sensitivity of [IL10] with respect to parameters used in simulations is shown in Figure 5. Fig. 5. Sensitivity Analysis of [IL-10] with respect to parameters It is clear from Figure 5 that IL-10 production depends heavily on 2and 1which is desirable as different promoter strengths (weak, medium and strong) would enable optimisation of the response. It can also be observed that sRNA binding rate kf has less signiﬁcant effect on the output compared to 2and hence, a stronger promoter would mean that a much longer sRNA would be required to inhibit IL-10 translation. The length of sRNA was optimised to 24 base pairs [10] yielding G=45:97Kcalmole1. The graph of Gagainst the position of base pairs is illustrated in Figure 6. Based on the calculated G, 1, 2andk1were set to 5 10 15 20 25 30-60-50-40-30-20-100sRNA Binding EnergyFig. 6. sRNA binding free energy variation with its base pair position 1nMmin1,10nMmin1and100nM1min1for com- plete inhibition of IL-10 translation in the positive pathway, in addition to a high IL-10 translation rate in negative pathway. IV. S TEADY STATE MODEL The simplest model is an output/input, time-independent plot which describes the system fate in response to varying input concentrations. Secretion of IL-10 has been ignored for steady state response as it is assumed that the system is given ample time to settle, secreting all translated IL-10. The steady state concentrations of species can be determined by equating (1), (2) and (3) to zero. The steady state values for [IL-10], [sRNA] and [IL-10 mRNA] are given below: [IL-10 mRNA ]= 2 1+k1(12)p! 2 2k1; (7) [sRNA ]=1 1+k1[IL10 mRNA ]; (8) [IL-10 ]=k2 4[IL-10 mRNA ]; (9) where!= k1(21) 1 22+ 4 1 2k12. It should be noted that idenotes Hill functions for adenine and NO. A surface can be plotted describing the output/input steady state relation, with contours representing the steady state characteristic of each of the pathways. The plots can then be used for experimental ﬁtting to obtain more realistic parameters. The steady state surface and it’s contours could be found in Figure 7. V. P ARAMETERS AND MODELLING PREDICTIONS FOR CIRCUIT DESIGN The parameters used for modelling are all listed below: Transcription level : 1= 10nM=min , 2= 1nM=min both optimised to give satisfactory results [4]. K1= 300nM [11],K2= 10nM [4]. Translation level : k1= 100min1optimised based on the number of base pairs in sRNA. k2= 0:3min1is the translation rate of IL-10 [4]. Degradation and Diffusion : 1= 0:03min1, 2= 0:14min1, 3= 4= 0:03min1all based on [4]. Fig. 7. The steady state surface and its contours to both of the inputs The diffusion rate parameter was calculated based on the number of transport proteins on the membrane of E. Coli to bek3= 0:13min1[12]. Body response : The value of b1andb2were calculated from Figure 8 and a chosen value of a= 1e4s.b1= 1:6690 andb2= 28:24. A. Parameters range and output dependency Three of the parameters used throughout the paper can be tuned for optimisation, namely: 1, 2andk1; some of which have been previously mentioned. These parameters were modiﬁed slightly about their optimised value to ﬁnd a range for which the system maintains a desired working state. The optimised range for each of the parameters is given below. 1) 1has a wide range of 8< 1<90for which a satisfactory results could be obtained. 2) As opposed to 1, 2has a very narrow range of 0:7< 2<1:3for which a satisfactory results could be obtained. 3)k1is almost saturated at the value chosen for simulation and hence does not have an upper bound; however, there is a lower bound of 1<k1for which a satisfactory result can be obtained. Equations 7 - 9 can be differentiated to ﬁnd the dependency of IL-10 with respect to all three tunable parameters. As [IL10] depends linearly on [IL10 mRNA ], (7) has been differentiated with respect to 1,2andk1for simplicity. It should be noted that since iis the only tunable parameter in Hill functions, i has been used as the parameter for simplicity in differentiation. d[IL-10 mRNA ] d1=d d1 2 1+k1(12) +p! 2 2k1 )=k1 2 2k11p! k1(21) 1 2 1 ; (10) d[IL-10 mRNA ] d2=d d2 2 1+k1(12) +p! 2 2k1 )=k1 2 2k1 11p! k1(21) 1 2 ; (11) d[IL10 mRNA ] dk1=d dk1 2 1+k1(12) +p! 2 2k1 )= (12) + (2 2k1)2 2 (2 2k1)2; (12) where= k1(12) 1 2 (12 +4 1 22p!and= 1 2+k1(12) +p! . Based on (10), (11) and (12), the dependency of IL-10 with respect to 1, 2andk1is clear. The expression k1(21) 1 2is less thanp!and therefore it could be concluded that: 1)d[IL-10 mRNA ] d1<0 2)d[IL-10 mRNA ] d2>0 3)d[IL-10 mRNA ] dk1<0 1andk1are both in the inhibitory pathway as opposed to 2and hence, all three statements above validate the model. B. Modelling predictions for circuit design Based on the range and dependency of [IL-10 ]on different parameters, a medium strength promoter is suggested to be used for negative pathway as opposed to positive pathway where either a strong promoter or a medium strength promoter with a high copy number is suggested. As for the length of sRNA, it was shown that as long as it has just a few base pairs, corresponding to k1= 1, the output would not be affected signiﬁcantly and hence, the most energetically feasible length is used. VI. B ODY DYNAMICS The response of the body is crucial to obtain comprehensive model reﬂecting systemic interactions between the engineered E. coli and the body. Body response mechanisms have been studied in detail however due to stochastic noise, signiﬁcant assumptions were made to simplify the response for modelling. These assumptions are summarised below; please note that even though some of the assumptions are considerable, they are negligible compared to the amount of variation present in bodily systems and between individuals:The model is based on bulk amount of bacteria translating or inhibiting IL-10. The body responds in about 3-4 hours [] and the response has been taken to be linear due to lack of data points. The body has been assumed to behave similarly to an LTI (Linear time-invariant) system and the model can be extrapolated, with care, for people with different body responses. Based on assumptions stated above, two lines can be plotted for elevated levels of NO and Adenosine versus IL-10 - the highest and lowest concentration of IL-10 [13] and its corre- sponding signalling molecule [6] [8] that has been measured in human body. The plot of steady-state concentrations is shown in Figure 8. 0 1 2 3 4 5 6 7 [IL10]/µM ×10-6020406080100[NO,Adenine ]/µMSteady state body response Fig. 8. Linear approximation for steady-state body response in both reaction pathways The dynamic response of the body is assumed to be linear with a response time of about 3 hours (104s)as very few data points were available for the steady state response of the body. Therefore, the differential equation describing the response is as follows: d[NO,Adenine ] dt=ai[NO,Adenine ] +bi[IL10] +ci;(13) Figure 8 can be used to ﬁndb aandc afor both positive and negative pathways. Then by choosing a= 104s1, for the response time, (13) is fully deﬁned and can be solved simultaneously with (1) - (4) to give an overall prediction on how the body and the probiotic bacteria would interact with each other to reach an equilibrium concentration of IL-10 in the body. It should be noted that for these two to be solved together, one should be aware of the correction factor required to approximate metabolite concentrations inside a bacterium relative to the body, in addition to an indication of the number of bacteria being used. Hence, the initial number of bacteria is optimised in a correction loop by looking at the settle time of IL-10 inside the body. The dynamics of IL-10 in the body for elevated level of NO and adenine are shown in Figures 9 and 10 respectively. It should be noted that for representation purposes, the concentration of IL-10 in the body - magenta dashed line - has been multiplied by 106. It is clear from Figures above that the body and E. Coli reach a stable concentration of IL-10 without oscillation in any of the metabolite concentrations. The accuracy of the response time 0 500 1000 1500 2000050100150200250300350Combined Dynamic Response to Elevated NO sRNA mRNA IL10 - E Coli IL10 - Body NO AdenineFig. 9. E. Coli - Body dynamic response to elevated level of NO 0 500 1000 1500 2000020406080100Combined Dynamic Response to Elevated Adenine sRNA mRNA IL10 - E Coli IL10 - Body NO Adenine Fig. 10. E. Coli - Body dynamic response to elevated level of adenine of104seconds is reﬂected by the similarity of the steady-state. The number of E Coli used for simulation would signiﬁcantly affect the settle time or the response time and hence, the time taken by the combined system to reach 90% of the steady-state concentration of IL-10 has been plotted against the number of bacteria in Figure 11. It is essential to understand the signiﬁcance and importance of the positive pathway in stability and robustness of the proposed probiotic. Hence, dynamics of the combined system has been plotted in absence of Adenine pathway to demonstrate that the system becomes unstable without a negative feedback loop. 10101015100101102103Settle Time - Number of E Coli Fig. 11. 90 % Settle Time - Number of E Coli used for Simulation 0 2 4 6 8 10 1040100200300400500600Unstable system without the proposed feedback loop sRNA mRNA IL10 - E Coli IL10 - Body NO AdenineFig. 12. Unstable system in the absence of the inhibiting pathway VII. N EGATIVE FEEDBACK AND FREQUENCY DOMAIN ANALYSIS As discussed earlier in the paper, both positive and negative feedback could potentially correct the concentration of IL- 10 in the body. However, it is important to note that body response should be included in the model to give an actual and meaningful prediction of how E. coli would behave in the body. For this purpose, control theory and the negative feedback loop have been used by transforming the body and the bacteria in a standard cascade compensation model and deriving their transfer functions. Hence, we developed a cascade compensation model where we have the body as the plant and the bacteria as the controller. Assumptions made for frequency domain analysis are listed below: The system has been linearised about an equilibrium point and therefore, it has been assumed that perturbations are small for the Taylor series to converge. The reaction pathways will be analysed separately as based on superposition, the response of the linearised model would be a sum of its response to inputs separately. Secretion has been ignored for simplicity in transfer function. Therefore, by using the correct transfer function, both models would fall into the same control design, just with different transfer functions. The proposed design for the feedback loop could be seen in Figure 13. Controller E coli G(s) Referrence(s) NO/Adenine + - Error(s) elevated NO/Adenine Plant Body F(s) IL10(s) NO/Adenine (s) Fig. 13. Cascade compensation design The aim of the controller is to set the reference signal, namely NO and adenosine, to nominal level as this would correspond to a healthy concentration of IL-10 within the body. Transfer functions for the controller- both positive and negative pathway- and the body are required for control analysis of system response in a human body. VIII. T RANSFER FUNCTIONS As stated in the assumptions, the model has been linearised due to non-linearity raised by the Hill functions and sRNA binding. Therefore, Hill functions have been replaced by  ias stated earlier in Equation 6. The block diagram and signal ﬂow graph for the linearised system are shown in Figures 14 and 15, where transcription of IL-10 mRNA, sRNA and translation of IL-10 can be seen clearly. γ⋅11 sk1[mRNA]⋅ γ⋅21 s α2 k1[mRNA]⋅k1[sRNA]⋅α1 k1[mRNA]⋅ k21 sAdenine(s) NO(s) + + + + + - + + - - - - - - sRNA(s) mRNA(s) IL10(s) α3 Fig. 14. Block diagram in frequency domain for linearised model NO(s)γ⋅21 sk2 IL10(s)α2 k1[sRNA]⋅ α3 γ⋅11 s sRNA(s)α1 k1[mRNA]⋅Adenine(s)mRNA(s) k1[mRNA]⋅k1[sRNA]⋅1 s Fig. 15. Signal Flow Graph(SFG) for the linearised model Transfer functions for both negative and positive pathways were found using Mason’s gain relation [14] and conﬁrmed against [4]. Transfer functions of both positive and negative pathways can be found below: G1(s) =IL10(s) Adenine(s)=k2  1(s+ 1+k1[IL10 mRNA ]) (s+ 3)(s+s+)(s+s); (14) G2(s) =IL10(s) NO(s)=k1[IL10 mRNA ]k2  2 (s+ 3)(s+s+)(s+s); (15) wheres=1 2( 2+k1[sRNA ]+ 1+k1[IL10 mRNA ]) 1 2p dandd= ( 2+k1[sRNA ] 1k1[IL10 mRNA ])2+ 4k2 1[IL10 mRNA ][sRNA ]. A. Body response transfer function Laplace transform can be taken from both sides of the Equa- tion 13 and by ignoring c, being very small and negligible, the transfer function for the body can be derived: s[NO,Adenine(s) ] =ai[NO,Adenine(s) ] +bi[IL10(s) ]=)Fi(s) =[NO,Adenine(s) ] [IL10(s) ]=bi s+ai(16) IX. S TABILITY The closed loop transfer function for cascade cascade con- troller design given in Figure 13 would contain important information regarding steady state error and relative stability. It should be mentioned that similar to most control problems where the controller is designed for speciﬁc speciﬁcations [15], the reaction pathway can be changed by altering promoter strength and sRNA length, making it almost possible to adopt desired structures for the controller. The closed loop transfer function for positive and negative pathways is as follows: Ti(s) =Fi(s)Gi(s) 1 +Fi(s)Gi(s)(17) Based on (17), the stability of the system can be determined by observing the location of the poles which are roots of the characteristic polynomial 1 +Fi(s)Gi(s) = 0 . The character- istic equation for adenine and NO reference signals is given in Equations 18 and 19. 1 +k2  1(s+ 1+k1[IL10 mRNA ]) (s+ 3)(s+s+)(s+s)b1 s+a1= 0 (18) 1 +k1[IL10 mRNA ]k2  2 (s+ 3)(s+s+)(s+s)b2 s+a2= 0 (19) Equations (18) and (19) were solved numerically [16] and smallest poles and zeros for both reaction pathways are summarised in the Argand diagram below, Figure 16. All poles and Zeros are summarised after the ﬁgure. Fig. 16. poles and zeros for NO and Adenine responses It should be noted that some poles and zeros were not plotted as they were signiﬁcantly larger(in magnitude) than those shown in Figure 16 and would have thus hindered comparison with those illustrated. All poles and zeros are given below. 1) NO response : Poles :s11:075104,s20:030 ,s3 0:030ands41:060104. Zeros :z1= 1:075104,z2=0:030,z3= 0:030,z4= 1:000104. 2) Adenine response : Poles :s150:656,s20:474,s30:474 and s4=0:262. Zeros :z1= 50:656,z2= 50:000,z3= 0:656, z4= 1:000104. Despite right-plane poles, it should be noted that both re- sponses can be considered stable as the dominant poles are signiﬁcantly larger- in magnitude- than all other poles. X. S TEADY STATE ERROR - DC GAIN E(s) =1 1 +Gi(s)Fi(s)R(s); R(s) =C sC2R FVT : lim t!1e(t) = lim s!0sE(s) )lim t!1e(t) = lim s!01 1 +Gi(s)Fi(s)(20) By substituting back Transfer functions for both E. Coli, Gi(s), and the body, Fi(s), the steady state error can be found as follows: For NO response limt!1e(t) =2:9501104. For adenine response limt!1e(t) =2:1912103. It is evident that the combined model has a very small DC gain and hence, a very small margin of error is estimated for the steady concentration of IL-10. ACKNOWLEDGEMENT Supervision and guidance of the following people were signiﬁcant in achieving and ﬁnishing this paper. Their con- tributions and help is greatly appreciated. Karandip Saini - Department of Chemistry, our team wet Lab coordinator Dr. Aivar Sootla - Department of Engineering sciences University of Oxford. Dr.Nicolas Delalez - Department of Engineering sciences University of Oxford. Prof.Antonis Papachristodoulou - Department of Engi- neering sciences University of Oxford. Prof.Christopher Macminn - Department of Engineering sciences University of Oxford. Dr.George Wadhams - Department of Biochemistry Uni- versity of Oxford.