Team:Oxford/Model/Frequency-Domain-Analysis
Model
Modeling Description
Various features of the project have been modelled to verify the validity of our proposed solution and to enable us to make predictions based on different aspects of the system. Firstly, nonlinear ordinary differential equations (ODEs) were written for all reactions taking place. These were then solved numerically over 250 seconds to characterise both the transient and steady-state response of the system. The time-domain analysis of the system was proceeded by looking at steady-state curves and output/input behaviour of the system in the steady state as well as body response dynamics, which gave an overall prediction about the fate of the combined system-body model. Secondly, optimisation and modelling were used to determine the optimum promoter strength and number of base pairs needed for sRNA binding. Modelling was then taken to the frequency domain for transfer function derivation and cascade controller design.
The detailed reaction pathway is shown in Figure 1.
Methodology
A list of general assumptions is summarised below. It should be noted that more specific assumptions are stated where they have been used for simplification.
- Adenosine substituted with adenine as the hydrolase reaction is believed to be much faster than the bodily response.
- Adenine and NO concentrations kept constant for dynamic analysis due to the intra- and extracellular abundance.
- Initial conditions used for time domain analysis correspond to maximal immunodeficiency.
- The dynamic model has been developed based on the assumption that all reactions are irreversible and they all assumed to behave as if they were taking place in a cell-free reaction vessel. It should be noted that later in our model we introduce correction factors to take into account concentration differences inside and outside of the cell.
- The stochastic response has been ignored due to the large number of E Coli that is used.
The following biochemical reactions have been modelled for both time domain and frequency domain analysis.
Parameters used in modelling are all listed in the table below.
| Parameters | Description | Value from Literature | Value used | Units | Optimised | Refference |
|---|---|---|---|---|---|---|
| \(\beta_1\) | Maximal Transcription rate | \(1\) | \(10\) | \(nMmin^{-1}\) | Yes | A |
| \(K_1\) | Dissociation constant | \(300\) | \(300\) | \(nM\) | No | B |
| \(n_1\) | Hill Coefficient | \(0.7-3.5\) | \(2\) | \(Dimensionless\) | No | A |
| \(\beta_2\) | Maximal Transcription rate | \(1\) | \(1\) | \(nMmin^{-1}\) | Yes | A |
| \(K_2\) | Dissociation constant | \(10\) | \(10\) | \(nM\) | No | A |
| \(n_2\) | Hill Coefficient | \(0.7-3.5\) | \(2\) | \(Dimensionless\) | No | A |
| \(\alpha_1\) | Degradation rate of sRNA | \(0.03\) | \(0.03\) | \(min^{-1}\) | No | A |
| \(\alpha_2\) | Degradation rate of mRNA | \(0.14\) | \(0.14\) | \(min^{-1}\) | No | A |
| \(\alpha_3\) | Degradation rate of intracellular IL10 | \(0.03\) | \(0.03\) | \(min^{-1}\) | No | A |
| \(\alpha_4\) | Degradation rate of extracellular IL10 | \(0.03\) | \(0.03\) | \(min^{-1}\) | No | A |
| \(k_1\) | sRNA-mRNA binding rate | \(100\) | \(100\) | \(nM^{-1}min^{-1}\) | Yes | C |
| \(k_2\) | Translation rate of IL10 | \(0.3\) | \(0.3\) | \(min^{-1}\) | No | A |
| \(k_3\) | secretion rate of IL10 | \(0.13\) | \(0.13\) | \(min^{-1}\) | No | D |
| \([NO]\) | Nominal concentration of NO | \(13.24\) | \(13.24\) | \(\mu M\) | No | E |
| \([NO]\) | Elevated concentration of NO | \(19.88\) | \(19.88\) | \(\mu M\) | No | E |
| \([Adenine]\) | Nominal concentration of Adenine | \(18\) | \(18\) | \(\mu M\) | Yes | F |
| \([Adenine]\) | Elevated concentration of Adenine | \(100\) | \(100\) | \(\mu M\) | Yes | F |
Reactions were modelled based on ODEs which were solved numerically using ODE15s with an absolute tollerance of \(10^{-30}\) and a relative tolerance of \(10^{-7}\). The differential equations and their response is summarised in the following section.
Dynamics of the system
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 in Figure 1 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 2.
Adenine response dynamics
The complexity of the positive pathway can be seen in Figure 1. 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 2.
As it could be seen from Figure - 2, \([IL10]\) reaches a concentration of \(14.16 \; \mu M\) after roughly 180 seconds. 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-2, 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.
Steady State Model
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- 5. Contours of the surface are shown in Figures 6.
Sensitivity analysis and Optimisation
Sensitivity and sRNA 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 3.
It is clear from Figure-3 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- 4(red curve).
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
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)$$
Body Dynamics and Combined Model
Body Dynamics
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 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 9.
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 9 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 [label] 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 [label]. 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 the response time of \(10^{4}\) seconds seem to be rather accurate as in both scenarios, steady-state is reached not long after that response time.
Frequency Domain analysis
Negative feedback and Assumptions
As discussed earlier on the page, both positive and negative feedback could potentially correct the concentration of IL10 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 if it were in the body. For this purpose, control theory and 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 are going to be analysed separately as based on superposition, the response of our linearised model would be a sum of its response to inputs separately.
- Secretion has been ignored for simplicity in the transfer function.
Therefore, by using the correct transfer function, both models will fall into the same control design with just different transfer functions. The proposed design for the feedback loop could be seen in Figure [label] .
The aim of the controller is to set the reference signal, namely NO and Adenosine, to nominal level as that would correspond to IL10 having a healthy concentration within the body. Transfer functions for both the controller- both positive and negative pathway- and the body are required for control analysis of system response in a human body.
Transfer Functions
As stated in the assumptions, model has been linearised due to non-linearity raised by the Hill functions and sRNA binding; Hence, Hill functions have been replaced by \(\gamma_i^*\) -described below. The block diagram and signal flow graph for the linearised system are shown in Figure [label] , where transcription of IL10 mRNA, sRNA and translation of IL10 could be seen clearly.
$$ \gamma_i^* = \beta_i \frac{n_iK_i^{n_i}[u_i]^{n_i-1}}{(K_i^{n_i}+u_i^{n_i})^2} $$
Transfer functions for both negative and positive pathway were found using Mason's gain relation (Bolton &Newnes, 1998) and were double checked against(Steel &Papachristodoulou, 2017) . Transfer functions of both positive and negative pathways could be found below.
$$ G_1(s) = \frac{\text{IL10(s)}}{\text{Adenine(s)}} = \frac{k_2\gamma_1^*(s+\alpha_1+k_1[\text{IL10 mRNA}]^*)}{(s+\alpha_3)(s+s_+)(s+s_-)}, $$ $$ G_2(s) = \frac{\text{IL10(s)}}{NO(s)} = \frac{-k_1[\text{IL10 mRNA}]^*k_2\gamma_2^*}{(s+\alpha_3)(s+s_+)(s+s_-)}, $$
$$s_\pm = \frac{1}{2}(\alpha_2 + k_1[\text{sRNA}]^* + \alpha_1 +k_1[\text{IL10 mRNA}]^*)\pm \frac{1}{2}\sqrt{d},$$ $$d = (\alpha_2 + k_1[\text{sRNA}]^*- \alpha_1-k_1[\text{IL10 mRNA}]^*)^2+4k_1^2[\text{IL10 mRNA}]^*[\text{sRNA}]^*$$
Body response transfer function
Laplace transformed could be taken from both sides of the Equation [label] and by ignoring \(c\), being very small and negligible, we would be able to derive the transfer function for the body as shown below.
$$ s[\text{NO,Adenine(s)}] = -a_i[\text{NO,Adenine(s)}]+b_i[\text{IL10(s)}] \implies F_i(s) =\frac{[\text{NO,Adenine(s)}]}{[\text{IL10(s)}]} = \frac{b_i}{s+a_i} $$
Stability
The closed loop transfer function for cascade cascade con-troller design given in Figure [label] would contain importantinformation regarding steady state error and relative stabil-ity. It should be mentioned that similar to most control prob-lems where the controller is designed for specific specifica-tions (Dorf & H.Bishop, 1998), our reaction pathway couldalso be changed by changing promoter strength and lengthof sRNA, making it almost possible to adopt desired struc-tures for the controller. The closed loop transfer function forpositive and negative pathways is as follow.
$$ T_i(s) = \frac{F_i(s)G_i(s)}{1+F_i(s)G_i(s)} $$
Based on Equation [label] , the stability of the system could befound by looking at the location of the poles which are theroots to the characteristic polynomial \(1+Fi(s)Gi(s)=0\). Thecharacteristic equation for Adenine and NO reference signalsis given in Equations [label] and [label].