Team:Oxford/Model/Time-Domain-Analysis
Time Domain Analysis
ODE based model
The dynamic model of each pathway was analysed using ODE15s in Matlab. Some parameters were sourced from "Frequency domain analysis of small non-coding RNAs shows summing junction-like behaviour" (Steel et al. 2017). The dynamic time domain analysis of the system was conducted using the combined model represented here with different initial conditions corresponding to negative and positive pathways separately. It is important to note that data fitting was used to determine some of the parameters. Differential equations describing the systems are as follow.
- $$\frac{d[\text{sRNA}]}{dt} = \Gamma_1 - \alpha_1[\text{sRNA}] - k_1[\text{IL10 mRNA}][\text{sRNA}], \;$$
- $$\frac{d[\text{IL10 mRNA}]}{dt} = \Gamma_2 - \alpha_2[\text{IL10 mRNA}] - k_1[\text{IL10 mRNA}][\text{sRNA}], \;$$
- $$ \frac{d[\text{IL10 - i}]}{dt} = k_2[\text{IL10 mRNA}] - \alpha_3[\text{IL10 - i}] - k_3[\text{IL10 - i}], \;$$
- $$\frac{d[\text{IL10 - s}]}{dt} = k_3[\text{IL10 - i}] - \alpha_4[\text{IL10 - s}], \;$$
- $$1.\; [\text{IL10 mRNA}]^* = \frac{\alpha_2 \alpha_1 +k_1(\Gamma_1 -\Gamma_2) \pm \sqrt{\omega}}{-2\alpha_2k_1},$$
- $$2.\; [\text{sRNA}]^* = \frac{\Gamma_1}{\alpha_1 +k_1[\text{IL10 mRNA}]^*},$$
- $$3.\; [\text{IL10}]^* = \frac{k_2}{\alpha_4}[\text{IL10 mRNA}]^*,$$
where \(\Gamma_i\) denotes the Hill functions describing the rate of transcription with inputs of Adenine and NO for \(\Gamma_1\) and \(\Gamma_2\) respectively as represented by \(u_i\) below. $$\Gamma_i = \frac{\beta [u_i]^n}{K^n+[u_i]^n}$$
In which \(\beta\) is the maximal transcription factor, \(K\) the dissociation coefficient and \(n\) the Hill coefficient. \(k_1\) is the rate of sRNA binding to IL10 mRNA which heavily depends on the length of the sRNA, making it easy to exploit this relation for a better control over the fate of the system. \(k_2\) is the translation rate of IL10 mRNA and \(k_3\) is the secretion and diffusion rate of IL10. \(\alpha_i\) also represent the dilution + degradation rate.
NO response dynamics
The negative pathway dynamics is fully described by the ODEs stated earlier and it should be noted that based on the second assumption, the concentration of NO would not change with time. An initial concentration of \(19.88 \mu M\) was used to simulate an elevated level of NO corresponding to IL10 deficiencies in patients with IBD, as well as the nominal concentration of \(18 \mu M\) for Adenine. The evolution of the system in response to NO stimuli is illustrated in Figure 1.
Adenine response dynamics
The complexity of the positive pathway can be seen here. The system was kept as detailed as possible for the dynamic analysis based on ODEs stated earlier. It should be noted that based on the second assumption, the concentration of Adenine would not change with time. An initial concentration of \(100 \mu M\) was used to represent an elevated level of Adenine corresponding to IL10 abundance in patients with IBD, as well as the nominal concentration of \(13.24 \mu M\) for NO. The evolution of the system to Adenine stimuli is illustrated below in Figure 1.
As it could be seen from Figure - 1, \([IL10]\) reaches a concentration of \(14.16 \; \mu M\) after roughly 180 Minutes. it is important to note that the initial amount of IL10 is set to 0 for simplicity in the model meaning that in reality, the model can correct deficiencies of IL10 to a margin of \(14.16 \mu M\). For bacteria to be able to follow and maintain a healthy level of IL10, a control mechanism and design is essential. To develop notions into how the system would respond and whether or not it would be able to stabilise the concentration of IL10, we need to construct a bigger picture of the model that has the body as well as the system embedded in a negative feedback loop. The control modelling is discussed in more details later on the page. Lastly, the model behaviour is as expected as the production rate is constant as opposed to the degradation rate and hence a bounded steady state value is expected.
It could be observed from Figure-1, that \([IL10]\) is inhibited significantly and it could be assumed that a high concentration of adenine could inhibit translation of IL10 completely. It should be noted that the length of sRNA and the promoter strength in both negative and positive pathways(\(k_1\, ,\, \beta_1 \, ,\, \beta_2\)) play a significant role and hence they have been optimised for maximum inhibition as well as maximum translation of IL10 in both positive and negative pathways respectively.
The simplest model is an output/input, a time-independent plot which describes the system fate in response to different concentrations of inputs. The secretion of IL10 has been ignored for the steady-state response as it is assumed that the system is given more than enough time to settle, secreting all translated IL10. The steady-state concentrations of species could be easily found by equating ODEs stated above to zero. The steady state values for [IL10], [sRNA] and [IL10 mRNA] are given below.
where \(\omega = \big(k_1(\Gamma_2-\Gamma_1)-\alpha_1\alpha_2\big)^2+4\alpha_1\alpha_2 k_1 \Gamma_2\). It should be noted that \(\Gamma_i\) are Hill functions for Adenine and NO. A surface could be plotted describing the output/input steady state relation as well as contours describing the steady state characteristic of each of the pathways. The plots then could be used for experimental fitting for more realistic parameters. The steady state surface could be found in Figure- 2. Contours of the surface are shown in Figures 3.
Sensitivity analysis and Optimisation
Sensitivity analysis was done separately for each reaction pathway to determine the most significant steps and parameters, 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 4.
It is clear from Figure-4 that IL10 production depends heavily on \(\beta_2\) and \(\beta_1\) which is desirable as they could be easily optimised by using different promoters with different strength(weak, medium and strong). It could also be observed that sRNA binding rate, \(k_f\), has less considerable effects on the output compared to \(\beta_2\) and hence, a stronger promoter would mean that a much longer sRNA is required to inhibit IL10 translation.
The length of sRNA was optimised to 24 base pairs[G], which yielded \(\Delta G = -45.97 Kcal\cdot mole^{-1}\). A graph of \(\Delta G \) against the position of base pairs is illustrated in Figure- 5
Based on the calculated \(\Delta G\), \(\beta_1\), \(\beta_2\) and \(k_1\) were set to \(1 nMmin^{-1}\), \(10 nMmin^{-1}\) and \(100 nM^{-1}min^{-1}\) for complete inhibition of IL10 in positive pathway as well as high translation rate of it in negative pathway. Parameters range, output dependancy and modelling predictions for circuit design is discussed below.
Parameters range and output dependency
Three of the parameters used throughout the model could be tuned for optimisation namely, \(\beta_1\), \(\beta_2\) and \(k_1\); some of which have been previously mentioned. These parameters were changed slightly about their optimised value to find a range for which our system keeps behaving in a desired manner. The optimised range for each of the parameters is given below.
- \(\beta_1\) has a wide range of \(8 < \beta_1 <90\; nMmin^{-1}\) for which a satisfactory results could be obtained.
- As opposed to \(\beta_1\), \(\beta_2\) has a very narrow range of \(0.7<\beta_2<1.3 nMmin^{-1}\) for which a satisfactory results could be obtained.
- \(k_1\) is almost saturated at the value chosen for simulation and hence does not have an upper bound; however, there is a lower bound of \( 1 < k_1 \) for which a satisfactory results could be obtained.
Steady-State Equations could be differentiated to find dependency of IL10 with respect to all three tunable parameters. As \([\text{IL10}]^*\) depends linearly on \([\text{IL10 mRNA}]^*\), Equation 2 has been differentiated with respect to \(\Gamma_1,\; \Gamma_2\)and \(k_1\) for simplicity. It should be noted that since \(\beta_i\) is the only tunable parameter in Hill functions, \(\Gamma_i\) has been used as the parameter for simplicity in differentiation.
- $$ \frac{d[\text{IL10 mRNA}]^*}{d\Gamma_1} = \frac{d}{d\Gamma_1}\big(\frac{\alpha_2 \alpha_1 +k_1(\Gamma_1 -\Gamma_2) + \sqrt{\omega}}{-2\alpha_2k_1}\big) \\ \therefore \quad = \frac{k_1}{2\alpha_2k_1}\bigg(\frac{1}{\sqrt{\omega}}\big(k_1(\Gamma_2 - \Gamma_1)-\alpha_1\alpha_2\big)-1\bigg), $$
- $$\frac{d[\text{IL10 mRNA}]^*}{d\Gamma_2} = \frac{d}{d\Gamma_2}\big(\frac{\alpha_2 \alpha_1 +k_1(\Gamma_1 -\Gamma_2) + \sqrt{\omega}}{-2\alpha_2k_1}\big), $$
- $$\therefore \quad = \frac{k_1}{2\alpha_2k_1}\bigg(1-\frac{1}{\sqrt{\omega}}\big(k_1(\Gamma_2 - \Gamma_1)-\alpha_1\alpha_2\big)\bigg), $$
- $$\frac{d[\text{IL10 mRNA}]^*}{dk_1} = \frac{d}{dk_1}\big(\frac{\alpha_2 \alpha_1 +k_1(\Gamma_1 -\Gamma_2) + \sqrt{\omega}}{-2\alpha_2k_1}\big), $$
- $$ \therefore \quad = \frac{\bigg((\Gamma_1 - \Gamma_2)+ \lambda \bigg)(-2\alpha_2k_1)-2\alpha_2\mu}{(-2\alpha_2k_1)^2} , $$
where \(\lambda = \frac{\big(k_1(\Gamma_1-\Gamma_2)-\alpha_1\alpha_2\big)(\Gamma_1 - \Gamma_2\big)+4\alpha_1\alpha_2\Gamma_2}{\sqrt{\omega}} \) and \(\mu =\big(\alpha_1\alpha_2+k_1(\Gamma_1 - \Gamma_2)+\sqrt{\omega}\big) \). Based on Equations , the dependency of IL-10 with respect to \(\beta_1\), \(\beta_2\) and \(k_1\) is clear. The expression \(k_1(\Gamma_2-\Gamma_1)-\alpha_1\alpha_2\) is less than \(\sqrt{\omega}\) and therefore it could be concluded that:
- $$\frac{d[\text{IL-10 mRNA}]^*}{d\Gamma_1} <0 $$
- $$\frac{d[\text{IL-10 mRNA}]^*}{d\Gamma_2} >0$$
- $$\frac{d[\text{IL-10 mRNA}]^*}{dk_1} <0$$
\(\beta_1\) and \(k_1\) are both in the inhibitory pathway as opposed to \(\beta_2\) and hence, all three statements above validate the model.
Based on the range and dependency of \([\text{IL-10}]\) on different parameters, a medium strength promoter is suggested to be used for the negative pathway as opposed to the positive pathway where either a strong promoter or a medium strength promoter with a high copy number plasmid is suggested. As for the length of sRNA, it was shown that as long as it has just a few base pairs, corresponding to \(k_1 =1\), the output would not be affected significantly and hence, the most energetically feasible length is used.
Body Dynamics and Combined Model
The response of the body is crucial for a comprehensive model and understanding how our engineered E.coli would interact with the body. The body response mechanisms have been studied in detailed however due to stochastic behaviour of the body as well as the complexity of each and every mechanism, the response has not been modelled very well in the literature and therefore, significant 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 randomness present in the human body and variations between individuals.
- The model is based on a bulk amount of bacteria translating or inhibiting IL10 .
- item 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 could be extrapolated, with care, for people with different body responses
Based on assumptions stated above, two lines could be plotted for elevated levels of Nitric Oxide and Adenine against IL10- the highest and lowest concentration of IL10 and it's corresponding signalling molecule that has been measured in the human body. The plot of steady-state concentrations is shown in Figure 6.
The dynamic response of the body is assumed to be linear with a response time of about 3 hours\((10^{4} s)\) as only very few data points were available for the steady-state response of the body. Therefore, the differential equation describing the response is as follow.
$$\frac{d[\text{NO,Adenine}]}{dt} = -a_i[\text{NO,Adenine}]+b_i[\text{IL10}]+c_i, $$Figure 6 could be used to find \(\frac{b}{a}\) and \(\frac{c}{a}\) for both positive and negative pathways. Then by choosing \( a = 10^{-4} s^{-1}\), for the response time, Equation above is fully defined and could be solved simultaneously with Equations 1-4 to give an overall prediction as to how the body and our engineered E.coli would interact with each other to reach a healthy steady concentration of IL10 in the body. It should be noted that for these two to be solved together, one should be aware of the correction required for concentrations inside one bacteria and inside the body as well as the number of bacteria that is being used. Hence the initial number of bacteria is optimised in a correction loop by looking at the settling time of IL10 inside the body. The dynamics of IL10 in the body for the elevated level of NO and adenine could be found in Figure 7. It should be noted that for the sake of presentation the concentration of IL10 in the body - magenta dashed line - has been multiplied by \(10^{6}\).
It is clear from Figures above that the body and E.coli reach a stable and steady concentration of IL10 without any instability and oscillation in any of the chemicals. It should also be noted that even though the response time seems to be rather long, it is comparable to existing medications for IBD and It should also be noted that one would be able to reach a much faster response time by increasing the number of E.coli used. To demonstrate the significance of the existence of both of the pathways for a controlled response, the following plot illustrates the response of the combined system in the absence of an inhibiting pathway. The response time is also plotted against the number of E.coli used to show how the speed of the response would depend on the number of bacteria used. Both Figures have been given below.
