Analysis of the system

Simulation of the bacteria distribution on the plate

The toggle switch behavior

At early stage, the goal of the modelling team was to confirm the behaviour of the whole circuit. We divided the the network into two main models, Toggle switch and Quorum Sensing. Very early the modelling results seemed promising and we could rapidly infer that our Toggle Switch design would be effective. Indeed, with the models described in chapter 3, we can see the behaviour of our bacteria on the plate. On the plate, one whole region features bacteria in the LacI way and the rest of the plate features bacteria in the TetR way.

The following simulation was realized for an IPTG gradient of $1x10^{-6} M$ to $1x10^{-2} M$ and an aTc concentration of $5x10^{-6} M$. The first graph present the logarithmic IPTG gradient in green and the homogeneous concentration of aTc in red. The second represent the concentration of both repressor on the plate.

Figure 1: Observation of the switch on the plate for an aTc concentration of $5x10^{-6}$

The first thing we could observe on this figure is that the switch doesn't appears at the equality of the concentration of aTc and IPTG but for an IPTG concentration of $1.5x10^{-4} M$. This is due to the value of the parameters in the ODE system presented previously. In fact, the dissociation constant of respective repressor and their inhibitor are not the same.

Figure 2: Observation of the switch on the plate for a higher aTc concentration of $5x10{-5}$

On the previous two figures, X axis represents physical points on the plate, form left to right of the plate. In each of these points the only difference is the IPTG concentration, as we will apply on our plate an IPTG gradient. The interface between the two regions depends on [aTc]. Higher aTc concentration will move the interface to the right edge of the plate as in figure. We therefore demonstrated that the Toggle switch behaviour was the one we wanted for our application.

With this model, we also demonstrated that the presence of degradation tags were necessary to get the appropriate behaviour. If the degradation rate of the LacI and TetR proteins were too long (typical half-time of 10 hours) the concentrations in each protein would be too high and the switching in one way or another would be way too long for our application. As a result we decided to use only LVA tagged LacI and TetR genes which impose their half-life time of 10 minutes.

Demonstration that the Toggle switch behaviour was the one we wanted for our application.
Use only LVA tagged LacI and TetR genes which impose the half-life time of 10 minutes.

Quorum Sensing

Our models for Quorum Sensing allowed us to simulate the behaviour of our whole system, confirm our expectations and finally have a visual representation of our entire device.

In a first step, we observed the distribution of the protein acting in the quorum sensing system and the concentration of internal and external quorum sensing. The objective is to show that coupling toggle switch and quorum sensing modelling works.

Figure 3: Observation of the Quorum sensing molecule distribution on the plate

On the first graph of this figure we see(in green) the cinI concentration (which follows the same equation as lacI) and (in red) the cinR concentration (which follows approximately the same equation as tetR). If cinR concentration is not as high as cinI concentration, it's because in cinR equation we needed to take into account the complexation of cinR with the quorum sensing molecule as a disparition term.

Moreover, the two other curve in the first figure show the concentration of the quorum sensing molecule inside and outside the cells. And we see that, because of the diffusion of the quorum sensing molecule in the medium (third graph), the internal concentration of quorum sensing is not equal to zero where cinI is absent. Which indicate that quorum sensing well diffused in the medium and was caught by receiving bacterias.

In the following graphs we show the complexation of cinR with the quorum sensing molecule.

Figure 4: Observation of the Quorum sensing complexation with cinR receptor

On the first graph of this figure, intern quorum sensing concentration (in green) and cinR concentration (in red) are plotted. We can well see that there is an area on the plate where cinR concentration and intern quorum sensing concentration are both not equal to zero. This is predicting that a complexation between both of them could happen.
That's what it's shown is the second graph of this figure, the concentration of the cinR/Quorum Sensing complex in the bacterias.

With the two previous figures, we can confirm that the quorum sensing is diffusing on the right side of the plate. This quorum sensing should be caught by the receiving bacteria. This would produce lycopene and activate a diffused coloration on the plate.

Figure 5: Observation of the red stripe on the plate

With modelling we show that the system should work as expected. But we also hightlighted a problem: the diffusion of the quorum sensing which is decreasing the accuracy of the measure. To fixe this problem, we needed to optimize our device

Stability Studies of the Toggle Switch

Nullclines studies

In order to predict the set point and the specifications of our system, we studied first the existence and the value of the steady state solutions of the set of ODE.

Isocline study is a classical study which implies a research of stationnary point in a system. These stationnary points are deduced from the equations of the differential system: when the variation of concentration of both repressors are equal to zero.

$\frac{d[TetR]}{dt} = \frac{k_{pLac}.[pLac]_{tot}}{1 + (\frac{[lacI_{total}]}{K_{pLac} + \frac{K_{pLac}.[IPTG]}{K_{lacI-IPTG}}.})^\beta} - \delta_{TetR}.[TetR] = 0$
$\frac{d[lacI]}{dt} = \frac{k_{pTet}.[pTet]_{tot}}{1 + (\frac{[tetR_{total}]}{K_{pTet} + \frac{K_{pTet}.[aTc]}{K_{TetR-aTc}}.})^\gamma} - \delta_{lacI}.[lacI] = 0$

To facilitate the manipulation of the equation and reduced the number of parameters, we posed:

• $E_{TetR}$ = $k_{pLac}.[pLac]_{tot}$
• $R_{TetR}$ = $\frac{1}{1 + (\frac{[TetR_{total}]}{K_{pMerT} + \frac{K_{pTet}.[aTc]}{K_{TetR-aTc}}.})^\gamma}$
• $E_{LacI}$ = $k_{pTet}.[pTet]_{tot}$
• $R_{LacI}$ = $\frac{1}{1 + (\frac{[LacI_{total}]}{K_{pLac} + \frac{K_{pLac}.[IPTG]}{K_{LacI-IPTG}}.})^\beta}$
• $[TetR]_{r}$ = $R_{TetR}.[TetR]$ the relative concentration of TetR
• $[LacI]_{r}$ = $R_{LacI}.[LacI]$ the relative concentration of LacI
• $K$ = $\frac{R_{TetR}.E_{TetR}}{\delta_{TetR}}$
• $K_{prime}$ = $\frac{R_{LacI}.E_{LacI}}{\delta_{LacI}}$

After manipulation with these reduced parameters, we get this following equations:

$[TetR]_{r} = \frac{K}{1 + ([LacI]_{r})^\beta}$ (1)
$[LacI]_{r} = \frac{K_{prime}}{1 + ([TetR]_{r})^\gamma}$ (2)

From this equation we could see that, if $[LacI]_r$ >> 1, $[TetR]_r = 0$ and $[LacI]_r \approx K_{prime}$. In the other case if $[TetR]_r$ >> 1, $[LacI]_r = 0$ and $[TetR]_r \approx K$

From these equations, we get this figures:

Figure 6: Solution of the equation and emergence of three steady state

On this figure, the red lines represent the solution of the equation (1) and the green line the solution of (2). This figure was realized with $[aTc] = 5x10^{-6} M$ and $[IPTG] = 1.55x10^{-4} M$. These parameters reflect the situation of our system in the center of the plate in the presence of a logarithmic gradient of IPTG of $1x10^{-6} M$ to $1x10^{-2} M$.

Three stationary points emerge from this graph. These are the three points of intersection of two curves and represent the steady state of the system.
However, there is one of the three points which is an unstable steady state: the point 2. It represents the point when both relative concentration are equal. In a Toggle Switch, it's impossible to have concentration of both repressors equal because one repressed the other. So one of these should take the avantage on the other.

Figure 7: Nullclines for the left side of the plate

Figure 8: Nullclines for the right side of the plate

These figures were realized with $[aTc] = 5x10^{-6} M$ and for the left curve with $[IPTG] = 1x10^{-6} M$ and for the right curve with $[IPTG] = 1x10^{1} M$. The left graph represents the left side of the plate where aTc concentration is dominant and the right graph represents the right side of the plate where IPTG concentration is dominant.
These figures show that when the concentration of one of the repressor is too high, the system is no longer bistable but monostable.

Stochastic analysis of the stability

By working on histograms, we get the distribution of bacteria's states(lacI or tetR pathway) along the plate.

Figure 9:Histogram for several runs on the same point of the plate. We are far from interface and only the LacI way is transcripted. X axis is normalized concentrations and the Y axis is number of runs that finished with the corresponding concentration (negative for LacI and positive for tetR)

This figure show the bacteria distribution in the left of the plate, where aTc is predominent. The green peak indicates bacterias in the lacI pathway. Which is showing to us that in the left of the plate, bacteria could only be in the lacI genetic pathway. The distribution is monomodal.

Figure 10:Histogram for several runs on the same point of the plate. It is on one point of the interface between LacI area and TetR area of the plate. (LacI = green; TetR = blue)

This figure show the bacteria distribution at the interface. The presence of two peaks indicates that bacterias are presents both in the lacI pathway and the tetR pathway as we were expecting. At this point the two ways are equally likely to be chosen in the cell, which is why we have an interface.

As we saw it with the nullcline study, stochastic modelling shows that on the edge of the plate, the toggle switch is monostable and at the interface it's bistable.

Conclusion about stability

According to the previous studies, we were able to predict(in fonction of aTc and IPTG concentration) where the system is monostable and where it's bistable.

Figure 11: Stability of the toggle switch on the plate

• On the extreme side of the plate, the system is monostable.
• On the switch area of the plate, the system is bistable.
• Bistability, in the switch area, allows us to obtain neighboring bacteria in different states. These bacterias could communicate together and give rise to the coloration

Device Specificities

Determination of the inferior quantification limit by hysteresis study

The goal of the hysteresis study is to examine the switch conditions when the toggle switch is already locked in a predefined pathway. In our mathematical study, we blocked the system in the lacI pathway with different preliminary aTc concentrations. Then the amount of IPTG was increased until the system switched. Then, we decreased IPTG concentration to see when the system switched back to the initial state. The blue curve shows the evolution of TetR concentration when IPTG concentration grows. The red curve shows the evolution of TetR concentration when IPTG concentration decreases.

Figure 12: Hysteresis curve for $[aTc] = 1x10^{-6} M$

To quantitatively exploit these curves we determined at which IPTG concentration the system switched. On this curve, we can get the switch up concentration: ~ $1x10^{-2} M$ of IPTG and the switch back concentration ~ $3x10^{-5} M$.

The switch back concentration is very similar to the dissociation constant between lacI and IPTG (which is $2.96x10^{-5} M$). It means that, when there is not enough IPTG in the bacteria, the IPTG-lacI complexe is faster degraded than produced. So the repression is no longer effective.

The following curves show different hysteresis for growing aTc concentrations:

Figure 13: Hysteresis for $[aTc] = 1x10^{-9} M$

Figure 14: Hysteresis for $[aTc] = 1x10^{-8} M$

Figure 15: Hysteresis curve for $[aTc] = 1x10^{-7} M$

The two first curves (for $[aTc] = 1x10^{-9} M$ and $[aTc] = 1x10^{-8} M$) show that the switch up and switch back concentrations are the same. This concentration is ~ $3x10^{-5} M$, the dissociation constant between lacI and IPTG. The last curve (for $[aTc] = 1x10^{-7} M$) shows that the switch back concentration stay the same. But the switch up concentration is higher. In fact, for aTc concentration superior to $1x10^{-7} M$, the switch up concentration is growing with aTc concentration.

The concentration of aTc $1x10^{-7} M$ appears to be the limit of sensibility to the toggle switch. This concentration represents the dissociation constant of aTc to TetR repressor.

Hysteresis permits to determined the inferior limit of quantification of our device: $1x10^{-7} M$ of aTc, which is the dissociation constant between aTc and TetR.

Mathematical calibration of the system

With the model we developed, you have to run a simulation if you want to know for a certain IPTG concentration what is the aTc concentration making the system switched. Simulation could take a lot of time and it's not the best way to calculate the aTc concentration.
So, we tried to developed an equation capable of giving the same result as the simulation but much faster.

The best solution would have been used the ODE system of the toggle switch and apply limited development on this system to get the equation we were looking for. But when the 2 equations are coupled it's more difficult. So we tried an other way to get this equation.

We also used the ODE system of the toggle system. From this system, we made the following hypothesis:

1. We are in a steady state: both of the equation are equal to zero
2. The synthesis rate of both promoters and the degradation terms of both repressors are the same
3. The cooperativity of repression number ($\beta$ and $\gamma$) are approximately equal and in the steady state [lacI] and [tetR] are the same

Using these hypothesis we get the following development:

$ \frac{d[TetR]}{dt} = \frac{k_{pLac}.[pLac]_{tot}}{1 + (\frac{[lacI]}{K_{pLac} + \frac{K_{pLac}.[IPTG]}{K_{lacI-IPTG}}.})^\beta} - \delta_{TetR}.[TetR] = 0$
$\frac{d[lacI]}{dt} = \frac{k_{pTet}.[pTet]_{tot}}{1 + (\frac{[tetR]}{K_{pTet} + \frac{K_{pTet}.[aTc]}{K_{TetR-aTc}}.})^\gamma} - \delta_{lacI}.[lacI] = 0$

By applying the first hypothesis we get:

$ \frac{k_{pTet}.[pTet]_{tot}}{1 + (\frac{[tetR]}{K_{pTet} + \frac{K_{pTet}.[aTc]}{K_{TetR-aTc}}.})^\gamma} - \delta_{lacI}.[lacI] = \frac{k_{pLac}.[pLac]_{tot}}{1 + (\frac{[lacI]}{K_{pLac} + \frac{K_{pLac}.[IPTG]}{K_{lacI-IPTG}}.})^\beta} - \delta_{TetR}.[TetR]$

By applying the second hypothesis we get:

$ \frac{1}{1 + (\frac{[tetR]}{K_{pTet} + \frac{K_{pTet}.[aTc]}{K_{TetR-aTc}}.})^\gamma} = \frac{1}{1 + (\frac{[lacI]}{K_{pLac} + \frac{K_{pLac}.[IPTG]}{K_{lacI-IPTG}}.})^\beta}$

Which become:

$ (\frac{[tetR]}{K_{pTet} + \frac{K_{pTet}.[aTc]}{K_{TetR-aTc}}.})^\gamma =(\frac{[lacI]}{K_{pLac} + \frac{K_{pLac}.[IPTG]}{K_{lacI-IPTG}}.})^\beta$

And finally by applying the third hypothesis, we get the following equation

$ [aTc] = K_{tetR-aTc}(\frac{K_{pLac}}{K_{pTet}} (1 + \frac{[IPTG]}{K_{lacI-IPTG}})$

We finally get an equation depending on parameters of the system and the IPTG concentration.
To test this equation, we compare aTc concentration obtained with simulation and with the equation.

[IPTG] M	[aTc] obtained by simulation M	[aTc] obtained by equation M	Deviation
0,0000255	0,0000005	2,89E-07	42,13%
8,90E-05	1,00E-06	6,40E-07	35,99%
1,55E-04	0,0000015	1,00E-06	33,02%
2,52E-04	2,00E-006	1,54E-06	22,98%
3,18E-04	2,50E-006	1,91E-06	23,80%
5,54E-04	4,00E-006	3,21E-06	19,78%
7,14E-04	5,00E-006	4,09E-06	18,15%
1,60E-03	1,00E-005	8,99E-06	10,14%
3,40E-03	2,00E-005	1,89E-05	5,36%
8,30E-03	5,00E-005	4,60E-05	8,01%
1,18E-02	7,00E-005	6,53E-05	6,67%
0,0172	1,00E-004	9,52E-05	4,84%
0,085	5,00E-004	4,70E-04	6,07%
0,1723	1,00E-003	9,52E-04	4,81%

Considering a deviation inferior to 10% acceptable. The equation demonstrated previously is applicable only for IPTG concentration higher than $1.6x10^{-3} M$ which correspond to an aTc concentration of $1x10^{-5} M$. For lower concentration, the equation wouldn't give an good estimation of the concentration, but it could give a good range of where the switch will appear.

Stochastic study for statistic determination of severals device specificities

Even though deterministic modelling predicted a promising behaviour for our system, we modelled our system with Stochastic algorithms in order to check the robustness of our predictions with a highly stochastic medium and to get statistical information on our system.
For biosensors the importance of stochastic modelling is clear, it gives a lot of information on the precision of the measure that is mainly caused by the inner randomness of the genetical network.

Stochastic simulation has been performed by many iGEM teams during the previous competitions. However, many of the results obtained by those previous teams were merely analysed quantitatively. The amount of information obtained via Gillespie simulation is therefore wasted. We wanted to exploit these results and set Gillespie simulation as an unavoidable modelling aspect of synthetic biology, especially in the case of biosensors.

In order to do this, we performed a statistical analysis of the results obtained. We used the results of this analysis for the sizing of our final device.

We started with an estimation of the mean of the $LacI$ and $TetR$ variables over the entire plate. Computing this simulation required 5 computers running for about 70 hours. On each of 200 points of the plate, we computed 1000 runs.

The estimators for $LacI_{cell}$, $TetR_{cell}$ and $LacI \times TetR_{cell}$ variables are simple, non-biased estimators:

$\displaystyle\hat{\mu}_{LacI} = \frac{1}{n}\sum_{i=1}^{n}LacI_{i}$
$\displaystyle\hat{\sigma}_{LacI}^{2} = \frac{1}{n-1}\sum_{i=1}^{n}(LacI_{i} - \overline{LacI})^{2}$

Of course similar estimators are used for $TetR$ and $LacI \times TetR$.

Figure 16: Curves of the means of TetR (red) and LacI (blue) variables over a normalised IPTG gradient

We can see on this figure that the interface between the LacI area of the plate and the TetR area is important. Its width will of course depend on the setting of the IPTG gradient ($\Delta IPTG$) between the channels on the plate. We draw the curve corresponding to $E(LacI \times TetR)$ :

Figure 17: Curve of the mean of $LacI \times TetR$(green)

The width of the bell-shaped green curve will set the $\Delta IPTG$ between each well. One would understand that a proper IPTG step will be smaller than the width of this curve. If it were bigger than this width there would be a chance that no channel turns red.

On the other hand, if the IPTG step between channel was too small, there would be too many channels turning red without any way to know which one is the actual center of the interface.

Here we decided to set the smallest $\Delta IPTG$ so that a maximum of 3 channels could possibly be in the top $10\%$ of the bell-shaped $E(LacI \times TetR)$ curve. We get the IPTG resolution range for this particular point: between $6x10^{-3} M$ and $3.4x10^{-6} M$. These results are of course simple estimations and need more experimental validations in order to get a precise knowledge of the levels of coloration, for example.

Another important aspect of the sensor was its Standard Error of Measure. We needed to know the variance of $LacI \times TetR$.

Figure 18: Estimated standard deviation of the $LacI \times TetR$ variable

$Var(LacI \times TetR) $ was computed on all different points on the plate and final result was not surprisingly higher at the interface. The maximal estimated value of standard deviation in the interface is here $1.9x10^{4}$ $(proteins^{2}/cell)^{2}$.

Knowing the number of bacteria we will put in the channels of the plate, we will then know the precision of the sensor. The $LacI \times TetR$ variable is a sum of all cells' proteins over the channel. Which makes it a sum of independant random variables ( $LacI_i \times TetR_i $ and $LacI_j \times TetR_j $ independant for $i \neq j$).

The precision can therefore be calculated with the Central Limit theorem :

$precision_{68\%} = \frac{\sigma_{LacI \times TetR_{cell}}}{\mu_{LacI \times TetR_{cell}}\sqrt{n_{cells}}}$

For example, if we want to be sure that $68\%$ of the bacteria will be within +/- $10\%$ around the mean of the channel, we can therefore state that 9025 bacteria are needed in the channels at least. Knowing that the number of bacteria per channel will be about millions, we can be sure that the precision will be much higher.

Justification of the need of a regulation system

In our project, we developed a translation regulation system in order to keep the system OFF until the pollutant is added. With modelling, we show that this system is really important and without it, the measure is really disturbed.

In all of our simulation, we considered the initial concentrations of both repressors equal to zero. This will be the case with the regulation system. Without this system, the initial concentrations of the repressors are higher. In the following figure, we performed a simulation for an aTc concentration of $1.10^{-6}$ with initial concentrations equal to zero and one with initial concentrations equal to $5\%$ of the concentrations in the steady state of the previous simulation.

Figure 19: Observation of the interface when the regulation system works and when it doesn't for an aTc concentration of $1x10^{-6}$.

When the regulation system doesn't work, the interface is not at the place it supposed to be. Because we can't have measures of the initial concentration of both repressors, to well predict where the interface will appear, we need to control these concentrations. That's why we developed the translation regulation system.