Team:TU Delft-Leiden/Modeling/Landmine

Landmine Module

An important part of our iGEM project is a promoter sensitive to DNT/TNT, first described by Yagur-Kroll et al. [1]. We will use two of the promoters described in this paper, ybiJ and ybiFB2A1, in our project, see the landmine detection module. Of these promoters, not much is known other than the fact that they have a DNT/TNT-dependent response curve (see figure 1). Our goal was to find a model which would be able to reproduce the response curves of both promoters.

To achieve this, we constructed two different models, both using deterministic modeling methods. One model is based on a simple binding model of DNT to the promoter, the other is based on cooperative binding of DNT to the promoter.
When based on the simple binding model, fits of promoter activation with respect to DNT concentration to the experimental data of [1] did not yield good results. However, when the fits were based on the cooperative binding model, we were able to match the experimental data in [1] really well, see figure 2.
From the fits shown in figure 3, it is clear that both the standard and cooperative activation model fail to describe the experimental data we obtained. This might have several reasons. Firstly, the response ratio of our measurements is very small compared to the measurements in [1]. Secondly, our experimental data suggest a detrimental effect of induction with high DNT concentration. This effect is not described in [1] and we don't know its reason. It is therefore not included in our model. Thirdly, both the data set from [1] and our experimental data set are small compared to the amount of parameters. This makes finding the right model difficult.
To improve the promoter activation model, a larger and more consistent data set has to be obtained. We need to conduct fluorescence measurements over a wide range of DNT concentrations. Besides that, possible detrimental effects of its solvent (acytonitrile) have to be investigated. Also, the induction of the N genes with rhamnose might influence the measurements in an unexpected way. Doing all this is unfortunately not possible in the time span of this iGEM project.

Simple Binding Model

Our first approach was to solve a system of Ordinary Differential Equations (ODEs) resembling the transcription and translation of a gene activated by the DNT-sensitive promoter. The ODEs were derived from the following system of reactions: $$ P_{R} + \ DNT \ \overset{k_{on}}{\underset{k_{off}}{\rightleftharpoons}} \ P_{A} \tag{1}$$ $$ P_{A} \ \xrightarrow{s_{A}} \ P_{A} + \ mRNA \tag{2}$$ $$ P_{R} \ \xrightarrow{s_{R}} \ P_{R} + \ mRNA \tag{3} $$ $$ mRNA \ \xrightarrow{s_{P}} \ mRNA + \ R \tag{4} $$ $$ mRNA \ \xrightarrow{d_{M}} \ \emptyset \tag{5} $$ $$ R\ \xrightarrow{d_{P}} \ \emptyset \tag{6} $$ In these reaction equations, \(P_{R}\) and \(P_{A}\) indicate respectively repressed and active promoters. DNT indicates DNT molecules, mRNA stands for an mRNA molecule transcribed from the gene behind the promoter and R is the reporter protein which is translated from the mRNA. \(k_{on}\) and \(k_{off}\) are the rates at which a promoter goes from the repressed to the active state and vice versa. \(s_{A}\) is the transcription rate from an active promoter. \(s_{R}\) is the transcription rate from a repressed promoter; this is also referred to as leakage. \(s_{P}\) is the translation rate. \(d_{M}\) and \(d_{P}\) are the mRNA and protein degradation rates.

Reaction (1) describes the activation and deactivation of a promoter in the presence of DNT. We assume that the activation mechanism can be described as the binding of the DNT to the repressed promoter, with the resulting complex being an active promoter. Reaction (2) describes transcription. Reaction (3) described transcription through leakage, i.e. transcription from a repressed promoter. Reaction (4) describes translation. Reactions (5) and (6) describe mRNA and reporter protein degradation respectively.

This system of reactions leads to the following system of ODEs: $$ \frac{d}{dt} [P_{R}] = \ -k_{on}[P_{R}][DNT] \ + \ k_{off}[P_{A}] \tag{7.1} $$ $$ \frac{d}{dt} [P_{A}] = \ k_{on}[P_{R}][DNT] \ - \ k_{off}[P_{A}] \tag{7.2} $$ $$ \frac{d}{dt} [DNT] = \ -k_{on}[P_{R}][DNT] \ + \ k_{off}[P_{A}] \tag{7.3} $$ $$ \frac{d}{dt} [mRNA] = \ s_{A}[P_{A}] \ + \ s_{R}[P_{R}] \ - \ d_{M}[mRNA] \tag{7.4} $$ $$ \frac{d}{dt} [R] = \ s_{P}[mRNA] \ - \ d_{P}[R] \tag{7.5} $$ A system of ordinary differential equations like this can be solved numerically, with the MATLAB function ode45. The result of this will be a set of curves describing the concentration of each compound in time. However, when choosing parameters with an realistic order of magnitude, the time needed to reach steady state, i.e. when the reporter protein concentration does not change anymore, is so long that the script turns extremely computationally intensive. Because time-dependence is not relevant when imitating the dose-response curves from the Yagur-Kroll paper [1], we decided to discard the ODE method and switch to analytic steady state methods.

The goal of this steady state approach is to find a relation between the concentration of reporter protein ([R]) and the amount of DNT ([DNT]), with as few parameters as possible. To obtain such a relation, we started from equation (7.2), which describes the change in time of the concentration of activated promoters. We used a steady state assumption to simplify this relationship, which means that we assumed that the concentration of active promoters does not change in time, i.e \(\frac{d}{dt} [P_{A}] = \ 0\). This assumption can be justified by examining the experimental data in figure 4, from which it can be clearly seen that after a certain amount of time, the reporter protein level reaches a constant value.

This yields the following: $$ [P_{A}] = \ \frac{k_{on}}{k_{off}}[P_{R}][DNT] = \ \frac{1}{K_{d}}[P_{R}][DNT] \tag{8} $$ Here we introduced the dissociation constant \(K_{d} = \ \frac{k_{on}}{k_{off}}\). The next assumption we make is that the total amount of promoter (\([P_{T}\)) is constant. This makes sense, since an active promoter can become a repressed promoter and vice versa, but it is not possible to make new promoters or let promoters vanish. We ignore and therefore neglect the fact that during cell replication, the promoter level will rise due to DNA replication. However, we assume that this is a minor contribution. This assumption yields the following: $$ [P_{T}] = \ [P_{A}] + \ [P_{R}] \ \rightarrow \ [P_{R}] = \ [P_{T}] - \ [P_{A}] \tag{9} $$ Combining equations (8) and (9), we obtain: $$ [P_{A}] = \ \frac{1}{K_{d}}([P_{T}] - \ [P_{A}])[DNT] \tag{10} $$ This can be rewritten to the following form: $$ \frac{[P_{A}]}{[P_{T}]} = \ \frac{[DNT]}{K_{D} + \ [DNT]} \tag{11} $$ Equation (11) thus gives us the number of active promoters as a fraction of the total amount of promoters.

The change of concentration of reporter protein ([R]) is given by the following ODE: $$ \frac{d}{dt} [R] = \ \frac{[P_{A}]}{[P_{T}]}k_{1} + \ \frac{[P_{R}]}{[P_{T}]}k_{2} - \ d_{P}[R] \tag{12} $$ In equation (12), \(k_{1}\) is the rate of protein production from an active promoter (this could be regarded as a combination of \(s_{A}\) and \(s_{P}\)) and \(k_{2}\) is the rate of protein production from a repressed promoter (combination of \(s_{R}\)and \(s_{P}\)). In this equation, the transcription/translation-process is wrapped up in one equation instead of two as in equation (7.4) and (7.5). This can be justified by the fact that in in vivo systems, transcription and translation are strongly coupled [3]. Plugging in equation (11) and again considering the steady state, we arrive at the following expression for the equilibrium concentration of the reporter protein: $$ [R]_{eq} = \ \frac{1}{d_{P}} \frac{k_{1}[DNT]+ \ K_{D} k_{2}}{K_{D} + \ [DNT]} \tag{13} $$ For the protein degradation rate we used the fixed value of \(4 \times 10^{-6} \ s^{-1}\), which is estimated from the data in [2]. To obtain the values for \(K_{D}\), \(k_{1}\) and \(k_{2}\), we fitted equation (13) to the data from the Yagur-Kroll paper [1]. To do this, we used the nlinfit function in MATLAB.

From figure 5, it is clear that the fits do not match the data. This is a strong indication that the model as described by equation (13) is not correct for the system at hand. We therefore went on to explore another model.

Cooperative Binding Model

Since one-on-one interaction between DNT and the promoter did not yield good results, we hypothesized that this might be a case of cooperative interaction. The Hill coefficient (n) is commonly used as a way to quantify the cooperative binding effect. Making use of the Hill coefficient, equation (7.2) can be rewritten for the case of cooperative binding: $$ \frac{d}{dt} [P_{A}] = \ k_{on}[P_{R}][DNT]^{n} - \ k_{off}[P_{A}] \tag{14} $$ This equation is almost completely similar to equation (7.2), except for the Hill coefficient appearing as an exponent. The analysis needed to find a equation relating the equilibrium concentration of the reporter protein to the amount of DNT is also completely similar to the analysis done to obtain equation (13). Therefore, we will suffice with giving the end result: $$ [R]_{eq} = \ \frac{1}{d_{P}} \frac{k_{1}[DNT]^{n}+ \ K_{D}^{n} k_{2}}{K_{D}^{n} + \ [DNT]^{n}} \tag{15} $$ Fitting equation (15) (again with the fixed value of \(4 \times 10^{-6} \ s^{-1}\) for protein degradation) to the data in the Yagur-Kroll paper gives the parameter values given in table (1) and the fits shown in figure 6. Common methods to determine confidence intervals for the fit parameters such as obtaining them from the covariance matrix or the Jacobian did not yield correct results. We therefore tried to determine confidence intervals for the promoter parameters with a bootstrapping procedure adapted from [4]. However, we were not able to determine confidence intervals due to non-convergence of some fits. This is a clear indication that the uncertainty of these values is quite high.

From figure 6, it can be seen that the model described by equation (15) fits remarkably well to the experimental data. It is therefore worthwhile to look into the found parameter values a bit more. The lower dissociation constant for the yqjFB2A1 promoter compared to the jbiJ promoter indicates that the DNT binds better to the yqjFB2A1 promoter. This is not surprising, since the yqjFB2A1 promoter is the result of a directed evolution experiment in which is DNT response was enhanced with respect to the original yqjF promoter, which had similar characteristics as the jbiJ promoter [1]. This is also the most obvious explanation of the lower values for the Hill coefficient. A large increase in both \(k_{1}\) and \(k_{2}\) can be seen between the jbiJ promoter and the yqjFB2A1 promoter. This indicates that the reporter protein production does not only increase for the activated promoter, which is as expected, but also for the repressed promoter. This corresponds to the fact that in [1], for the yqjFB2A1 promoter, a higher background luminescence was recorded than for the non-mutated promoters, see figure 7.

Since we plan to use the found values of \(k_{1}\) in other parts of the modelling, we performed a quick back-of-the-envelope calculation to check if the order of magnitude of \(k_{1}\) is realistic. In [1], an \(OD_{600}\) of 0.2 is recorded. This corresponds to a cell concentration of \(1.6 \times 10^{8} \ mL^{-1}\) or approximately \(2.7 \times 10^{-13} \ M\). With a value for \(k_{1}\) of \(7.3 \times 10^{-12} \ M\), this corresponds to a protein production of approximately 27 proteins per cell per second, which seems fairly reasonable.

Although the model described by equation (15) is able to match the experimental data in [1] really well, assuming that the cooperative binding of DNT to the promoter is the actual mechanism that happens in nature, is too short-sighted. It was already hypothesized by Jagur-Kroll et al. [1] that "the induction of the yqjF and ybiJ promoters is probably not caused by our target compounds but rather by their metabolites or degradation products, possibly a quinol or quinol derivative." To obtain a more suitable model description of the landmine promoters, further research in their activation mechanism is needed. However, this is beyond the scope of iGEM and therefore we will suffice with the cooperative binding model described by equation (15).

Fitting to Experimental Data

To further test our promoter activation models, we fitted both models to experimental data found with a plate reader experiment. We decided to fit only the data sets of the yqjF and ybiJ promoters combined with the N genes, since these are the only ones to show a clear increase in fluorescence. The results are shown in figure 8.

From the fits shown in figure 8, it is clear that both the standard and cooperative activation model fail to describe the experimental data we obtained. This might have several reasons. Firstly, the response ratio of our measurements is very small compared to the measurements in [1]. Secondly, our experimental data suggest a detrimental effect of induction with high DNT concentration. This effect is not described in [1] and we don't know its reason. It is therefore not included in our model. Thirdly, both the data set from [1] and our experimental data set are small compared to the amount of parameters. This makes finding the right model difficult.

To improve the promoter activation model, a larger and more consistent data set has to be obtained. We need to conduct fluorescence measurements over a wide range of DNT concentrations. Besides that, possible detrimental effects of its solvent (acytonitrile) have to be investigated. Also, the induction of the N genes with rhamnose might influence the measurements in an unexpected way. Doing all this is unfortunately not possible in the time span of this iGEM project.

References

[1] S. Yagur-Kroll, S. Belkin et al., “Escherichia Coli bioreporters for the detection of 2,4-dinitrotoluene and 2,4,6-trinitrotoluene”, Appl. Microbiol. Biotechnol. 98, 885-895, 2014.

[2] M.R. Maurizi, “Proteases and protein degradation in Escherichia coli”, Experientia 48(2), 178-201. p.181, 1992, table 2.

[3] M. Iskakova, K. Nierhaus et al., “Troubleshooting coupled in vitro transcription-translation system derived from Escherichia coli cells: synthesis of high-yield fully active proteins”, Nucleic Acid Research 34(19), 2006.

[4] D. Das, “Experimental data analysis Lecture 3: Confidence intervals”, Mathematical Cell Biology Graduate Summer Course Lecture conducted from University of British Columbia, Vancouver, May 2012.
Available at: http://www.math.ubc.ca/~keshet/MCB2012/SlidesDodo/DataFitLect3.pdf [Accessed 27 Sept. 2014].