Team:TU Delft-Leiden/Modeling/Curli

From 2014.igem.org

(Difference between revisions)

Latest revision as of 22:56, 17 October 2014

Curli Module

The goal of our project for the conductive curli module is to produce a biosensor that consists of E. coli that are able to build a conductive biofilm, induced by any promoter, see the gadget section of our wiki and the extracellular electron transport (EET) module. The biofilm consists of curli containing His-tags that can connect to gold nanoparticles, see the conductive curli module. When the curli density is sufficiently high, a dense network of connected curli fibrils is present around the cells. Further increasing the amount of curli results in a conductive pathway connecting the cells, thereby forming conductive clusters. Increasing the amount of curli even further, sufficiently curli fibrils are present to have a cluster that connects the two electrodes and thus have a conducting system.
The goal of the modeling of the curli module is to prove that our biosensor system works as expected and to capture the dynamics of our system. So, we want to answer the question: "Does a conductive path between the two electrodes arise at a certain point in time and at which time does this happen?" However, we not only want to answer the question if our system works as expected qualitatively, but we also want to make quantitative predictions about the resistance between the two electrodes of our system in time.

Curli Module

To capture the dynamics of our system, we have implemented a three-layered model, consisting of the gene level layer, the cell level layer and the colony level layer:

At the gene level, we calculate the curli subunits production rates and curli subunit growth that will be used in the cell level.
At the cell level, we use these production and growth rates to calculate the curli growth in time, which we will use at the colony level.
At the colony level layer, we determine if our system works as expected, ie. determine if a conductive path between the two electrodes arises at a certain point in time and at which time this happens. We also determine the change of the resistance between the two electrodes of our system in time.

A figure of our three-layered model is displayed below.

Click in the figure to move to the corresponding page.

Figure 1: A schematic view of our model, which is a three-layered model. Each layer determines characteristic parameters for the layer above it. At the gene level, we calculate the curli subunits production rates. At the cell level, we use these production rates to calculate the curli growth in time. At the colony level, we use the curli growth in time to determine the change of the resistance between the two electrodes of our system in time.

We start with the modeling of the gene expression of proteins involved in the curli formation pathway at the gene level. In the constructs we made in the wet lab, CsgA is continuously being produced and the CsgB gene is placed under the control of a landmine promoter, activated by either TNT or DNT, see the Landmine Detection Module. So, when the cells get induced by TNT or DNT, CsgB protein production will get started and CsgA will already be present in the system, as CsgA is continuously being produced. We first modeled this system by constructing an extensive gene expression model of the curli formation pathway. Subsequently, we simplified this model, so less parameters were needed.
Based on the simplified model, we made a plot of the curli growth as function of time for different initial concentrations of $CsgA_{free}$, see figure 2. We conclude the following from this figure:

Firstly, as expected, curli growth stabilizes to a rate equal to $p_{A}$ after approximately 2 hours, independent of the initial concentration of $CsgA_{free}$. The width of this peak is determined by the product $ k p_B$, where k is the production rate of curli and $ p_B$ is the production rate of CsgB proteins.
Secondly, increasing the initial concentration of $CsgA_{free}$, increases the height of the peak. Even with zero initial $CsgA_{free}$ concentration, a small peak can be found at one hour. This is a consequence of $CsgA_{free}$ build-up when the CsgB concentration is still very small.
Thirdly, during the first two hours, few CsgB proteins are present in the system. We therefore expect that the length of the curli fibrils that started in the first few hours are much longer than the fibrils that started at later times.

Figure 2: The curli subunit growth in units per second for various initial concentrations $ A_0 $ of CsgA as function of time. Initial concentrations that equal 0, 5, 10 or 15 hours of CsgA production are shown.

Using the growth rate of curli and production of CsgB protein as function of time obtained from the gene level model, the conductance as a function of time can be computed for the cell. We obtained an analytical expression for $\rho_{curli}$, which represents the density of curli fibrils around the cell. We have fitted the function $$ \rho_n = C_{1_n} e^{-\frac{r}{C_{2_n}}} + C_{3_n} e^{-\frac{r}{C_{4_n}}} \tag{1}$$ to our curli density curves at each time $ n $, see figure 3 the red line. At first, we tried to fit our data to only one exponential term (green line). It can clearly be seen in the figure that this does not adequately capture the dynamics of the curve. However, equation 1 gives a very good fit to the curli density curves at each time $ n $. The reason for fitting such a simple function is that, in the colony level, we need to quantify the conductance between the cells. The integral for this rather complicated and we need an analytical function for $\rho_{curli}$ to analytically solve this integral.
We also calculated the conductive radius of the cell as a function of the radius, see figure 4. The conductive radius is the largest radius where $\rho_{curli}$ is bigger than a certain threshold of curli density. We use the conductive radius in the colony level to determine when a conductive path between the two electrodes of our system arises.

Figure 3: Blue line: Right behind the red line, at t=5 hr the mean of all density curves. Green line: a weighted fit of $ \rho_n = C_{1_n} e^{-\frac{r}{C_{2_n}}} $. Red line: A fit $ \rho_n = C_{1_n} e^{-\frac{r}{C_{2_n}}} + C_{3_n} e^{-\frac{r}{C_{4_n}}} $ to the blue line.

Figure 4: The green lines are the conductive radius plotted versus the time for 100 cells with a critical density of $ \rho_{crit}=1204 $ curli subuntis $ \mu m ^{-3} $. The orange red represents the mean conductive radius and the dark blue lines represent two standard deviations from the mean.

An illustrative view of what our cell looks like during the adding of curli subunits is shown in figure 5.

meh — Figure 5: Schematic view of our cell (black sphere centred at x=y=z=0) with growing curli fibrils. The wires represent the curli fibrils. **Click to play!**

meh — Figure 5: Schematic view of our cell (black sphere centred at x=y=z=0) with growing curli fibrils. The wires represent the curli fibrils. **Click to play!**

Now that we determined values for $\rho_{curli}$ and $r_{cond}$ at the cell level, we can finally predict if our system works as expected and capture the dynamics of our system. Firstly, we want to prove that our system works as expected. So, we want to predict if a conductive path between the two electrodes arise at a certain point in time and at which time this happens. Secondly, we not only want to answer the question if our system works as expected with a yes or no answer, but we also want to make quantitative predictions about the resistance between the two electrodes of our system in time. We do this by modeling the curli growth on the colony level; each cell is now visualized and has curli growth. Now, we have come up with two different approaches:

Firstly, we let the cells increase their conductive radius in time, according with our findings on the cellular level (figure 4). A connection is created from one electrode to the other electrode when there is a conductive path between them. Conductive paths consists of cells that have a connection between each other; cells connect when there is an overlap between their conductive radius. This problem is very similar to problems in percolation theory. From this, we can make conclusions about how our system works in an experimental setting.
Secondly, we also want to make quantitative predictions about the resistance between the two electrodes of our system in time. Therefore, we used graph theory to translate the cells on the chip to a graph and used an algorithm from graph theory to calculate the resistance between the two electrodes. The conductance between the cells is computed from an integral and is equal to: $$ \sigma (y) = \ 2 \pi \left( C_{1}^2 C_2 e^{-\frac{y}{C_2}} + C_{3}^2 C_4 e^{-\frac{y}{C_4}} + \frac{4 C_{1}C_3 C_2^2 C_4^2}{y \left( C_2^2 - C_4^2 \right)} \left( e^{-\frac{y}{C_2}} -e^{-\frac{y}{C_4}} \right) \right) \tag{1}$$

Using the first approach, a simulation of our resulting model is shown in figure 6. Percolation is computed by applying an algorithm that can find clusters of connected cells. When one of the clusters connects both electrodes, there is percolation.

meh — Figure 6: NorthWest: A visual representation of our cells on the plate. The circles represent the cells with an increasing conductive radius. In this simulation there are 5000 cells present on a chip of 500µmx500µm. NorthEast: A spy matrix of 5000x5000 where the blue dots represent connections between the individual cells. A blue dot on position x,y means that cell x is connected with y. Each cell is connected to itself (diagonal). At the point of percolation, $ \approx 0.1 \% $ of the matrix is connected, meaning that each cell is on average connected to 5 others. SouthWest: Each square of nxn represents a cluster of n connected cells. The squares are sorted from small to large. SouthEast: This figure shows the largest cluster of cells in different colors. **Click to play!**

meh — Figure 6: NorthWest: A visual representation of our cells on the plate. The circles represent the cells with an increasing conductive radius. In this simulation there are 5000 cells present on a chip of 500µmx500µm. NorthEast: A spy matrix of 5000x5000 where the blue dots represent connections between the individual cells. A blue dot on position x,y means that cell x is connected with y. Each cell is connected to itself (diagonal). At the point of percolation, $ \approx 0.1 \% $ of the matrix is connected, meaning that each cell is on average connected to 5 others. SouthWest: Each square of nxn represents a cluster of n connected cells. The squares are sorted from small to large. SouthEast: This figure shows the largest cluster of cells in different colors. **Click to play!**

We simulated our resulting model 100 times and for each point in time we checked the chance of percolation without variation in $r_{cond}$ (yellow line in figure 7) and with cellular variation in $r_{cond}$ (blue line in figure 7). Fortunately, the blue and yellow lines are very similar. This means that cellular variation has little influence on the chance of percolation at each point in time. Therefore, the results of our model are robust to cellular variation and it is likely that many factors that could increase the cellular variation, e.g. different CsgA or CsgB protein production rates, are relatively unimportant.

Figure 7: The chance of percolation with 5000 cells on a 500x500 $\mu m $ chip. as function of time. The results are from 100 simulations. The yellow line represents the chance of percolation where all the cells have the same conductive radius. The blue line is the same simulation, but all cells have slightly different conductive radii. Note how there is no notable difference between the two.

As we want our system to be usable as a biosensor, it has to be strongly dependent on the analyte concentration, and therefore the CsgB production rate. From figure 8, we see that there are very distinguishable difference for different CsgB production rates. First of all, the moment of percolation differs a lot. Equally important, the transition from no percolation to percolation is much less sharper.

Figure 8: The change of induction for t=0:10 hours, when cellular differences are included in the cell level for different induction strengths. The orange line is created by reducing the promoter strength of the cyan line ($p_B = 1.3 \cdot 10^{-13} M/s $ ) by 50%.

Using the second approach, we have calculated the conductance as function of time for different dimensions of our plate. We varied the distance between the electrodes between 210, 460 and 960 $ \mu m $. The length of the electrodes is kept at $ 500 \mu m$. The result is shown in figure 9.

Figure 9: The conductance of our system as function of time, using the results from the gene and the colony level. Different simulations have been done with different chip sizes. The length of the chips is in all cases $ 500 \mu m $. The distance between the electrodes is varied from $ 210 \mu m $ (green lines) to $ 460 \mu m $ (red lines) and $ 960 \mu m $ (blue lines). There are multiple lines for the green and red lines as we repeated the simulations multiple times. However the simulation for when the electrodes were $ 960 \mu m $ apart was only performed once.

Another parameter we investigated is the influence of cell density on the conductivity of the system. From percolation theory we have learned that this has a strong influence on the moment of percolation. The results are shown in figure 10.

Figure 10: The conductance of our chip ( $ 500\mu m \ by \ 500\mu m $ as function of time. Two different cell densities have been simulated. The red lines show simulations with a cell density of $ 2 \cdot 10^4 cells / mm^2$. The blue lines show the results with a cell density of $ 2.4 \cdot 10^4 cells / mm^2$

The most important question of all is: Will our system work, and under what conditions? To answer this question we have done the exact same simulation for the gene level and the cell level as elsewhere, except with the only difference that we halved the promoter strength of CsgB. This simulates less CsgB protein production, which corresponds to less TNT/DNT present in the system. The results are shown in figure 11.

Figure 11: The conductance of our chip ( $ 500\mu m \ by \ 500\mu m $ as function of time. Two different activation strengths have been simulated. The red lines (same as in figure 6 and 7) show a 100% activated CsgB promoter. The blue lines show the results for a promoter that is activated for 50%.

Based on the above three figures, there are a couple of recommendations we have for the design of the chip:

Based on figure 9, we recommend to decrease the distance between the electrodes. This will increase the conductance without increasing the uncertainty. From the shape of the curve, we expect that we do not have to wait very long (only a couple of hours) to do a measurement. After a steep increase in conductance, the increase in conductance will be relatively small in time. Furthermore, the uncertainty due to the placement of our cells increases in time.
Based on figure 10, we recommend high cell density for increased signal strength of our system. However, it is crucial that the cell density is kept constant throughout the measurement. Cell density influences the conductance significantly and when the cell density is not kept constant, it will influence the conductance measurements in such a way that an increase or decrease in conductance can not only be attributed to an increase or decrease in DNT any longer.
From figure 11, we learned that we have to be very careful with our measurement time in order to be able to draw conclusions about the concentration of TNT/DNT. The influence of the induction strength of the CsgB promoter is very small and its influence decreases after the initial amount of CsgA protein have formed curli. A recommendation for the wet lab is to increase the waiting time prior to induction of the promoter and thereby increasing the initial amount of CsgA protein. This will increase the relative differences in conductance between different induction strengths of the CsgB promoter. In our model, this increase in relative differences in conductance happens between 0 and 2.5 hours.
The time this increase in relative differences in conductance happens, is strongly dependent on the parameter k in our simplified gene model. For our model, we could fit this parameter to conductance curves obtained in the wet lab. Measuring conductance until observing a linear increase in time will give the most useful information for obtaining the inducer concentration, as the time of this transition depends on the strength of induction.

@@ Line 3: / Line 3: @@
 {{:Team:TU_Delft-Leiden/Templates/Start}}
-{{:Team:TU_Delft-Leiden/Templates/stylemod}}
 <html>
 <body>
 <div class='grid_12'>
-<h2> Curli Model</h2>
-</body>
+<h2> Curli Module</h2>
-</html>
-{{:Team:TU_Delft-Leiden/Templates/tocmodnonfloat}}
+<p>
+The goal of our project for the conductive curli module is to produce a biosensor that consists of <i>E. coli</i> that are able to build a conductive biofilm, induced by any promoter, see the <a href="https://2014.igem.org/Team:TU_Delft-Leiden/Project/Gadget">gadget section</a> of our wiki and the <a href="https://2014.igem.org/Team:TU_Delft-Leiden/Project/Life_science/EET">extracellular electron transport (EET) module</a>. The biofilm consists of curli containing His-tags that can connect to gold nanoparticles, see the <a href="https://2014.igem.org/Team:TU_Delft-Leiden/Project/Life_science/curli">conductive curli module</a>. When the curli density is sufficiently high, a dense network of connected curli fibrils is present around the cells. Further increasing the amount of curli results in a conductive pathway connecting the cells, thereby forming conductive clusters. Increasing the amount of curli even further, sufficiently curli fibrils are present to have a cluster that connects the two electrodes and thus have a conducting system. <br>
+The goal of the modeling of the curli module is to prove that our biosensor system works as expected and to capture the dynamics of our system. So, we want to answer the question: "Does a conductive path between the two electrodes arise at a certain point in time and at which time does this happen?" However, we not only want to answer the question if our system works as expected qualitatively, but we also want to make quantitative predictions about the resistance between the two electrodes of our system in time.
+</p>
-<html>
-<body>
-<a name="Gene Level"></a>
-<h3> Gene Level Modeling </h3>
-<p>We will start with modelling the expression of Curli on the gene level. Proteins that are dedicated to curli formation are the CsgA/B/D/E/F/G [1]. CsgA is the main building block of the curli. When produced this is secreted out of the cell by the csgEFG complex. In the absence of CsgB there is no curli formation, since the CsgA proteins remain unpolymerized. CsgB is the starting block of the curli fibrils, likes piles under a house, that connect the cell membrane to the first CsgA. Once this is placed on the cell surfaced, the CsgA can polymerize onto the starting curli fibril. In the constructs we produce, CsgA is continuously produced. The CsgB promoter, on the other hand is in our design activated by a landmine promoter, activated by either TNT or DNT. In our design, the curli are first grown in a medium. During this time, they produce CsgA. Hence, there is an initial concentration CsgA present in the system when the CsgB promoter is activated. The csgEFG operon encodes four accessory proteins required for the curli assembly. Overexpression of these genes is not needed [bron?] for curli formation. </p>
+<div class="tableofcontents">
+<ul> <p> Curli Module </p>
+               <ul>
+                    <li>
+                    <a href="https://2014.igem.org/Team:TU_Delft-Leiden/Modeling/Curli/Gene">
+                    <p>Gene Level Modeling</p>
+                    </a>
+                          <ul>
+                               <li>
+                               <a href="https://2014.igem.org/Team:TU_Delft-Leiden/Modeling/Curli/Gene#extendedgenelevel">
+                               <p>Extensive Gene Level Modeling</p>
+                               </a>
+                               </li>
+                               <li>
+                               <a href="https://2014.igem.org/Team:TU_Delft-Leiden/Modeling/Curli/Gene#simplifiedgenelevel">
+                               <p>Simplified Gene Level Modeling</p>
+                               </a>
+                               </li>
+                          </ul>
+                    </li>
+                    <li>
+                    <a href="https://2014.igem.org/Team:TU_Delft-Leiden/Modeling/Curli/Cell">
+                    <p>Cell Level Modeling</p>
+                    </a>
+                          <ul>
+                               <li>
+                               <a href="https://2014.igem.org/Team:TU_Delft-Leiden/Modeling/Curli/Cell#discretization">
+                               <p>Discretization of Gene Level Model</p>
+                               </a>
+                               </li>
+                               <li>
+                               <a href="https://2014.igem.org/Team:TU_Delft-Leiden/Modeling/Curli/Cell#building">
+                               <p>Building the Curli Fibrils</p>
+                               </a>
+                               </li>
+                               <li>
+                               <a href="https://2014.igem.org/Team:TU_Delft-Leiden/Modeling/Curli/Cell#FittingCurliDensity">
+                               <p>Fitting the Curli Density</p>
+                               </a>
+                               </li>
+                               <li>
+                               <a href="https://2014.igem.org/Team:TU_Delft-Leiden/Modeling/Curli/Cell#conductiveradius">
+                               <p>Conductive Radius of the Cell</p>
+                               </a>
+                               </li>
+                          </ul>
+                    </li>
+                    <li>
+                    <a href="https://2014.igem.org/Team:TU_Delft-Leiden/Modeling/Curli/Colony">
+                    <p>Colony Level Modeling</p>
+                    </a>
+                          <ul>
+                               <li>
+                               <a href="https://2014.igem.org/Team:TU_Delft-Leiden/Modeling/Curli/Colony#percolation">
+                               <p>Percolation</p>
+                               </a>
+                               </li>
+                               <li>
+                               <a href="https://2014.igem.org/Team:TU_Delft-Leiden/Modeling/Curli/Colony#resistivity">
+                               <p>Resistance</p>
+                               </a>
+                               </li>
+                               <li>
+                               <a href="https://2014.igem.org/Team:TU_Delft-Leiden/Modeling/Curli/Colony#recommendations">
+                               <p>Recommendations for product design and wet lab</p>
+                               </a>
+                               </li>
+                          </ul>
+                               <li>
+                               <a href="https://2014.igem.org/Team:TU_Delft-Leiden/Modeling/Curli/Reflection">
+                               <p>Critical Reflection on our Model</p>
+                               </a>
+                               </li>
+                    </li>
+               </ul>
+</ul>
+</div>
 <br>
-<h4> Extended gene level modeling</h4>
-<p>to be written</p>
+<p>
+To capture the dynamics of our system, we have implemented a three-layered model, consisting of the gene level layer, the cell level layer and the colony level layer:
+<ul>
+<li>
+At the gene level, we calculate the curli subunits production rates and curli subunit growth that will be used in the cell level.
+</li>
+<li>
+At the cell level, we use these production and growth rates to calculate the curli growth in time, which we will use at the colony level.
+</li>
+<li>
+At the colony level layer, we determine if our system works as expected, ie. determine if a conductive path between the two electrodes arises at a certain point in time and at which time this happens. We also determine the change of the resistance between the two electrodes of our system in time.
+</li>
+</ul>
+</p>
+<p>
+A figure of our three-layered model is displayed below.
+</p>
+<h3>Click in the figure to move to the corresponding page.</h3>
+<figure>
+<img src="https://static.igem.org/mediawiki/2014/f/fa/TU_Delft_2014_ModellingFigureCurli_small.png" width="800" height="920" usemap="#ModelMap">
+<map name="ModelMap">
+  <area shape="rect" coords="0,0,800,296" href="https://2014.igem.org/Team:TU_Delft-Leiden/Modeling/Curli/Gene" alt="Gene level">
+  <area shape="rect" coords="0,296,800,592" href="https://2014.igem.org/Team:TU_Delft-Leiden/Modeling/Curli/Cell" alt="Cell level">
+  <area shape="rect" coords="0,592,800,920" href="https://2014.igem.org/Team:TU_Delft-Leiden/Modeling/Curli/Colony" alt="Colony level">
+</map>
+<figcaption>
+Figure 1: A schematic view of our model, which is a three-layered model. Each layer determines characteristic parameters for the layer above it. At the gene level, we calculate the curli subunits production rates. At the cell level, we use these production rates to calculate the curli growth in time. At the colony level, we use the curli growth in time to determine the change of the resistance between the two electrodes of our system in time.
+</figcaption>
+</figure>
 <br>
-<h4> Analytical gene level modeling</h4>
-<p>Though the model as described above, providing that all rates are known, has a more accurate (though still simplified) representation of the curli assembly system, we have chosen decrease the complexity further to the bare essentials. Most of the production rates and kinetic constants cannot be found in literature. Measuring the accurate rates is, within the scope of this project, unfeasible for many reasons. First of all, we think that the secretion of the proteins is not a rate limiting step. We are interested in the CsgA/B production and curli formation. The scope in which we expect the system to respond is in the order of hours. If secretion were the rate limiting step, it would mean that many CsgA/B units are piled up inside the bacteria. One might raise the question whether it is justified that the growth of new subunits is independent of the length of the curli, since the CsgA is secreted at the edge of the cell distance and has to ‘travel’ towards the end of the curli fibrils. This would suggest that the assumption used in deterministic modelling that all concentrations are homogeneous is not justified. A quick calculation shows that after one second, the displacement of a spherical particle with radius r=10nm due to Brownian motion in liquid water at room temperature using equation 0 is 6.6 μm; many times the bacterial radius! Hence the diffusion is not rate limiting. [4] </p>
-$$ \bar{x}^2  = \frac{k_b T t}{3 \pi \eta r} \tag{0}$$
+<p>
-<p> What we do expect to be rate limiting for biofilm formation is that a large amount of curli subunits has to be created, hence the expression rate of the proteins. Instead of including the intermediate steps, we have replaced the production rate of the CsgA and CsgB with one expression rate each. These rates have to be measured in the lab. What we did not include is the maturation time needed for the production. Luckily, for CsgB speeds up the CsgA formation, the lag phase of CsgA is decreased to the same lag phase as CsgB when they’re both expressed. [2][3] </p>
+We start with the modeling of the gene expression of proteins involved in the curli formation pathway at the gene level. In the constructs we made in the wet lab, CsgA is continuously being produced and the CsgB gene is placed under the control of a landmine promoter, activated by either TNT or DNT, see the <a href="https://2014.igem.org/Team:TU_Delft-Leiden/WetLab/landmine">Landmine Detection Module</a>. So, when the cells get induced by TNT or DNT, CsgB protein production will get started and CsgA will already be present in the system, as CsgA is continuously being produced. We first modeled this system by constructing an extensive gene expression model of the curli formation pathway. Subsequently, we simplified this model, so less parameters were needed. <br>
+Based on the simplified model, we made a plot of the curli growth as function of time for different initial concentrations of \(CsgA_{free}\), see figure 2. We conclude the following from this figure: </p>
+<ul>
+<li>
+Firstly, as expected, curli growth stabilizes to a rate equal to \(p_{A}\) after approximately 2 hours, independent of the initial concentration of \(CsgA_{free}\). The width of this peak is determined by the product \( k p_B\), where k is the production rate of curli and \( p_B\) is the production rate of CsgB proteins.
+</li>
+<li>
+Secondly, increasing the initial concentration of \(CsgA_{free}\), increases the height of the peak. Even with zero initial \(CsgA_{free}\) concentration, a small peak can be found at one hour. This is a consequence of \(CsgA_{free}\) build-up when the CsgB concentration is still very small.
+</li>
+<li>
+Thirdly, during the first two hours, few CsgB proteins are present in the system. We therefore expect that the length of the curli fibrils that started in the first few hours are much longer than the fibrils that started at later times.
+</li>
+</ul>
+<figure>
+<img src="https://static.igem.org/mediawiki/2014/e/e9/Delft2014_DifferentA0.png" width="60%" height="60%">
+<figcaption>
+Figure 2: The curli subunit growth in units per second for various initial concentrations \( A_0 \) of CsgA as function of time. Initial concentrations that equal 0, 5, 10 or 15 hours of CsgA production are shown.
+</figcaption>
+</figure>
 <br>
-<p> We will use the following system:</p>
-$$ \emptyset \xrightarrow{p_{A}} \ CsgA \tag{1} $$
-$$ \emptyset \xrightarrow{p_{B}} \ CsgB \tag{2} $$
-$$ CsgA + CsgB \xrightarrow{K} C_{sub} + CsgB \tag{3} $$
-<p>\(C_{sub}\) is a folded CsgA protein that has connected to one of the curli fibrils. \(P_{A}\) and \(P_{B}\) are the production rates of CsgA and CsgB respectively. Equation (3) may seem confusing, since in reality the growing of curli takes place between a \(C_{sub}\) and a CsgA. However, this happens only at the end of the curli fibrils. In deterministic modelling, the system is a homogeneous batch and we cannot discriminate between the outer and the other curli subunits. It is theoretically possible to model the system, as an infinite amount of possible curli reactions that can take place that increase a curli with length i to length i+1 at rate K. [7] However, we are merely interested in the growth rates of the curli, since the distribution of the curli length will follow from the model of <a href="url">the cell level</a>. The reactions will only take place at the end points of the curli fibrils. The reaction rate is therefore proportional to the fibril concentration, hence the amount of CsgB units. </p>
+<p>
+Using the growth rate of curli and production of CsgB protein as function of time obtained from the gene level model, the conductance as a function of time can be computed for the cell. We obtained an analytical expression for \(\rho_{curli}\), which represents the density of curli fibrils around the cell. We have fitted the function
+$$ \rho_n = C_{1_n} e^{-\frac{r}{C_{2_n}}} + C_{3_n} e^{-\frac{r}{C_{4_n}}} \tag{1}$$
+ to our curli density curves at each time \( n \), see figure 3 the red line. At first, we tried to fit our data to only one exponential term (green line). It can clearly be seen in the figure that this does not adequately capture the dynamics of the curve. However, equation 1 gives a very good fit to the curli density curves at each time \( n \). The reason for fitting such a simple function is that, in the colony level, we need to quantify the conductance between the cells. The integral for this rather complicated and we need an analytical function for \(\rho_{curli}\) to analytically solve this integral. <br>
+We also calculated the conductive radius of the cell as a function of the radius, see figure 4. The conductive radius is the largest radius where \(\rho_{curli}\) is bigger than a certain threshold of curli density. We use the conductive radius in the colony level to determine when a conductive path between the two electrodes of our system arises.
+</p>
+<figure>
+<img src="https://static.igem.org/mediawiki/2014/3/3f/TUDelft_2014_FittingCurli.png" width="60%" height="60%">
+<figcaption>
+Figure 3: Blue line: Right behind the red line, at t=5 hr the mean of all density curves. Green line: a weighted fit of \( \rho_n = C_{1_n} e^{-\frac{r}{C_{2_n}}} \). Red line: A fit \( \rho_n = C_{1_n} e^{-\frac{r}{C_{2_n}}} + C_{3_n} e^{-\frac{r}{C_{4_n}}} \) to the blue line.
+</figcaption>
+</figure>
 <br>
-<p> Writing equations 1-3 into differential equations results in:</p>
-$$ \frac{d}{dt} [CsgA] = p_{A} - K [CsgA][CsgB] \tag{4.1} $$
-$$ \frac{d}{dt} [CsgB] = p_{B} \tag{4.2} $$
-$$ \frac{d}{dt} [C_{sub}] = K [CsgA][CsgB] \tag{4.3} $$
-<p> Fortunately, this system can be solved analytically. To do this, we need the initial conditions. Say the CsgB promoter is activated at \(t=0\) (or maybe at \(t=-T_{lagtime}\) if you want to include this as well). At this time there are no fibrils present, so \([CsgB]|_{t=0} =[Csub]|_{t=0}=0\). However, the CsgA promoter is continuously active, so we expect that at time \(t_{0}\), there is an initial concentration \(A_0\) of CsgA.</p>
+<figure>
+<img src="https://static.igem.org/mediawiki/2014/e/e1/TUDelft_2014_cellradius.png" width="60%" height="60%">
+<figcaption>
+Figure 4: The green lines are the conductive radius plotted versus the time for 100 cells with a critical density of \( \rho_{crit}=1204 \) curli subuntis \( \mu m ^{-3} \).  The orange red represents the mean conductive radius and the dark blue lines represent two standard deviations from the mean.
+</figcaption>
+</figure>
 <br>
-<p> The solution to equation 4.2 is trivial: </p>
-$$ [CsgB] = p_B t \tag{5}$$
-<p>Substituting this into equation 4.1 results in:</p>
-$$ \frac{d}{dt} [CsgA] = p_{A} - K p_B [CsgA]t \tag{6} $$
-<p>It can easily be proven that a first order differential equation of the form \(y(t)'+f(t)y(t)=g(t)\)   has solution of the form \(y(t)=e^{-F(t) }\int{g(t) e^{F(t)}  dt}+y_0 e^{-F(t) }\), where \(F(t)= \int{f(t) dt}\). Replacing \(f(t)\) with \(K p_B t\) and \(g(t)\) with \(p_A\) yields equation 7 </p>
-$$ [csgA]= p_A e^{\frac{-Kp_B t^2}{2}}\int{e^{\frac{Kp_B t^2}{2}  } dt} + C e^{\frac{-K p_B t^2}{2}}=p_A e^{\frac{-Kp_B t^2}{2}}\int_{\tau=0}^{t}{e^{\frac{Kp_B \tau^2}{2} }  d\tau}+\tilde{C} ̃e^{\frac{-Kp_B t^2}{2}  } \tag{7} $$
-<p>One with a keen eye may recognize the Dawson function (equation 8): </p>
-$$ D_+ (x)= e^{-x^2 } \int_{y=0}^x{e^{y^2} dy} \tag{8}  $$
+<p>
+An illustrative view of what our cell looks like during the adding of curli subunits is shown in figure 5.
+</p>
-<p>Substituting \(x^2\) and \(y^2\) with \(K p_B t^2\) and \(K p_B \tau^2\), equation 9 obtained.</p>
-$$ [csgA]=\frac{p_A D_+ (t\sqrt{\frac{K p_B}{2}})}{\sqrt{\frac{K p_B}{2}}}  + \tilde{C} e^{\frac{-K p_B t^2}{2 } }  \tag{9}$$
-<p>Using the boundary condition \([CsgB]|_{t=0}=A_0\), the concentration becomes:</p>
+<figure>
-$$ [csgA]=\frac{p_A D_+ (t\sqrt{\frac{K p_B}{2}})}{\sqrt{\frac{K p_B}{2}}}  + A_0 e^{\frac{-K p_B t^2}{2 } }  \tag{10}$$
+<img class="theimage" src="https://static.igem.org/mediawiki/2014/e/e3/TU_Delft_2014_Curli_picturestill.png" alt="meh" width="60%" height="60%"/>
+<img class="theimage2 hidden" src="https://static.igem.org/mediawiki/2014/3/37/TU_Delft_2014_Curli_picture.gif" alt="meh" width="60%" height="60%"/>
+<figcaption>
+Figure 5: Schematic view of our cell (black sphere centred at x=y=z=0) with growing curli fibrils. The wires represent the curli fibrils. <b> Click to play!</b>
+</figcaption>
+</figure>
+<script>
+$(".theimage").click(function() {
+$(this).addClass('hidden');
+$(this).next().removeClass('hidden');
+});
+$(".theimage2").click(function() {
+  $(this).addClass('hidden');
+$(this).prev().removeClass('hidden');
+});
+</script>
+<br>
-<p> Since we’re mainly interested in the curli production rate, rather than the total amount of curli, we can fill in equations 10 and 5 into equation 4.3, which gives us equation 11. </p>
-$$ \frac{d}{dt} [C_{sub}] = K p_B t \frac{p_A D_+ (t\sqrt{\frac{K p_B}{2}})}{\sqrt{\frac{K p_B}{2}}}  + A_0 e^{\frac{-K p_B t^2}{2 } }  \tag{11} $$
+<p>
+Now that we determined values for \(\rho_{curli}\) and \(r_{cond}\) at the cell level, we can finally predict if our system works as expected and capture the dynamics of our system. Firstly, we want to prove that our system works as expected. So, we want to predict if a conductive path between the two electrodes arise at a certain point in time and at which time this happens. Secondly, we not only want to answer the question if our system works as expected with a yes or no answer, but we also want to make quantitative predictions about the resistance between the two electrodes of our system in time. We do this by modeling the curli growth on the colony level; each cell is now visualized and has curli growth. Now, we have come up with two different approaches:<ul>
+<li>
+Firstly, we let the cells increase their conductive radius in time, according with our findings on the cellular level (figure 4). A connection is created from one electrode to the other electrode when there is a conductive path between them. Conductive paths consists of cells that have a connection between each other; cells connect when there is an overlap between their conductive radius. This problem is very similar to problems in <a href="https://2014.igem.org/Team:TU_Delft-Leiden/Modeling/Techniques#PercolationTheory">percolation theory</a>. From this, we can make conclusions about how our system works in an experimental setting.
+</li>
+<li>
+Secondly, we also want to make quantitative predictions about the resistance between the two electrodes of our system in time. Therefore, we used <a href="https://2014.igem.org/Team:TU_Delft-Leiden/Modeling/Techniques#GraphTheory">graph theory</a> to translate the cells on the chip to a graph and used an algorithm from graph theory to calculate the resistance between the two electrodes. The conductance between the cells is computed from an integral and is equal to:
+$$  \sigma (y) = \ 2 \pi \left( C_{1}^2 C_2 e^{-\frac{y}{C_2}} + C_{3}^2 C_4 e^{-\frac{y}{C_4}} +  \frac{4 C_{1}C_3 C_2^2 C_4^2}{y \left( C_2^2 - C_4^2 \right)} \left( e^{-\frac{y}{C_2}} -e^{-\frac{y}{C_4}} \right) \right)  \tag{1}$$
+</li>
+</ul>
+</p>
-<p> For the parameters in equation 11, we have estimated the following values [explain!]: </p>
-</br>
-<table style="width:100%">
-<caption>
-Table 1: Parameters used to obtain quantitative results from the analytical solution for the curli production.
-</caption>
- <tr>
-    <th><b>Parameters</b></th>
-    <th><b>Value</th>
-    <th><b>Unit</th>
- </tr>
- <tr>
-    <td>\(\boldsymbol{P_{A}}\)</td>
-    <td>\(1.0 \cdot 10^{-10}\)</td>
-    <td>\(\frac{1}{Ms}\)</td>
- </tr>
- <tr>
-    <td>\(\boldsymbol{P_{B}}\)</td>
-    <td>\(3.4 \cdot 10^{-12}\)</td>
-    <td>\(\frac{M}{s}\)</td>
- </tr>
- <tr>
-    <td>\(\boldsymbol{K}\)</td>
-    <td>\(4.0 \cdot 10^{4}\)</td>
-    <td>\(\frac{1}{Ms}\)</td>
- </tr>
- <tr>
-    <td>\(\boldsymbol{A_0}\)</td>
-    <td>\(6.0 \cdot 10^{-6}\)</td>
-    <td>\(M\)</td>
- </tr>
-</table>
 <br>
-<p> Plotting equation 11 with the parameter values in table 1 yields the graph shown in figure 1. [insert picture!]</p>
+<p>
+Using the first approach, a simulation of our resulting model is shown in figure 6. Percolation is computed by applying an algorithm that can find clusters of connected cells. When one of the clusters connects both electrodes, there is percolation.
+</p>
+<figure>
+<img class="theimage" src="https://static.igem.org/mediawiki/2014/1/1d/TU_Delft_2014_model_curli_colony_static.png" alt="meh" width="60%" height="60%"/>
+<img class="theimage2 hidden" src="https://static.igem.org/mediawiki/2014/9/9f/TU_Delft_2014_model_curli_colony.gif" alt="meh" width="60%" height="60%"/>
+<figcaption>
+Figure 6: NorthWest: A visual representation of our cells on the plate. The circles represent the cells with an increasing conductive radius. In this simulation there are 5000 cells present on a chip of 500µmx500µm. NorthEast: A spy matrix of 5000x5000 where the blue dots represent connections between the individual cells. A blue dot on position x,y means that cell x is connected with y. Each cell is connected to itself (diagonal). At the point of percolation, \( \approx 0.1 \% \) of the matrix is connected, meaning that each cell is on average connected to 5 others. SouthWest: Each square of nxn represents a cluster of n connected cells. The squares are sorted from small to large. SouthEast: This figure shows the largest cluster of cells in different colors. <b> Click to play!</b>
+</figcaption>
+</figure>
+<script>
+$(".theimage").click(function() {
+$(this).addClass('hidden');
+$(this).next().removeClass('hidden');
+});
+$(".theimage2").click(function() {
+ $(this).addClass('hidden');
+$(this).prev().removeClass('hidden');
+});
+</script>
 <br>
+<p>
+We simulated our resulting model 100 times and for each point in time we checked the chance of percolation without variation in \(r_{cond}\) (yellow line in figure 7) and with cellular variation in \(r_{cond}\) (blue line in figure 7). Fortunately, the blue and yellow lines are very similar. This means that cellular variation has little influence on the chance of percolation at each point in time. Therefore, the results of our model are robust to cellular variation and it is likely that many factors that could increase the cellular variation, e.g. different CsgA or CsgB protein production rates, are relatively unimportant.
+</p>
+<img src="https://static.igem.org/mediawiki/2014/d/dc/TUDelft_2014_500x500ChipResults.png" width="60%" height="60%">
+<figcaption>
+Figure 7: The chance of percolation with 5000 cells on a 500x500 \(\mu m \) chip. as function of time. The results are from 100 simulations. The yellow line represents the chance of percolation where all the cells have the same conductive radius. The blue line is the same simulation, but all cells have slightly different conductive radii. Note how there is no notable difference between the two.
+</figcaption>
+</figure>
+<br>
+<p>
+As we want our system to be usable as a biosensor, it has to be strongly dependent on the analyte concentration, and therefore the CsgB production rate. From figure 8, we see that there are very distinguishable difference for different CsgB production rates. First of all, the moment of percolation differs a lot. Equally important, the transition from no percolation to percolation is much less sharper.
+</p>
 <figure>
-<img src="" width="100%" height="100%">
+<img src="https://static.igem.org/mediawiki/2014/0/02/TUDelft2014_Percolation_InductionDifferences.png" width="60%" height="60%">
 <figcaption>
-Figure 1: hoi!
+Figure 8: The change of induction for t=0:10 hours, when cellular differences are included in the cell level for different induction strengths. The orange line is created by reducing the promoter strength of the cyan line (\(p_B = 1.3 \cdot 10^{-13} M/s \) ) by 50%.
 </figcaption>
 </figure>
 <br>
-<p> Figure 1 shows a steady production of CsgB. The curli subunit production at \(t=0\) is as expected zero, since there is no CsgB at that point. In the next few hours, the curli subunit production peaks. We think that this is due to the high concentration of CsgA that is present at \(t=0\). In figure 2, the curli production as function of time is plotted for different initial concentrations of CsgA.
-<br>
+<p>
+Using the second approach, we have calculated the conductance as function of time for different dimensions of our plate. We varied the distance between the electrodes between 210, 460 and 960 \( \mu m \). The length of the electrodes is kept at \( 500 \mu m\). The result is shown in figure 9.
+</p>
 <figure>
-<img src="" width="100%" height="100%">
+<img src="https://static.igem.org/mediawiki/2014/c/c5/TUDelft2014_Colony_DifferentSizes.png" width="60%" height="60%">
 <figcaption>
-Figure 2: The initial CsgA concentration equals x hours of CsgA production of the bacteria
+Figure 9: The conductance of our system as function of time, using the results from the gene and the colony level. Different simulations have been done with different chip sizes. The length of the chips is in all cases \( 500 \mu m \). The distance between the electrodes is varied from \( 210 \mu m \) (green lines) to \( 460 \mu m \) (red lines) and \( 960 \mu m \) (blue lines). There are multiple lines for the green and red lines as we repeated the simulations multiple times. However the simulation for when the electrodes were \( 960 \mu m \) apart was only performed once.
 </figcaption>
 </figure>
- <br>
-<p>A few conclusions can be drawn regarding the CsgA production. Firstly, as expected, the curli growth stabilizes to a rate equal to pA after approximately 2 hours, independent of the initial concentration \(A_0\). Increasing the concentration \(A_0\) increases the height of the peak. During the first two hours, few CsgB subunits are present. We therefore expect that the length of the curli fibrils that started in the first few hours are much longer than fibrils that started at later times. Even with zero initial curli concentration, a small peak can be found at one hour. This is a consequence of the curli build-up near \(t=0\) when the CsgB concentration is very small.</p>
 <br>
-<a name="Cell Level"></a>
-<h3>Cell Level Modeling</h3>
+<p>
+Another parameter we investigated is the influence of cell density on the conductivity of the system. From percolation theory we have learned that this has a strong influence on the moment of percolation. The results are shown in figure 10.
+</p>
+<figure>
+<img src="https://static.igem.org/mediawiki/2014/f/ff/TUDelft2014_Colony_CellDens.png" width="60%" height="60%">
+<figcaption>
+Figure 10: The conductance of our chip ( \( 500\mu m \ by \ 500\mu m \) as function of time. Two different cell densities have been simulated. The red lines show simulations with a cell density of \( 2 \cdot 10^4 cells / mm^2\). The blue lines show the results with a cell density of \( 2.4 \cdot 10^4 cells / mm^2\)
+</figcaption>
+</figure>
+<br>
+<p>
+The most important question of all is: Will our system work, and under what conditions? To answer this question we have done the exact same simulation for the gene level and the cell level as elsewhere, except with the only difference that we halved the promoter strength of CsgB. This simulates less CsgB protein production, which corresponds to less TNT/DNT present in the system. The results are shown in figure 11.
+</p>
+<figure>
+<img src="https://static.igem.org/mediawiki/2014/6/6d/TUDelft2014_Colony_ActivationStrength.png" width="60%" height="60%">
+<figcaption>
+Figure 11: The conductance of our chip ( \( 500\mu m \ by \ 500\mu m \) as function of time. Two different activation strengths have been simulated. The red lines (same as in figure 6 and 7) show a 100% activated CsgB promoter. The blue lines show the results for a promoter that is activated for 50%.</figcaption>
+</figure>
 <br>
-<a name="Colony Level"></a>
-<h3>Colony Level Modeling</h3>
+<p>
+Based on the above three figures, there are a couple of recommendations we have for the design of the chip:
+<ul>
+<li>
+Based on figure 9, we recommend to decrease the distance between the electrodes. This will increase the conductance without increasing the uncertainty. From the shape of the curve, we expect that we do not have to wait very long (only a couple of hours) to do a measurement. After a steep increase in conductance, the increase in conductance will be relatively small in time. Furthermore, the uncertainty due to the placement of our cells increases in time.
+</li>
+<li>
+Based on figure 10, we recommend high cell density for increased signal strength of our system. However, it is crucial that the cell density is kept constant throughout the measurement. Cell density influences the conductance significantly and when the cell density is not kept constant, it will influence the conductance measurements in such a way that an increase or decrease in conductance can not only be attributed to an increase or decrease in DNT any longer.
+</li>
+<li>
+From figure 11, we learned that we have to be very careful with our measurement time in order to be able to draw conclusions about the concentration of TNT/DNT. The influence of the induction strength of the CsgB promoter is very small and its influence decreases after the initial amount of CsgA protein have formed curli. A recommendation for the wet lab is to increase the waiting time prior to induction of the promoter and thereby increasing the initial amount of CsgA protein. This will increase the relative differences in conductance between different induction strengths of the CsgB promoter. In our model, this increase in relative differences in conductance happens between 0 and 2.5 hours. <br>
+The time this increase in relative differences in conductance happens, is strongly dependent on the parameter k in our <a href="https://2014.igem.org/Team:TU_Delft-Leiden/Modeling/Curli/Gene#simplifiedgenelevel">simplified gene model</a>. For our model, we could fit this parameter to conductance curves obtained in the wet lab. Measuring conductance until observing a linear increase in time will give the most useful information for obtaining the inducer concentration, as the time of this transition depends on the strength of induction.
+</li>
+</ul>
+</p>
 </div>