Our equations and how we obtained them

The Two levels of Modelling

In order to facilitate the construction of the model, we divided the model into two levels :

• Toggle Switch modelling: to predict the state bacteria will be.
• Quorum Sensing and Coloration: to predict where the red line will appear on our device.

The first level of the model allows us to define the switch conditions in bacteria, depending on pollutant and IPTG concentration. From these results, we can see the switch zone on the device and so the red line position. Thanks to these simulations, we can resize our device for a better detection efficiency.

The second level gives us the width of the red line and so the precision of the system. From these results, we could improve the precision of our device. For example it allows us to choose between a device with a continuous baterial layer and a device with strips containing bacteria. Strips have different concentration of IPTG and the same concentration of Mercury.

However, in our model, we didn't consider Mercury but anhydro-tetracycline (aTc). It's a set of well characterized genes in E. coli and it's really useful for the model because it allows us to do a proof of concept for our system. If the system works with aTc, it could also work with mercury, copper or lead for example.

So in the following chapter, we are going to use TetR instead of MerR and aTc instead of Mercury.

Toggle switch

1. Biological Models

A toggle switch consists of two genes, each coding for a protein that represses the expression of the other gene. This double repression system ensures the basic function of a toggle switch: a bistable system which can be switched from one state to the other by putting, in our example, some IPTG or aTc molecules in the medium.

The biological system we are trying to implement is more complex, on both the biological and physical side. However, the toggle switch model is basically the same. In the study of the toggle switch itself, the system can be reduced to a simple two-state subsystem that we will then use for the rest of our modelling. The toggle switch itself is not influenced by the rest of the system, if we do not consider the RsmA regulatory system that will be modeled at the very end of our work. Thus we will be able to model the toggle switch independently and then build the rest of the model upon this basis.

2. Mathematical Models

We use a common model for the toggle switch that we demonstrate below. The differential equation describing the production of TetR is as follows

$$\frac{d\left[\mathrm{TetR}\right]}{dt}=\frac{k_{\mathrm{plac}}\left[P_{\mathrm{lac\; total}}\right]}{\mathrm{1\; +}{\left(\frac{\left[\mathrm{LacI}\right]}{K_{\mathrm{plac}}+\frac{K_{\mathrm{plac}}\left[\mathrm{IPTG}\right]}{K_{\mathrm{LacI\; -\; IPTG}}}}\right)}^{n_{\mathrm{plac}}}}-{\mathrm{\&;}}_{\mathrm{TetR}}\left[\mathrm{TetR}\right]$$

For better understanding of this model we demonstrated it. The differential equation describing the production of xR is as follow :

$$\frac{d\left[\mathrm{TetR}\right]}{dt}=k_{\mathrm{plac}}\left[P_{\mathrm{lac\; free}}\right]-{\delta }_{\mathrm{TetR}}\left[\mathrm{TetR}\right]$$

with [Plac free ] being the concentration of available binding sites - i.e. not repressed by LacI molecule. Plac_free is of course related to the total number of promoters Plac : $$\left[P_{\mathrm{lac\; free}}\right]+\left[P_{\mathrm{lac\; -\; LacI}}\right]=\left[P_{\mathrm{lac\; total}}\right]$$ with [Plac − LacI ] being the concentration of promoters repressed by LacI. If we set $K_{\mathrm{plac}}=\frac{\left[P_{\mathrm{lac\; free}}\right]\left[\mathrm{LacI}\right]}{\left[P_{\mathrm{lac\; -\; LacI}}\right]}$ we get : $$\left[P_{\mathrm{lac\; free}}\right]=\frac{\left[P_{\mathrm{lac\; total}}\right]}{\mathrm{1\; +}\frac{\left[\mathrm{LacI}\right]}{K_{\mathrm{plac}}}}$$ We then try to get [LacI] : $$\left[\mathrm{LacI\; -\; IPTG}\right]+\left[{\mathrm{LacI}}_{\mathrm{free}}\right]+\left[\mathrm{LacI\; -\; plac}\right]=\left[\mathrm{LacI}\right]$$ with [LacI - IPTG] the concentration of LacI bound to IPTG and [LacI - plac ] the concentration of the complex of LacI and the promoter. If we set $K_{\mathrm{LacI\; -\; IPTG}}=\frac{\left[\mathrm{LacI}\right]\left[\mathrm{IPTG}\right]}{\left[\mathrm{LacI\; -\; IPTG}\right]}$ we get : $$\left[\mathrm{LacI}\right]=\frac{\left[{\mathrm{LacI}}_{\mathrm{total}}\right]}{\mathrm{1\; +}\frac{\left[\mathrm{IPTG}\right]}{K_{\mathrm{LacI\; -\; IPTG}}}}$$ which finally gives our differential equations for TetR : $$\frac{d\left[\mathrm{TetR}\right]}{dt}=\frac{k_{\mathrm{plac}}\left[P_{\mathrm{lac\; total}}\right]}{\mathrm{1\; +}{\left(\frac{\left[\mathrm{LacI}\right]}{K_{\mathrm{plac}}+\frac{K_{\mathrm{plac}}\left[\mathrm{IPTG}\right]}{K_{\mathrm{LacI\; -\; IPTG}}}}\right)}^{n_{\mathrm{plac}}}}-{\delta }_{\mathrm{TetR}}\left[\mathrm{TetR}\right]$$

The equation can be generalized to account for cooperativity arising from multimerization of the transcription factors, here represented by the cooperativity constant n_plac

In a similar way, we obtain a differential equation for LacI: $$\frac{d\left[\mathrm{LacI}\right]}{dt}=\frac{k_{\mathrm{pTet}}\left[P_{\mathrm{Tet\; total}}\right]}{\mathrm{1\; +}{\left(\frac{\left[\mathrm{TetR}\right]}{K_{\mathrm{pTet}}+\frac{K_{\mathrm{pTet}}\left[\mathrm{aTc}\right]}{K_{\mathrm{TetR\; -\; aTc}}}}\right)}^{n_{\mathrm{pTet}}}}-{\delta }_{\mathrm{LacI}}\left[\mathrm{LacI}\right]$$

These two equations can be easily computed with a differential solver. We can precisely estimate the effects of each parameter.

Our equations and how we obtained them

Quorum sensing

1. Mathematical Models

Our work mainly refers to the models set up by the 2007 iGEM Bangalore team. On the basis of their work we set up models adapted to our own system. For these reasons we strongly recommend getting familiar with the works of the 2007 Bangalore team for an easier understanding of the models we used.
The main difference between our models is that their model is designed for a whole medium, in which the concen- trations of quorum sensing molecules are considered for a whole fixed volume of a medium. Our system, however, is supposed to describe the spatial diffusion of quorum sensing molecules as well, and therefore needs to be designed for an infinitesimal volume of medium containing bacteries and outside medium. A few other differences exist between our model and theirs, mainly due to the fact that the system we intend to describe is made of other different parts. For example the production rate of our Quorum Sensing enzymes are directly linked to the previously described toggle switch model.

For these reasons we strongly recommend getting familiar with the works of the 2007 Bangalore team for an easier understanding of the models we used.

Bangalore 07 modeled the behaviour of quorum sensing for a simple quorum sensing system. With the input of the toggle switch model taken into account, we can adapt their equations to our system. With our toggle switch system the production would be ruled by the regulatory network of LacI and TetR :

$$\frac{d\left[\mathrm{CinI}\right]}{dt}=\frac{k_{\mathrm{pTet}}\left[P_{\mathrm{Tet\; total}}\right]}{\mathrm{1\; +}{\left(\frac{\left[\mathrm{TetR}\right]}{K_{\mathrm{pTet}}+\frac{K_{\mathrm{pTet}}\left[\mathrm{aTc}\right]}{K_{\mathrm{TetR\; -\; aTc}}}}\right)}^{n_{\mathrm{pTet}}}}-{\delta }_{\mathrm{CinI}}\left[\mathrm{CinI}\right]$$ $$\frac{d\left[{\mathrm{CinR}}_{\mathrm{free}}\right]}{dt}=\frac{k_{\mathrm{plac}}\left[P_{\mathrm{lac\; total}}\right]}{\mathrm{1\; +}{\left(\frac{\left[\mathrm{LacI}\right]}{K_{\mathrm{plac}}+\frac{K_{\mathrm{plac}}\left[\mathrm{IPTG}\right]}{K_{\mathrm{LacI\; -\; IPTG}}}}\right)}^{n_{\mathrm{plac}}}}-{\delta }_{\mathrm{CinR}}\left[{\mathrm{CinR}}_{\mathrm{free}}\right]-V_{\mathrm{complexation}}$$

We can therefore describe the production of CinI and CinR_free inside the cells. V_complexation is the complexation rate of CinR_free with AHL molecule. As a matter of fact CinR will be transformed into CinR* after being complexed with the Quorum Sensing molecules entering the cell. It is now taken into account via this complexation rate.
A simple way to write this rate would be as follow : $$V_{\mathrm{complexation}}=k_{\mathrm{comp}}{Q}_i^n\left[{\mathrm{CinR}}_{\mathrm{free}}\right]$$ with Qi being the concentration in QS molecule inside the cell. If we consider that only one QS molecule would bind to a CinR molecule, we obtain the following equation for CinR : $$\frac{d\left[{\mathrm{CinR}}_{\mathrm{free}}\right]}{dt}=\frac{k_{\mathrm{plac}}\left[P_{\mathrm{lac\; total}}\right]}{\mathrm{1\; +}{\left(\frac{\left[\mathrm{LacI}\right]}{\mathrm{1\; +}\frac{\left[\mathrm{IPTG}\right]}{K_{\mathrm{LacI\; -\; IPTG}}}}\right)}^{n_{\mathrm{plac}}}}-{\delta }_{\mathrm{CinR}}\left[{\mathrm{CinR}}_{\mathrm{free}}\right]-k_{\mathrm{comp}}{Q}_i^{}\left[{\mathrm{CinR}}_{\mathrm{free}}\right]$$

For the following equations the physical volume considered is an infinitesimal volume of medium along x - i.e. we only consider an l ∗ dx volume of cell, l being the width of our plate and dx an infinitesimal portion of length. In this infinitesimal volume we set a fixed number of non-growing cells and take into account the diffusion from one portion to the next ones.

$$\frac{d\left[\mathrm{Qi}\right]}{dt}=\eta \left(\left[\mathrm{Qe}\right]-\left[\mathrm{Qi}\right]\right)-{\delta }_{\mathrm{Qi}}\left[\mathrm{Qi}\right]+\mathrm{f([CinI])}$$ $$\frac{d\left[\mathrm{Qe}\right]}{dt}=\rho v_c\eta \left(\left[\mathrm{Qi}\right]-\left[\mathrm{Qe}\right]\right)-{\delta }_{\mathrm{Qe}}\left[\mathrm{Qe}\right]+D_{\mathrm{diff}}\frac{{\∂}^2\left[\mathrm{Qe}\right]}{\∂x^2}$$
• With f ([CinI]) being a mathematical function describing the production of QS molecule by CinI enzyme. Basically this fonction would be as follow :
$$\mathrm{f([CinI])\; =}{k}_{{\mathrm{QS}}_p}{\left[\mathrm{substrate}\right]}^n\left[\mathrm{CinI}\right]$$ But if this reaction is assumed first-order, we obtain : $$\frac{d\left[\mathrm{Qi}\right]}{dt}=\eta \left(\left[\mathrm{Qe}\right]-\left[\mathrm{Qi}\right]\right)-{\delta }_{\mathrm{Qi}}\left[\mathrm{Qi}\right]+k_{{\mathrm{QS}}_p}^{\text{'}}\left[\mathrm{CinI}\right]$$
• With Ddiff being the diffusion coefficient for our Quorum sensing molecule in our medium along spatial dimension x.
• In our case ρvc is a constant (we consider the cells do not grow in our time scale)

With the equations set (3.1); (3.3); (3.5); (3.6) we have, we can not use solvers like matlab ODE because of their space and time dependancies. To solve our problem we have to use a space-time derivation matrix we will describe in the next chapter.

2. Solvation ot the set of equations (3.1); (3.3); (3.5); (3.6)

3. To solve this set of equations we have to use a matrix that will describe our system in both space and time. for example for the QS molecule outside of the cell :

   $$M_{Q_e}\left(\mathrm{m,n}\right)=\left[Q_e\right]\left(\mathrm{x,t}\right)$$ $$M_{Q_e}\left(\mathrm{m+1,n+1}\right)=\left[Q_e\right]\left(\mathrm{x+dx,t+dt}\right)$$
   • On the spatial point of view, we only consider the x dimension, as the IPTG gradient will be only evolving along this dimension. Thus we consider the state of our cells is the same along the width of our plate.
   • With this Matrix, and after computation of all the terms, we can get the entire behaviour of CinI, CinR, QS inside and outside the cells.
   • The first line of the Matrix equals 0. These are the initial conditions we set to 0 at time t = 0s.
   • On the borders of the plate (x = 0 and x = L) the model used has to be different, limit conditions will be set.
   • Of course, Qi, CinI and CinR matrices will be similarly implemented.

   With our continuous equations set, we want to obtain discrete definition of each of the matrices. The interdepen- dancies of the equations imply that the computation of the matrices will be performed on the entire CinI matrix first, then each line of the Qi and Qe matrices will be computed alternatively. Finally CinR matrix computation will be performed.

   Parallel computation of all the matrices without proper control is not possible indeed, as the terms of Qi matrix will depend on the Qe terms of the preceding line (and vice-versa). Discretization is obtained with first order taylor series :

   $$M_{Q_i}\left(\mathrm{m,n+1}\right)=\Delta t\left(\eta \left(M_{Q_e}\mathrm{(m,n)}-M_{Q_i}\left(\mathrm{m,n}\right)\right)-{\delta }_{Q_i}{M}_{Q_i}\left(\mathrm{m,n}\right)+k_{{\mathrm{QS}}_p}{M}_{\mathrm{CinI}}\left(\mathrm{m,n}\right)+M_{Q_i}\left(\mathrm{m,n}\right)\right)$$ $$M_{Q_e}\left(\mathrm{m,n+1}\right)=\Delta t\left(D_m+D_{\mathrm{diff}}\frac{M_{Q_e}\left(\mathrm{m+1,n}\right)-{\mathrm{2M}}_{Q_e}\left(\mathrm{m,n}\right)+M_{Q_e}\left(\mathrm{m-1,n}\right)}{}\right)-M_{Q_e}\left(\mathrm{m,n}\right)$$

   with $D_m=\rho v_c\eta M_{Q_i}\left(\mathrm{m,n}\right)-M_{Q_e}\left(\mathrm{m,n}\right)\left({\delta }_{Q_e}+\rho v_c\eta \right)$

   $$M_{\mathrm{CinR}}\left(\mathrm{m,n+1}\right)=\Delta t\left(\frac{k_{\mathrm{plac}}\left[P_{\mathrm{lac\; total}}\right]}{\mathrm{1\; +}{\left(\frac{\left[\mathrm{LacI}\right]}{K_{\mathrm{plac}}+\frac{K_{\mathrm{plac}}\left[\mathrm{IPTG}\right]}{K_{\mathrm{LacI\; -\; IPTG}}}}\right)}^{n_{\mathrm{plac}}}}-M_{\mathrm{CinR}}\left(\mathrm{m,n}\right)\left({\delta }_{\mathrm{CinR}}-k_{\mathrm{comp}}{M}_{Q_i}\left(\mathrm{m,n}\right)\right)\right)+M_{\mathrm{CinR}}\left(\mathrm{m,n}\right)$$

With these discrete equations the 4 matrices can be computed through simple calculation loops over each line. The CinI matrix does not depend on space dimension, it is then possible to compute it without discretization with a differential solver.

Our algorithms

In the MATLAB archive that can be found here containing our matlab scripts for deterministic modelling (file Deterministic_archive.tar.gz) you can launch an ODE based simulation (see our ODEs in the two previous sections) with the file biosenseur1Dmain.m.

Several dialog boxes will pop up to enter the specificities of the simulation : (physical specificities of the device, chemical species concentrations and IPTG gradient )

At the end of the simulation you obtain 3 matrices named M_stock, M_QS and M_comp containing the concentrations in each protein species at each time point and on each physical point of the plate. We wrote three MATLAB scripts that display the concentration in proteins dynamically that you can call with DynamicplottingTS, DynamicplottingQS and DynamicplottingCP. For a good understanding of the models and of our results we also wrote a script to illustrate the coloration of our plate through time according to our models named Imageshow.m.

Isoclines and Hysteresis

In order to predict the stability of our system, mathematical studies are necessary.

Isoclines

In a first time, an isocline study is realized. Isocline study is a classical study which implies a research of stationnary point in a system. This stationnary point are deduced from the equations of the differential system: when the variation of concentration of both repressors are equal to zero.

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

• E_TetR = kplac * [pLac]tot
• $$RtetR=\frac{1_{}}{\mathrm{1\; +}{\left(\frac{\left[\mathrm{TetR}\right]}{K_{\mathrm{pTet}}+\frac{K_{\mathrm{pTet}}\left[\mathrm{aTc}\right]}{K_{\mathrm{TetR\; -\; aTc}}}}\right)}^{n_{\mathrm{pTet}}}}$$
• E_LacI = kpTet * [pTet]tot
• $$RLacI=\frac{1_{}}{\mathrm{1\; +}{\left(\frac{\left[\mathrm{LacI}\right]}{K_{\mathrm{pLac}}+\frac{K_{\mathrm{pLac}}\left[\mathrm{IPTG}\right]}{K_{\mathrm{LacI\; -\; IPTG}}}}\right)}^{n_{\mathrm{pLac}}}}$$
• [TetR]r = RTetR * [TetR] the relative concentration of TetR
• [LacI]r = RLacI * [LacI] the relative concentration of LacI
• $K=\frac{\mathrm{RtetR\; *}\mathrm{EtetR}}{\delta \mathrm{TetR}}$
• $K\text{'}=\frac{\mathrm{RlacI\; *}\mathrm{ElacI}}{\delta \mathrm{LacI}}$

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

From these equations, we get these figures:

On this figure, the red lines represent the equation(1) of the system and the green lines the second equation(2). This figure was obtained with [aTc] = 1E-6 M and [IPTG] = 1,55 E-4 M. These parameters reflect the situation of our system in the center of the plate to a logarithmic gradient of IPTG of 1E-6 M to 1E-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 middle one. It represents the equality of relative concentration. In a toggle switch, it's impossible to have the equality of concentration of both repressors because one repressed the other. So one of these should take advantage of the other.

The left graph represents the left side of the plate: aTc concentration is superior to IPTG concentration and the right graph represents the right side of the plate where IPTG concentration prevails. These figures show that when the concentration of one of the repressor is majoritary, the system is no longer bistable but monostable.

Points to Remember :

• 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. Thereby bacteria could communicate together and active the expression of coloration.

Hysteresis

The goal of the hysteresis study is to established the switch conditions when the toggle switch is already blocked in a predefined way. In our mathematical studies, we blocked the system in the lacI way with different preliminary aTc concentrations. To get hysteresis curves, after blocking the system, we grow the amount of IPTG to see from which concentration the system will switch. Then, we decrease IPTG concentration to see if the system switch back to the initial state at the same IPTG concentration. 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.

Important values to exploit from these curves are the two switch concentration of IPTG. On this curve, we can get the switch up concentration: ~ 1E-2 M of IPTG and the switch back concentration ~ 3E-5 M.

The switch back concentration is very similar to the dissociation constant between lacI and IPTG (which is 2.96E-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.

The two first curves (for aTc = 1E-9 M and aTc = 1E-8 M) show that the switch up and switch back concentrations are the same. This concentration is ~ 3E-5 M, the dissociation constant between lacI and IPTG. The last curve (for aTc = 1E-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 1E-7 M, the switch up concentration is growing with aTc concentration.

The concentration of aTc 1E-7 M appears to be the limit of sensibility to the toggle switch. This concentration represents the dissociation constant of aTc to TetR repressor.

Points to Remember :

From this mathematical study, we can get the limit of detection of our system: which is close to the dissociation constant of aTc and is repressor.