Team:Paris Liliane Bettencourt/Project/Population counter/model

From 2010.igem.org

(Difference between revisions)
 
(34 intermediate revisions not shown)
Line 7: Line 7:
</a>
</a>
  <font size=4>Population counter</font>
  <font size=4>Population counter</font>
 +
<a href="https://2010.igem.org/Team:Paris_Liliane_Bettencourt/Project/Synbioworld">
 +
<img src="https://static.igem.org/mediawiki/2010/2/25/SBW.jpg" width="129" height="107" align=right title="SynBioWorld">
 +
</a>
<a href="https://2010.igem.org/Team:Paris_Liliane_Bettencourt/Project/SIP">
<a href="https://2010.igem.org/Team:Paris_Liliane_Bettencourt/Project/SIP">
  <img src="https://static.igem.org/mediawiki/2010/4/4c/SIP.png" width="132" height="107" align=right title="SIP">
  <img src="https://static.igem.org/mediawiki/2010/4/4c/SIP.png" width="132" height="107" align=right title="SIP">
Line 28: Line 31:
   <li><a href="https://2010.igem.org/Team:Paris_Liliane_Bettencourt/Project/Population_counter/Results" target="_self">Results</a></li>
   <li><a href="https://2010.igem.org/Team:Paris_Liliane_Bettencourt/Project/Population_counter/Results" target="_self">Results</a></li>
   <li><a href="https://2010.igem.org/Team:Paris_Liliane_Bettencourt/Project/Parts" target="_self">Parts</a></li>
   <li><a href="https://2010.igem.org/Team:Paris_Liliane_Bettencourt/Project/Parts" target="_self">Parts</a></li>
 +
  <li><a href="https://2010.igem.org/Team:Paris_Liliane_Bettencourt/Project/Population_counter/References" target="_self">References</a></li>
</ul>
</ul>
</div>
</div>
Line 40: Line 44:
<h2> The process </h2>
<h2> The process </h2>
-
Here we will give the bigger picture, discussing how every component of the process reacts and interacts with each other. First, we send a pulse of arabinose through the tunnel. The bacteria are then expected to react in a way that is described in the next section. Those bacteria that recombinated start producing LuxI, which in turns promotes the production of AHL. This AHL will play two roles: <br />
+
To be able to model the system, we first need to have a feeling of how the components of the system react and interact with one another. First, we send a pulse of arabinose through the tunnel which triggers, with a certain probability, the recombination of some bacteria. Those bacteria are then expected to become red (RFP expression) and start producing LuxI, which in turns promotes the production of AHL. AHL will play two roles: <br />
<ul>
<ul>
-
<li> it can bind with LuxR in the cell. LuxR is a protein produced so that its concentration inside the cell will be considered as constant. We call it LuxR<sup>f</sup> as long as it is free from AHL, LuxR* otherwise;</li>
+
<li> it can bind with intracellular LuxR, which is a constitutively expressed protein so that its concentration inside the cell will be considered as constant. We call it LuxR<sup>f</sup> as long as it is free from AHL, LuxR* otherwise;</li>
<li> it can cross the cell membrane to bind with LuxR of other bacteria in the device.</li><br />
<li> it can cross the cell membrane to bind with LuxR of other bacteria in the device.</li><br />
</ul>
</ul>
-
Every time a pulse of arabinose is injected, the production of AHL increases. Once a certain level of AHL has been reached, all the cells produce GFP, which is the event putting an end to the counting.<br />
+
Every time a pulse of arabinose is injected, the production of AHL increases as a result of more recombinations. Once a certain level of AHL has been reached, all the cells respond to the quorum sensing signal and start producing GFP, which is the event putting an end to the counting.<br />
<center>
<center>
<img src="https://static.igem.org/mediawiki/2010/4/40/Mf_device.jpg" alt="Mf device" title="Scheme of the microfluidic device" align="center" height="300"/><br />Figure 1: Scheme of the microfluidic device</center><br />
<img src="https://static.igem.org/mediawiki/2010/4/40/Mf_device.jpg" alt="Mf device" title="Scheme of the microfluidic device" align="center" height="300"/><br />Figure 1: Scheme of the microfluidic device</center><br />
-
<h2> Model </h2>
+
<h2> The physical model </h2>
We can develop the information we have in a system of equations:<br />
We can develop the information we have in a system of equations:<br />
<ul>
<ul>
-
<li> <h5>LuxI</h5> In a first approach, we will consider that the dynamics of production of LuxI is represented by a step function: [LuxI] is 0 before the recombination, and 1 after. If the time between the two states cannot be neglected compared to the average time between two recombinations, we may have to reconsider this simple model. (which is probably the case since it is roughly the expected time of a cell recombinating eg close to 20 minutes, way larger than the diffusion significant times)</li>
+
<li> <h5>Recombination</h5>
-
<li> <h5> LuxR </h5>Let's call <em>R</em> the constant concentration of LuxR in a cell (how  <em>"constant"</em> is it?). Some of it is binded with AHL and we called it LuxR*, so this gives:<br />
+
In all our simulations, we assumed that recombination would happened after 20 minutes in contact of arabinose, plus a certain random time following an exponential law of average 10 minutes. After this random time, it is not certain that the bacteria will recombinate. In our simulations, we decided that our cell would recombinate with a probability of 0.05
 +
<li> <h5>LuxI</h5>
 +
Once a bacteria has recombinated, there is necesarily a latent time before LuxI is actually synthetised. There is then some more time before a bacteria produces full rate. The amount of LuxI produced inside a bacteria stays in this very bacteria, and it is commonly accepted that the concentration of LuxI reaches a plateau and keep this value. According to the data found in the literature, we chose to set the latent time to 10 minutes. Then the concentration increases during 10 minutes linearly until the maximum concentration is reached.
 +
</li>
 +
<li> <h5> LuxR </h5>Let's call <em>R</em> the constant concentration of LuxR in a cell. Some of it binded with AHL and we call it LuxR*, so this gives:<br />
<center>
<center>
Line 75: Line 83:
<li> <h5>AHL</h5> Following the paper "A synchronized quorum of genetic clocks" by Hasty and al., we define the intracellular concentration of AHL, [AHL<sub>i</sub>], and the external one, [AHL<sub>e</sub>].<br />
<li> <h5>AHL</h5> Following the paper "A synchronized quorum of genetic clocks" by Hasty and al., we define the intracellular concentration of AHL, [AHL<sub>i</sub>], and the external one, [AHL<sub>e</sub>].<br />
-
The concentration of AHL in the room [AHL<sub>e</sub>] decreases as the microfluidic device evacuates it, which gives birth to a diffusion term in the equation caracterized by the kinetic constant kdiff. In parallel, the internal and external AHL concentrations change as AHL is produced by the cells and can go throuh the cell membrane (kinetic constant kI and kmembrane respectively). For now, we shall assume the diffusion of AHL inside the device is instantaneous.<br/>
+
Internal and external AHL concentrations change as AHL is produced by the cells and can go throuh the cell membrane (kinetic constant k<sub>I</sub> and k<sub>membrane</sub> respectively).<br/>
<center>
<center>
<img src="https://static.igem.org/mediawiki/2010/2/22/Eq_AHLi.jpg" alt="Eq_AHLi" title="Equation for the intra-cellular AHL" align="center" height="50"/></center><br />
<img src="https://static.igem.org/mediawiki/2010/2/22/Eq_AHLi.jpg" alt="Eq_AHLi" title="Equation for the intra-cellular AHL" align="center" height="50"/></center><br />
-
The first term says that if the cell produces LuxI, then it produces AHL at a certain rate. The second term concerns the exchange between the cell and the medium, and is proportional to the difference of concentration.<br />
+
The first term says that if the cell produces LuxI, then it produces AHL at a certain rate. The second term concerns the exchange between the cell and the medium, and is proportional to the difference of concentration. The concentration of AHL in the room [AHL<sub>e</sub>] decreases as the microfluidic device evacuates it, which gives birth to a diffusion term in the equation caracterized by the kinetic constant k<sub>diff</sub>.<br />
<center>
<center>
<img src="https://static.igem.org/mediawiki/2010/a/a9/Eq_AHLe.jpg" alt="Eq_AHLe" title="Equation for the extra-cellular AHL" align="center" height="50"/></center><br />
<img src="https://static.igem.org/mediawiki/2010/a/a9/Eq_AHLe.jpg" alt="Eq_AHLe" title="Equation for the extra-cellular AHL" align="center" height="50"/></center><br />
-
where k<sub>diff</sub> is the constant of diffusion, related to parameters of the fluid and of the device, k<sub>membrane</sub> is the same as previously, but is now compensated by <em>d</em>, the cell density, so that the exchanges between cells and the medium is balanced. A note on diffusion further.<br/>
+
where k<sub>diff</sub> is related to parameters of the fluid and of the device, k<sub>membrane</sub> is the same as previously, but is now compensated by <em>d</em>, the cell density, so that the exchanges between cells and the medium is balanced.<br/>
-
Here we consider only the case where the concentration of AHL is homogeneous, ie the concentration inside and outside the cell are the same. This is justified since the caracteristic time of diffusion of AHL between the outside and the inside of a cell is short compared k<sub>I</sub>. The contrary would lead to a discussion we started in paragraph \ref{modif_model}. We end up with this equation:<br/>
+
Here we consider only the case where the concentration of AHL is homogeneous, ie the concentration inside and outside the cell are the same. This is justified since the caracteristic time of diffusion of AHL between the outside and the inside of a cell is short compared k<sub>I</sub> and k<sub>diff</sub>. We end up with this equation:<br/>
<center>
<center>
<img src="https://static.igem.org/mediawiki/2010/5/51/Eq_AHL.jpg" alt="AHL" title="Diffusion equation, AHL" align="center" height="60"/></center><br />
<img src="https://static.igem.org/mediawiki/2010/5/51/Eq_AHL.jpg" alt="AHL" title="Diffusion equation, AHL" align="center" height="60"/></center><br />
 +
In Hasty et al., k<sub>I</sub> is taken to be <em>0.06</em>, which value we kept. Then, k<sub>diff</sub> would have had to be determined experimentally. Instead, we followed Hasty et al.'s demarche, making the computations with some "reasonable" values.
 +
</li>
</li>
</ul>
</ul>
<h3> The synthesis of LuxI.</h3>
<h3> The synthesis of LuxI.</h3>
-
  Above, we introduced directly the concentration of LuxI without giving more details. However, we need to model the production of LuxI to calibrate the counter. This means modelling how the bacteria evolve in the device and how they react to arabinose.
+
  Above, we introduced directly the concentration of LuxI without giving more details. However, we need to model the production of LuxI to calibrate the counter. This means modelling how the bacteria evolve in the device and how they react to arabinose.<br/><br/>
<ul>
<ul>
-
<li> Arabinose is modeled simply in the model with a diffusion equation: <br/>
+
<li> <em>Arabinose</em> is modeled simply in the model with a diffusion equation: <br/><br/>
<center>
<center>
<img src="https://static.igem.org/mediawiki/2010/f/ff/Eq_arabinose.jpg" alt="arabinose" title="Diffusion equation, arabinose" align="center" height="53"/></center><br />
<img src="https://static.igem.org/mediawiki/2010/f/ff/Eq_arabinose.jpg" alt="arabinose" title="Diffusion equation, arabinose" align="center" height="53"/></center><br />
-
with the boundary condition at the border between the device and the tunnel corresponding to a positive value if a pulse is being sent, and 0 otherwise. Thus we get the profile of arabinose independently of everything else, except the pulse itself.
+
with the boundary condition at the border between the device and the tunnel corresponding to: a positive value if a pulse is being sent, and 0 otherwise. Thus we get the profile of arabinose independently of everything else, except the pulse itself. k<sub>A</sub> was chosen arbitrarily (but close to values for similar molecules found in the literature) to be <em>10<sup>-5</sup> mm<sup>2</sup>.min<sup>-1</sup></em> . The pulse of arabinose in the channel corresponds to a concentration of <em>0.04%</em>, which concentration was used in the lab.<br/>
 +
To trigger the production of LuxI, a bacteria needs to be exposed to a minimum concentration <em>C<sub>c</sub></em> of arabinose for at least <em>T<sub>c</sub></em> minutes. For the sake of the simulation, we chose values for which we knew that recombinations would happen, here <em>C<sub>c</sub> = 0.001%</em> during <em>T<sub>c</sub> = 1.5 min</em>.
</li>
</li>
-
 
+
<br/>
-
<li> The device may have several sizes, however one device will not ever host more than 2.10<sup>5</sup> cells. We can indeed approximate the maximum number of bacteria as the ratio between the volume of a device and the volume of a cell. Moreover, considering both populations, switched and unswitched, we will consider the bacteria are in exponential growth, and their excess is evacuated through the path to the tunnel. Therefore, we need a model for the behaviour of the bacteria (movement and birth/death) that allows us to deal with a few cells \footnote{at the beginning, when we only have a few recombinated cells.} up to more than 10<sup>5</sup> cells. An individually-centered model would be computationnaly inefficient, while a continuous model would perform wrong computations for small numbers of bacteria.  
+
<li> Our devices have several sizes, however one device will not ever host more than 2.10<sup>5</sup> cells. We can indeed approximate the maximum number of bacteria as the ratio between the volume of a device and the volume of a cell. Moreover, considering both populations, switched and unswitched, we will consider the bacteria are in exponential growth, and their excess is evacuated through the path to the tunnel. Therefore, we need a model for the behaviour of the bacteria (movement and birth/death) that allows us to deal with a few cells (when we have a few recombinated cells for instance) up to more than 10<sup>5</sup> cells. An individually-centered model would be computationnaly inefficient, while a continuous model would perform wrong computations for small numbers of bacteria.<br/>
-
I have implemented something relying on a binomial distribution: among n people, we pick a random number of them to move, or give birth; this random number being ruled by a probability p proportional to the gradient of concentration of bacteria. Discretising time and space, one can draw at every step time a random binomial number\footnote{X being non-negative if the first argument in the binomial law below is non-negative, negative otherwise}:<br/>
+
Birth and death happen randomly, but we know how often they happen in average. We will describe in the next section the algorithm we used to model this.<br/>
-
<img src="https://static.igem.org/mediawiki/2010/f/f7/Bin_distrib.jpg" alt="bin_distrib" title="Binary distribution, ruling the movement of the bacteria" align="center" height="35"/><br />
+
As for diffusion of bacteria, a model relying on a binomial distribution was implemented: among n cells, we pick a random number of them to move; this random number being ruled by a probability p proportional to the gradient of concentration of bacteria. Discretising time and space, one can draw at every step time a random binomial number:<br/><br/>
-
where...
+
<center>
 +
<img src="https://static.igem.org/mediawiki/2010/f/f7/Bin_distrib.jpg" alt="bin_distrib" title="Binary distribution, ruling the movement of the bacteria" align="center" height="35"/></center><br />
 +
where <em>f</em> and <em>f<sub>2</sub></em> are described in the next section.
</li>
</li>
 +
<br/><br/>
<h2> Algorithm and simulations </h2>
<h2> Algorithm and simulations </h2>
In this paragraph, we use several constants suggested in the <em>paper of Hasty et al.</em> (see references), or sometimes some approximations. We shall mention these numerical values along this paragraph.<br/><br/>
In this paragraph, we use several constants suggested in the <em>paper of Hasty et al.</em> (see references), or sometimes some approximations. We shall mention these numerical values along this paragraph.<br/><br/>
 +
We use a discrete lattice to make the computations, considering that the cavity of the microfluidic device can be divided in several slices of constant width <em>dx</em>. It has been further assumed (for simplicity of computations essentially) that diffusion only takes place along the <em>x</em> axis (see figure 1 above). The discrete step time wil be denoted <em>dt</em>. Thus our equations will depend on a one-dimensional space parameter and on time.<br/><br/>
<h3> Solving the diffusion equation using an implicit scheme </h3>
<h3> Solving the diffusion equation using an implicit scheme </h3>
-
mainly refers to NRC
+
One can simply solve numerically a partial differential equation by discrete differentiation. The simplest way gives directly the solution at time step n+1 thanks to what we know at time n. However, this discrete explicit scheme gives a stable solution only if <em>k<sub>diff</sub>.dt/dx<sup>2</sup></em> is less than 1. Therefore, we considered the Crank-Nicholson discretisation scheme. Unconditional to <em>k<sub>diff</sub></em>, <em>dt</em> and <em>dx</em>, this scheme produces stable computations.<br/><br/>
 +
 
 +
<center>
 +
<img src="https://static.igem.org/mediawiki/2010/6/60/Discrete_scheme.jpg" alt="bin_distrib" title="Binary distribution, ruling the movement of the bacteria" align="center" height="100"/></center><br />
 +
 
 +
However, it appears to be more complicated to solve since <em>u(x,t+dt)</em> cannot be obtained explicitly knowing the solution <em>u</em> at time <em>t</em>. Fortunately it can be expressed as a tridiagonal system of equations that we solved numerically.<br/><br/>
<h3> Stochastic Gillespie algorithm </h3>
<h3> Stochastic Gillespie algorithm </h3>
 +
The basic idea lying here is to compute random times for a population rather than computing them for each individual (which is much more computer intensive). Those random times appear in the recombinations and the birth-death calculus; Gillespie's idea, which we are about to describe, allows to compute the time of the next recombination (resp. birth or death) within the population instead of making the calculus and storing the result for each cell in the device.<br/>
 +
Under the assumption that cells' recombination times are independent from one bacteria to the other, and that these times follow an exponential law of parameter &lambda;, Gillespie tells us that the next recombination time (resp. birth or death time) follows an exponential law of parameter <em>N<sub>pop</sub>.&lambda;</em>.<br/><br/>
-
<h3> else </h3>
 
-
<img src="https://static.igem.org/mediawiki/2010/e/e5/Pic_soft1.jpg" alt="soft1" title="Our simulation tool, stage 1" align="center" height="440"/><br /><br />
+
<h3> Stochastic movement of the cells </h3>
 +
As promised before, we give the two functions describing the movement of the bacteria in the microfluidic device:
 +
<br/><br/>
 +
<center>
 +
<img src="https://static.igem.org/mediawiki/2010/f/f0/F.jpg" alt="bin_distrib" title="Binary distribution, ruling the movement of the bacteria" align="center" height="34"/></center><br />
-
<img src="https://static.igem.org/mediawiki/2010/b/be/Pic_soft2.jpg" alt="soft2" title="Our simulation tool, stage 2" align="center" height="400"/><br /><br />
+
<center>
 +
<img src="https://static.igem.org/mediawiki/2010/5/52/F2.jpg" alt="bin_distrib" title="Binary distribution, ruling the movement of the bacteria" align="center" height="36"/></center><br />
 +
where <em>N<sub>sw/unsw</sub></em> denotes either the number of switched bacteria (or recombinated bacteria) or unswitched bacteria. The idea here is to have, on average at each step time, <em>dt.k<sub>diff</sub>.f(x,t)</em> bacteria that will move, a result we would get from usual diffusion. As previously mentioned, we cannot use regular diffusion since our algorithm must apply to discrete data lying on a wide range of integers.
 +
<br/><br/>
 +
<h3> Simulations </h3>
 +
Now let's present the graphical display of our simulation tool. Once the user has selected a refreshment speed x, the software will draw every x minutes the unswitched bacteria in grey, the switched ones in red and the arabinose in small blue points.<br/>
 +
<br/>
 +
<center>
 +
<img src="https://static.igem.org/mediawiki/2010/e/e5/Pic_soft1.jpg" alt="soft1" title="Our simulation tool, stage 1" align="center" height="460"/></center><br /><br />
 +
We see on the following picture that after a certain amount of time, both cells and arabinose are distributed rather uniformly, and that the concentration of arabinose decreases near the tunnel as it is washed away by the flow. <br/><br/>
 +
<center>
 +
<img src="https://static.igem.org/mediawiki/2010/b/be/Pic_soft2.jpg" alt="soft2" title="Our simulation tool, stage 2" align="center" height="430"/></center><br /><br />
-
For now, we use a simple spatial discretisation scheme, assuming lateral diffusion is instantaneous\footnote{this may make simulations inaccurate and thus yields a piece of the program to be improved}, so that diffusion is only considered along one axis (the one from the tunnel to the top of the device (see the first figure).<br/>
+
The software also gives us the dynamics of recombined bacteria, LuxI and AHL slice by slice. This is what we will exploit here.<br/>
-
<img src="https://static.igem.org/mediawiki/2010/3/31/Chart_AHL.jpg" alt="AHL_chart" title="Number of switched cells, LuxI and AHL against time" align="center" height="500"/><br /><br />
+
In the following simulations, the device was divided into 15 slices. We exhibit some results on the fifth slice. To count, one can either look at the number of recombined bacteria thanks to RFP, or set up the timer thanks to AHL. We will see that depending on the parameters, one or other may not be possible. For the first one, it is necessary to have clear steps, implying that the time between the pulses will be large enough to allow the system to stabilise. For the second, distinct steps on the AHL graph are necessary to set up the timer. Indeed, this allows to design the bacteria in the way that they will produce GFP only once a certain treshold in AHL was reached. However, admitting that this timer is set up, then it is possible to count more frequent events even without steps as we shall see in figure 4. In the sequel, the time between two pulses will be denoted <em>T<sub>2</sub></em>.
 +
<br/><br/>
-
<img src="https://static.igem.org/mediawiki/2010/5/5d/Chart_AHL3.jpg" alt="AHL_chart3" title="Number of switched cells, LuxI and AHL against time" align="center" height="600"/><br /><br />
+
<center>
 +
<img src="https://static.igem.org/mediawiki/2010/e/ec/Chart_AHL4.jpg" alt="AHL_chart" title="Number of switched cells, LuxI and AHL against time" align="center" height="500"/><br/> Figure 2: Switched bact., AHL, LuxI for T<sub>2</sub> = 160 and k<sub>diff</sub> = 8.10<sup>-2</sup></center><br /><br />
 +
In the following figure, we show that it is still possible to count the number of pulses thanks to the number of recombined bacteria, although it might become uncertain. Nevertheless, it seems that we can still choose a threshold of AHL beyond which bacteria should be designed to luminesce. For instance, if we wanted our timer to stop at 4, then a threshold of 0.70 in standardised unit seems reasonable.
 +
<center>
 +
<img src="https://static.igem.org/mediawiki/2010/5/5d/Chart_AHL3.jpg" alt="AHL_chart3" title="Number of switched cells, LuxI and AHL against time" align="center" height="600"/><br/> Figure 3: Switched bact., AHL, LuxI for T<sub>2</sub> = 120 and k<sub>diff</sub> = 2.10<sup>-2</sup></center><br /><br />
-
Discuss the existence of two phases, the critical parameters being the diffusion coefficient of AHL and the time between two pulses.<br />
 
-
Although we have been able to count six pulses thanks to the previous graph of AHL, it is not always possible to do so. For instance, if we shorten the time between two pulses, keeping the same diffusion coefficient, then it is no longer possible to distinguish pulses:<br />
 
 +
Although we have been able to count six pulses thanks to the previous graph of AHL, it is not always possible to do so. For instance, if we shorten the time between two pulses, decreasing the diffusion coefficient, then it is no longer possible to distinguish pulses:<br />
-
<img src="https://static.igem.org/mediawiki/2010/1/14/Chart_AHL2.jpg" alt="AHL_chart2" title="AHL against time" align="center" height="400"/><br /><br />
+
<center>
 +
<img src="https://static.igem.org/mediawiki/2010/b/be/Chart_AHL5.jpg" alt="AHL_chart2" title="AHL against time" align="center" height="400"/><br/> Figure 4: AHL for T<sub>2</sub> = 100 and k<sub>diff</sub> = 10<sup>-3</sup></center><br /><br />
 +
<h2> Achievements of our model</h2>
 +
 +
Our model gave us a better understanding of the dynamics and helped us explain and suggest experiments <em>a posteriori</em>. Together with the results of the experiments, we may be able:
 +
<ul>
 +
<li> to have a better determination of some parameters such as caracteristic times for recombination and their statistical distribution, which we can then plug in;</li><br/>
 +
<li> to set up parameters in a way that allow us to count with RFP, or to make a timer thanks to AHL and GFP. Those parameters are for instance the speed of the flow in the microfluidic device, the dimensions of the microfluidic device, the initial number of bacteria or the viscosity of the medium fluid.</li>
 +
</ul>
 +
<br /><br />
 +
 +
 +
<h2> Sources of the simulation software </h2>
 +
 +
You may download the sources by clicking <a href="https://2010.igem.org/Image:Simu_paris.zip">here</a>. It is written in C and all the interesting parameters are in the file "constants.h". The files "nrutil.c" and "nrutil.h" are taken from the internet.<br /><br />
<h2> References </h2>
<h2> References </h2>
-
- Hasty...<br/>
+
<ul>
-
- Numerical recipes in C, ??, CUP  <br/>
+
<li> T.Danino, O.Mondragon-Palomino, L.Tsimring, J.Hasty: A synchronized quorum of genetic clocks. Nature 2009, 463:326 - 330</li><br/>
-
- Gillespie... <br/>
+
-
-     <br/><br/>
+
 +
<li> T.Lu, D.Volfson, L.Tsimring, J. Hasty: Cellular growth and division in the Gillespie algorithm. Systems biology, IEEE proceedings 2004, 1:121 - 128 </li><br/>
 +
 +
<li> D.Gillespie: A general method for numerically simulating the stochastic time evolution of coupled chemical reactions. The Journal of Computational Physics 1976, 22:403 - 434 </li> <br/>
 +
 +
<li> D.Gillespie: Exact stochastic simulation of coupled chemical reactions. The Journal of Physical Chemistry 1977, 81:2340 - 2361 </li><br/>
 +
 +
<li> S.Lampoudi, D.Gillespie, L.Petzold: The multinomial simulation algorithm for discrete stochastic simulation of reaction-diffusion systems. 2009 </li><br/>
 +
 +
<li> R.Bustin, H.Messer: An equivalent Markov model for Gillespie’s stochastic simulation algorithm for biochemical systems. 2006 </li><br/>
 +
 +
<li> W.H.Press, S.A.Teukolsky, W.T.Vetterling, B.P.Flannery: Numerical recipes in C. Cambridge University Press, 1992.</li> <br/>
 +
</ul>

Latest revision as of 03:58, 28 October 2010


Population counter





Introduction

We aim at setting up a biological device counting random events occuring in a microfluidic device (see figure 1). These events are recombinations in the cells of the device, stimulated by pulses of arabinose coming from the tunnel. Once a certain number of pulses have been triggered, we expect our device to start emitting green fluorescence. First, we shall describe the dynamics of the population of bacteria in our device as for instance its recombinations, which further lead to the rise of the concentration of produced AHL. AHL is a quorum sensing molecule that, once a certain concentration threshold is reached, triggers in our device the production of GFP. It then alerts us of the end of counting. To count how many recombinations (so how many events) are needed to get the quorum sensing response, we shall determine the role and assess the value of a few critical parameters in our model.

The process

To be able to model the system, we first need to have a feeling of how the components of the system react and interact with one another. First, we send a pulse of arabinose through the tunnel which triggers, with a certain probability, the recombination of some bacteria. Those bacteria are then expected to become red (RFP expression) and start producing LuxI, which in turns promotes the production of AHL. AHL will play two roles:
  • it can bind with intracellular LuxR, which is a constitutively expressed protein so that its concentration inside the cell will be considered as constant. We call it LuxRf as long as it is free from AHL, LuxR* otherwise;
  • it can cross the cell membrane to bind with LuxR of other bacteria in the device.

Every time a pulse of arabinose is injected, the production of AHL increases as a result of more recombinations. Once a certain level of AHL has been reached, all the cells respond to the quorum sensing signal and start producing GFP, which is the event putting an end to the counting.
Mf device
Figure 1: Scheme of the microfluidic device

The physical model

We can develop the information we have in a system of equations:
  • Recombination
    In all our simulations, we assumed that recombination would happened after 20 minutes in contact of arabinose, plus a certain random time following an exponential law of average 10 minutes. After this random time, it is not certain that the bacteria will recombinate. In our simulations, we decided that our cell would recombinate with a probability of 0.05
  • LuxI
    Once a bacteria has recombinated, there is necesarily a latent time before LuxI is actually synthetised. There is then some more time before a bacteria produces full rate. The amount of LuxI produced inside a bacteria stays in this very bacteria, and it is commonly accepted that the concentration of LuxI reaches a plateau and keep this value. According to the data found in the literature, we chose to set the latent time to 10 minutes. Then the concentration increases during 10 minutes linearly until the maximum concentration is reached.
  • LuxR
    Let's call R the constant concentration of LuxR in a cell. Some of it binded with AHL and we call it LuxR*, so this gives:
    Eq_luxR

    The complexation between LuxR and AHL gives a chemical equilibrium determined by the constant Kreac:
    Eq_luxR2

    so at the equilibrium:
    Eq_luxR3

    which gives:
    Eq_luxR4

  • AHL
    Following the paper "A synchronized quorum of genetic clocks" by Hasty and al., we define the intracellular concentration of AHL, [AHLi], and the external one, [AHLe].
    Internal and external AHL concentrations change as AHL is produced by the cells and can go throuh the cell membrane (kinetic constant kI and kmembrane respectively).
    Eq_AHLi

    The first term says that if the cell produces LuxI, then it produces AHL at a certain rate. The second term concerns the exchange between the cell and the medium, and is proportional to the difference of concentration. The concentration of AHL in the room [AHLe] decreases as the microfluidic device evacuates it, which gives birth to a diffusion term in the equation caracterized by the kinetic constant kdiff.
    Eq_AHLe

    where kdiff is related to parameters of the fluid and of the device, kmembrane is the same as previously, but is now compensated by d, the cell density, so that the exchanges between cells and the medium is balanced.
    Here we consider only the case where the concentration of AHL is homogeneous, ie the concentration inside and outside the cell are the same. This is justified since the caracteristic time of diffusion of AHL between the outside and the inside of a cell is short compared kI and kdiff. We end up with this equation:
    AHL

    In Hasty et al., kI is taken to be 0.06, which value we kept. Then, kdiff would have had to be determined experimentally. Instead, we followed Hasty et al.'s demarche, making the computations with some "reasonable" values.

The synthesis of LuxI.

Above, we introduced directly the concentration of LuxI without giving more details. However, we need to model the production of LuxI to calibrate the counter. This means modelling how the bacteria evolve in the device and how they react to arabinose.

  • Arabinose is modeled simply in the model with a diffusion equation:

    arabinose

    with the boundary condition at the border between the device and the tunnel corresponding to: a positive value if a pulse is being sent, and 0 otherwise. Thus we get the profile of arabinose independently of everything else, except the pulse itself. kA was chosen arbitrarily (but close to values for similar molecules found in the literature) to be 10-5 mm2.min-1 . The pulse of arabinose in the channel corresponds to a concentration of 0.04%, which concentration was used in the lab.
    To trigger the production of LuxI, a bacteria needs to be exposed to a minimum concentration Cc of arabinose for at least Tc minutes. For the sake of the simulation, we chose values for which we knew that recombinations would happen, here Cc = 0.001% during Tc = 1.5 min.

  • Our devices have several sizes, however one device will not ever host more than 2.105 cells. We can indeed approximate the maximum number of bacteria as the ratio between the volume of a device and the volume of a cell. Moreover, considering both populations, switched and unswitched, we will consider the bacteria are in exponential growth, and their excess is evacuated through the path to the tunnel. Therefore, we need a model for the behaviour of the bacteria (movement and birth/death) that allows us to deal with a few cells (when we have a few recombinated cells for instance) up to more than 105 cells. An individually-centered model would be computationnaly inefficient, while a continuous model would perform wrong computations for small numbers of bacteria.
    Birth and death happen randomly, but we know how often they happen in average. We will describe in the next section the algorithm we used to model this.
    As for diffusion of bacteria, a model relying on a binomial distribution was implemented: among n cells, we pick a random number of them to move; this random number being ruled by a probability p proportional to the gradient of concentration of bacteria. Discretising time and space, one can draw at every step time a random binomial number:

    bin_distrib

    where f and f2 are described in the next section.


  • Algorithm and simulations

    In this paragraph, we use several constants suggested in the paper of Hasty et al. (see references), or sometimes some approximations. We shall mention these numerical values along this paragraph.

    We use a discrete lattice to make the computations, considering that the cavity of the microfluidic device can be divided in several slices of constant width dx. It has been further assumed (for simplicity of computations essentially) that diffusion only takes place along the x axis (see figure 1 above). The discrete step time wil be denoted dt. Thus our equations will depend on a one-dimensional space parameter and on time.

    Solving the diffusion equation using an implicit scheme

    One can simply solve numerically a partial differential equation by discrete differentiation. The simplest way gives directly the solution at time step n+1 thanks to what we know at time n. However, this discrete explicit scheme gives a stable solution only if kdiff.dt/dx2 is less than 1. Therefore, we considered the Crank-Nicholson discretisation scheme. Unconditional to kdiff, dt and dx, this scheme produces stable computations.

    bin_distrib

    However, it appears to be more complicated to solve since u(x,t+dt) cannot be obtained explicitly knowing the solution u at time t. Fortunately it can be expressed as a tridiagonal system of equations that we solved numerically.

    Stochastic Gillespie algorithm

    The basic idea lying here is to compute random times for a population rather than computing them for each individual (which is much more computer intensive). Those random times appear in the recombinations and the birth-death calculus; Gillespie's idea, which we are about to describe, allows to compute the time of the next recombination (resp. birth or death) within the population instead of making the calculus and storing the result for each cell in the device.
    Under the assumption that cells' recombination times are independent from one bacteria to the other, and that these times follow an exponential law of parameter λ, Gillespie tells us that the next recombination time (resp. birth or death time) follows an exponential law of parameter Npop.

    Stochastic movement of the cells

    As promised before, we give the two functions describing the movement of the bacteria in the microfluidic device:

    bin_distrib

    bin_distrib

    where Nsw/unsw denotes either the number of switched bacteria (or recombinated bacteria) or unswitched bacteria. The idea here is to have, on average at each step time, dt.kdiff.f(x,t) bacteria that will move, a result we would get from usual diffusion. As previously mentioned, we cannot use regular diffusion since our algorithm must apply to discrete data lying on a wide range of integers.

    Simulations

    Now let's present the graphical display of our simulation tool. Once the user has selected a refreshment speed x, the software will draw every x minutes the unswitched bacteria in grey, the switched ones in red and the arabinose in small blue points.

    soft1


    We see on the following picture that after a certain amount of time, both cells and arabinose are distributed rather uniformly, and that the concentration of arabinose decreases near the tunnel as it is washed away by the flow.

    soft2


    The software also gives us the dynamics of recombined bacteria, LuxI and AHL slice by slice. This is what we will exploit here.
    In the following simulations, the device was divided into 15 slices. We exhibit some results on the fifth slice. To count, one can either look at the number of recombined bacteria thanks to RFP, or set up the timer thanks to AHL. We will see that depending on the parameters, one or other may not be possible. For the first one, it is necessary to have clear steps, implying that the time between the pulses will be large enough to allow the system to stabilise. For the second, distinct steps on the AHL graph are necessary to set up the timer. Indeed, this allows to design the bacteria in the way that they will produce GFP only once a certain treshold in AHL was reached. However, admitting that this timer is set up, then it is possible to count more frequent events even without steps as we shall see in figure 4. In the sequel, the time between two pulses will be denoted T2.

    AHL_chart
    Figure 2: Switched bact., AHL, LuxI for T2 = 160 and kdiff = 8.10-2


    In the following figure, we show that it is still possible to count the number of pulses thanks to the number of recombined bacteria, although it might become uncertain. Nevertheless, it seems that we can still choose a threshold of AHL beyond which bacteria should be designed to luminesce. For instance, if we wanted our timer to stop at 4, then a threshold of 0.70 in standardised unit seems reasonable.
    AHL_chart3
    Figure 3: Switched bact., AHL, LuxI for T2 = 120 and kdiff = 2.10-2


    Although we have been able to count six pulses thanks to the previous graph of AHL, it is not always possible to do so. For instance, if we shorten the time between two pulses, decreasing the diffusion coefficient, then it is no longer possible to distinguish pulses:
    AHL_chart2
    Figure 4: AHL for T2 = 100 and kdiff = 10-3


    Achievements of our model

    Our model gave us a better understanding of the dynamics and helped us explain and suggest experiments a posteriori. Together with the results of the experiments, we may be able:
    • to have a better determination of some parameters such as caracteristic times for recombination and their statistical distribution, which we can then plug in;

    • to set up parameters in a way that allow us to count with RFP, or to make a timer thanks to AHL and GFP. Those parameters are for instance the speed of the flow in the microfluidic device, the dimensions of the microfluidic device, the initial number of bacteria or the viscosity of the medium fluid.


    Sources of the simulation software

    You may download the sources by clicking here. It is written in C and all the interesting parameters are in the file "constants.h". The files "nrutil.c" and "nrutil.h" are taken from the internet.

    References

    • T.Danino, O.Mondragon-Palomino, L.Tsimring, J.Hasty: A synchronized quorum of genetic clocks. Nature 2009, 463:326 - 330

    • T.Lu, D.Volfson, L.Tsimring, J. Hasty: Cellular growth and division in the Gillespie algorithm. Systems biology, IEEE proceedings 2004, 1:121 - 128

    • D.Gillespie: A general method for numerically simulating the stochastic time evolution of coupled chemical reactions. The Journal of Computational Physics 1976, 22:403 - 434

    • D.Gillespie: Exact stochastic simulation of coupled chemical reactions. The Journal of Physical Chemistry 1977, 81:2340 - 2361

    • S.Lampoudi, D.Gillespie, L.Petzold: The multinomial simulation algorithm for discrete stochastic simulation of reaction-diffusion systems. 2009

    • R.Bustin, H.Messer: An equivalent Markov model for Gillespie’s stochastic simulation algorithm for biochemical systems. 2006

    • W.H.Press, S.A.Teukolsky, W.T.Vetterling, B.P.Flannery: Numerical recipes in C. Cambridge University Press, 1992.