Device specificities and optimization

Optimization of the device

As the modelling part shows our models for Quorum Sensing allowed us to have a visual representation of our entire device. This modelling highlighted the diffusion of the quorum sensing molecule.

Figure 1: Animation generated through MATLAB for visual representation of our models and the complete deterministic simulation

Decreasing the accuracy of the measure, quorum sensing diffusion is a problem for our system. To fixe this problem, we decided to developed a device with channel containing each of them a different IPTG concentration but the same aTc concentration.

Figure 2: First prototype of the device

This kind of device allows us to limit the diffusion of quorum sensing at one channel. Of course the measure won't be a unique concentration, but a range of concentration. However, the problem posed by quorum sensing diffusion would disrupt the measure more than the channel device (if there are a sufficent number of channel and a small enough IPTG step).

Determination of the limit of quantification

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 3: 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 4: Hysteresis for $[aTc] = 1x10^{-9} M$

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

Figure 6: 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 7: 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 8: 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 9: 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.