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 Module</h2>
 <p>
-<font color="red">introduction curli model, explain three different modules</font>
+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>
-<br>
-<p>
-<font color="red">summary of the conclusions</font>
-</p>
-</body>
+<div class="tableofcontents">
-</html>
-{{:Team:TU_Delft-Leiden/Templates/tocmodnonfloat}}
+<ul> <p> Curli Module </p>
-<html>
+               <ul>
-<body>
+                    <li>
-<a name="Gene Level"></a>
+                    <a href="https://2014.igem.org/Team:TU_Delft-Leiden/Modeling/Curli/Gene">
-<h3> Gene Level Modeling </h3>
+                    <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>
-<p>
+                               <li>
-We will start with the modeling of the expression of curli on the gene level. Proteins that are dedicated to the curli formation are CsgA/B/D/E/F/G [1]. CsgA is the main building block of the curli. When produced, this protein 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 and connect the cell membrane to the first CsgA protein in the curli fibril. Once CsgB is located on the outside of the cell surface, the CsgA can polymerize onto the starting curli fibril. <br>
+                               <a href="https://2014.igem.org/Team:TU_Delft-Leiden/Modeling/Curli/Gene#simplifiedgenelevel">
-In the constructs we made in the wet lab, CsgA is continuously being produced. However, in our constructs the CsgB gene is placed under the control of a landmine promoter, activated by either TNT or DNT <font color="red">reference to landmine</font>. 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.
+                               <p>Simplified Gene Level Modeling</p>
-</p>
+                               </a>
+                               </li>
+                          </ul>
+                    </li>
-<br>
+                    <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>
-<a name="extendedgenelevel"></a>
+                               <li>
-<h4> Extensive Gene Level Modeling</h4>
+                               <a href="https://2014.igem.org/Team:TU_Delft-Leiden/Modeling/Curli/Cell#FittingCurliDensity">
+                               <p>Fitting the Curli Density</p>
+                               </a>
+                               </li>
-<p>
+                               <li>
-<font color="red">to be written, low priority</font>
+                               <a href="https://2014.igem.org/Team:TU_Delft-Leiden/Modeling/Curli/Cell#conductiveradius">
-</p>
+                               <p>Conductive Radius of the Cell</p>
+                               </a>
+                               </li>
+                          </ul>
+                    </li>
-<br>
+                    <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>
-<a name="simplifiedgenelevel"></a>
+                               <li>
-<h4> Simplified Gene Level Modeling</h4>
+                               <a href="https://2014.igem.org/Team:TU_Delft-Leiden/Modeling/Curli/Colony#resistivity">
+                               <p>Resistance</p>
-<p>
+                               </a>
-Though the model described above, providing that all rates are known, has a more accurate (though still simplified) representation of the curli assembly system, we have chosen to decrease the complexity further to the bare essentials, as most of the production rates cannot be found in literature. Measuring the accurate rates in the wet lab is, within the scope of this project, infeasible and therefore, we constructed a model that only includes the rate limiting step of the system as this will mostly determine the dynamics of the system.<br>
+                               </li>
-First of all, we investigated if the diffusion of the CsgA and CsgB proteins to their final destination is the rate limiting step in curli formation. From the literature and the wet lab, we know that the system response to the induction by TNT or DNT is in the order of hours font color="red">[reference]</font>. If diffusion is the rate limiting step, it would mean that CsgA and CsgB proteins would pile up inside and outside the cell, because it takes a long time for them to travel to their final destination, the end of a growing curli fibril and the outer membrane, respectively. A quick calculation shows that after one second, the displacement of a spherical particle with radius \(r = \ 10 \ nm\) is 6.6 μm due to Brownian motion in liquid water at room temperature using equation 1; many times the bacterial radius! Hence, we conclude that diffusion is not rate limiting [4].
+                               <li>
-</p>
+                               <a href="https://2014.igem.org/Team:TU_Delft-Leiden/Modeling/Curli/Colony#recommendations">
+                               <p>Recommendations for product design and wet lab</p>
-$$ \bar{x}^2  = \ \frac{k_b T t}{3 \pi \eta r} \tag{1}$$
+                               </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>
 <p>
-What we do expect to be the rate limiting step for curli formation is the large amount of CsgA and CsgB proteins that have to be produced. Hence, we expect the production rate of one of these proteins to be the rate limiting step. Instead of including the intermediate steps, we have implemented the production of the CsgA and CsgB proteins with one reaction and associated production rate each. These rates have to be measured in the lab. We will use the following system of equations:
+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>
-$$ \emptyset \xrightarrow{p_{A}} \ CsgA_{free} \tag{2} $$
-$$ \emptyset \xrightarrow{p_{B}} \ CsgB \tag{3} $$
-$$ CsgA_{free} + \ CsgB \xrightarrow{k} \ CsgA_{curli} + \ CsgB \tag{4} $$
 <p>
-Reactions 2 and 3 represent the production of CsgA and CsgB proteins, respectively. Equation 4 represents the growing of a curli fibril, where a curli fibril reacts with a free CsgA protein to become part of the curli. In reality, this reaction only happens at the end of the curli fibrils. In our model, we assume a homogeneous concentration of all the substances and we cannot discriminate between curli subunits. It is theoretically possible to model the system as an infinite amount of possible reactions that can take place to increase a curli fibril 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 at <a href="url">the cell level</a>. Therefore, we decided to model the growing of curli at the gene level as reaction 4. We assume that each CsgB protein is the start of a curli fibrils, thus the concentration of CsgB equals the concentration of curli. We can do this, because we showed that the diffusion of CsgA and CsgB proteins to their final destination is not the rate limiting step. Therefore, nearly all the CsgB proteins will be the beginning of a curli fibril in reality and our assumption is valid. <br>
+A figure of our three-layered model is displayed below.
-So, in reaction 4 we let a free CsgA protein react with a curli fibril to a CsgA protein that is part of that curli and the curli itself again, as it is immediatily again availible for the next reaction with a free CsgA protein to grow even more. Therefore, curli growth is dependent on the rate k and the concentration of \(CsgA_{free}\) and CsgB.
 </p>
-<br>
-<p>
-Writing reactions 2-4 into differential equations results in:
-</p>
-$$ \frac{d}{dt} [CsgA_{free}] = \ p_{A} - \ k [CsgA_{free}][CsgB] \tag{5.1} $$
-$$ \frac{d}{dt} [CsgB] = \ p_{B} \tag{5.2} $$
-$$ \frac{d}{dt} [CsgA_{curli}] = \ k [CsgA_{free}][CsgB] \tag{5.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\). At this time there are no curli present, so \([CsgB]|_{t=0} = \ [CsgA_{curli}]|_{t=0}= \ 0\). However, the CsgA promoter is continuously active, so we expect to have an initial concentration \(A_0\) of free CsgA proteins at time  \(t= \ 0\).
-</p>
-<br>
-<p>
-The solution to equation 5.2 is trivial:
-</p>
-$$ [CsgB] = \ p_B t \tag{6}$$
-<p>
-Substituting this into equation 5.1 results in:
-</p>
-$$ \frac{d}{dt} [CsgA_{free}] = \ p_{A} - \ K p_B [CsgA]t \tag{7} $$
-<p>
-It can easily be proven that a first order differential equation of the form
-</p>
-$$ y(t)' + \ f(t)y(t) = \ g(t) $$
-<p>
-has a solution of the form
-</p>
-$$ y(t) = \ e^{-F(t)} \int{g(t) e^{F(t)}  dt} + \ y_0 e^{-F(t)} $$
-<p>
-where \(F(t)= \int{f(t) dt}\). In our case, \(f(t) = \ k p_B t\) and \(g(t) = \ p_A\). This yields equation 8.
-</p>
-$$ [CsgA_{free}] = \ p_A e^{\frac{-k \ p_B t^2}{2}} \int{e^{\frac{k \ p_B t^2}{2}} dt} + \ C_{1} e^{\frac{-k \ p_B t^2}{2}} = \ p_A e^{\frac{-k \ p_B t^2}{2}} \int_{0}^{t}{e^{\frac{k \ p_B \tau^2}{2}} d\tau} + \ C_{2} e^{\frac{-k \ p_B t^2}{2}} \tag{8} $$
-<p>
-One with a keen eye may recognize the Dawson function (equation 9):
-</p>
-$$ D_+ (x) = \ e^{-x^2 } \int_{0}^x{e^{y^2} dy} \tag{9}  $$
-<p>
-As in our case, \(x^2 = \ k p_B t^2 \) and \(y^2 = k p_B \tau^2 \) and equation 10 obtained.
-</p>
-$$ [CsgA_{free}] = \ \frac{p_A D_+ (t\sqrt{\frac{k \ p_B}{2}})}{\sqrt{\frac{k \ p_B}{2}}}  + \ C_{2} e^{\frac{-k \ p_B t^2}{2}}  \tag{10}$$
-<p>
-Using the boundary condition \([CsgA_{free}]|_{t=0}= \ A_0\), the expression for the concentration of free CsgA proteins becomes:
-</p>
-$$ [CsgA_{free}] = \ \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, we can fill in equations 11 and 6 into equation 5.3, which gives us equation 12.
-</p>
-$$ \frac{d}{dt} [CsgA_{curli}] = \ 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{12} $$
-<p>
-For the parameters \(p_{A}\), \(p_{B}\), \(k\) and \(A_0\), we have estimated the following values <font color="red">explain</font>:
-</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 12 with the parameter values in table 1 yields the graph shown in figure 1. <font color="red">insert caption</font>
-</p>
-<br>
+<h3>Click in the figure to move to the corresponding page.</h3>
 <figure>
-<img src="https://static.igem.org/mediawiki/2014/7/74/Delft2014_ExactSol31.png" width="100%" height="100%">
+<img src="https://static.igem.org/mediawiki/2014/f/fa/TU_Delft_2014_ModellingFigureCurli_small.png" width="800" height="920" usemap="#ModelMap">
-<figcaption>
-Figure 1: hoi!
-</figcaption>
-</figure>
-<br>
+<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">
-<p>
+</map>
-Figure 1 shows a steady production of CsgB. \(CsgA_{curli}\) concentration at \(t= \ 0\) is zero as expected, since there is no CsgB at that point. In the next few hours, \(CsgA_{curli}\) concentration peaks. We think that this is due to the high concentration of \(CsgA_{free}\) that is present at \(t= \ 0\). In figure 2, curli growth as function of time is plotted for different initial concentrations of \(CsgA_{free}\). <font color="red">insert caption</font>
-</p>
-<figure>
-<img src="https://static.igem.org/mediawiki/2014/e/e9/Delft2014_DifferentA0.png" width="100%" height="100%">
 <figcaption>
-Figure 2: por
+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>
@@ Line 210: / Line 146: @@
 <br>
-<p>
-We conclude the following from figure 2: <br>
-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}\), \(A_0\). <br>
-Secondly, increasing the initial concentration of \(CsgA_{free}\), \(A_0\), 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. <br>
-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.
-</p>
-<br>
 <p>
-<font color="red">write some conclusion about rate limiting step and which information we gained from this level and what will be used in the subsequent levels</font>
+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>
-</p>
+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>
-<br>
-<a name="Cell Level"></a>
-<h3>Cell Level Modeling</h3>
-<p>
-Now that the growth rate of curli and production of CsgB protein as function of time is obtained, the conductivity as a function of time can be computed. The relevant length scale is the cell length, or the micrometre scale. The approach we used for this is relatively simple: <br>
 <ul>
 <li>
-We discretize the amount of curli subunits (\(CsgA_{curli}\) in the gene level model) and CsgB proteins that have to be added for each time step.
+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>
-At each time step, we add more curli subunits to growing curli fibrils. Also, we add more new curli fibrils to the model.
+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>
-From the density of the curli fibrils around the cell as a function of the radius, we calculate the conductive radius of the cell.
+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>
-</p>
-<br>
+<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>
-<a name="discretization"></a>
-<h4>Discretization of Gene Level Model</h4>
-<p>
-We have discretized equations 5.2 and 12 in N time steps. These give the expected number of new CsgB proteins and curli subunits for each time step, as we plotted the solution of these two equations in figures 1 and 2. From these figures we determine the expected number of new CsgB proteins and curli subunits for each time step. However, a fundamental assumption in deterministic modeling is that the concentration is continuous. In reality, the amount of added curli subunits is discrete, since we cannot add half a curli subunit. <br>
-Furthermore, in the gene level model we did not take into account the statistical variation of gene transcription and adding of curli subunits; sometimes less and some times more curli subunits are added with respect to the expected value. To include this in the cell level model, we drew the amount of new curli subunits from a Poisson distribution where λ equals the expected amount of added subunits. <br>
-So, for each time step we now have \(B_n\) new CsgB proteins and \(C_n\) new curli subunits, where \(C_n\) varies for each time step, as it is drawn from a Poisson distribution. A assumption of this distribution is that the time at which a new curli subunit is added, is uncorrelated to the time at which the previous curli subunit was added, we think this is a fair assumption. Note that the cell level model we made, accounts for the stochasticity of adding curli subunits, but not for the stochasticity of gene expression, so for the production of CsgB protein. The value \(B_n\) and the Poisson distribution are determined from figures 1 and 2.
-<font color="red">say something about the time steps, how much time represents each step and determine Bn and Cn from figures</font>
-</p>
 <br>
-<a name="building"></a>
-<h4>Building the Curli Fibrils</h4>
 <p>
-Firstly, \(B_n\)  CsgB proteins are added to our model that mark the starting points for new curli fibrils. These new curli fibrils are located at random points on a sphere with radius r, which represents the cell. The radius r is chosen such that the volume of the cell is\(\ \sim 1.1 \ \mu m^3\) [5]. A CsgB protein is modeled by a line of length 4 nm that points radially outward, perpendicular to the cell surface <font color="red">[source]</font>. In reality, the distribution of CsgB on the cell surface is not uniformly distributed [6]. However, we assumed uniformly distributed CsgB to keep our model prehensile. This is a point that may be used to further improve the model.
+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
-</p>
+$$ \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>
-<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>
-Next, \(C_n\), which is drawn from the Poisson distribution, where λ equals the expected amount of added curli subunits, new curli subunits are added to curli fibrils by repeating the following process \(C_n\) times:
-<ul>
-<li>
-Firstly, a random curli fibril is selected, e.g. curli number k. A curli fibril is represented by a 3 (the x, y and z coordinates) by l+1 matrix, where l is the amount of curli subunits of the curli fibril and the origin is chosen to be the center of the sphere. Thus, by storing the ending coordinates of each curli subunit, we know the starting and end coordinates of each curli subunit. The curli subunits are modeled by a line of length 4 nm <font color="red">[source]</font>.
-</li>
-<li>
-Secondly, the polar angle in spherical coordinates of the last curli subunit is computed, \(\theta_{1}\).
-</li>
-<li>
-Thirdly, the new curli subunit has a small angular deviation with respect to the previous one. This polar angle \(\theta_{2}\) is chosen from a Gaussian distribution with parameters N(0,σ). σ is chosen such that the persistence length, the distance over which a fibril has bend by \(90^{\circ}\) and has ‘lost’ its directional information, is 4 µm. The azimuthal angle ϕ is completely random between 0 and 2π radians, and chosen from an uniform distribution.
-</li>
-<li>
-Fourthly, for the new curli subunit for which we determined \(\theta_{2}\) and ϕ, the polar angle is determined to be \(\theta_{1} + \theta_{2}\). We now know the length of the new curli subunit (4 nm), its polar angle and its azimuthal angle. Subsequently, we add it to the previous curli subunit of the fibril and calculate the ending coordinate of the added curli subunit from its length, polar angle and azimuthal angle and the ending coordinate of the previous curli subunit. This calculated ending coordinate of the added curli subunit is stored in the matrix that represents the curli fibril.
-</li>
-</ul>
-</p>
-<br>
-<p>
-The angular deviation σ is a critical parameter in our model. Increasing this value increases the flexibility of our curli, where decreasing this value increases the stiffness of the curli. This is shown in figure 3. If the length of one subunit is 4 nm and the total persistence length is 4 µm, then \(\sigma = \ 3.47^{\circ}\). Furthermore, we think that it is justified to add the curli subunits one at a time to a random curli. We expect no discrimination of the CsgA proteins for binding to a large or small curli or one that has recently gotten a new curli subunit. <font color="red">figure caption</font>
 </p>
 <figure>
-<img src="https://static.igem.org/mediawiki/2014/a/ab/TUDelft_2014_AnglePersistence.png" width="100%" height="100%">
+<img src="https://static.igem.org/mediawiki/2014/3/3f/TUDelft_2014_FittingCurli.png" width="60%" height="60%">
 <figcaption>
-Figure 3: por
+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>
-An illustrative view of what our cell looks like during the adding of curli subunits is shown in figure 4. This figure is created when just a few curli were added (\( \sim 1/2 \ hour\)). A similar figure after \(t = \ 10 \ hr\) would look like a fuzzy ball of curli. <font color="red">figure caption</font>
-</p>
 <figure>
-<img src="https://static.igem.org/mediawiki/2014/5/5b/TUDelft_2014_fuzzycurlicell.png" width="100%" height="100%">
+<img src="https://static.igem.org/mediawiki/2014/e/e1/TUDelft_2014_cellradius.png" width="60%" height="60%">
 <figcaption>
-Figure 4: por
+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>
@@ Line 312: / Line 200: @@
 <p>
-Now that we have a model of a cell with growing curli, we want to extract relevant data for the colony level modeling. Ideally, the resistance as function of radius and time would be calculated by looking at connections between the curli fibrils. However, this requires insight of the behavior of the curli on the nanoscopic scale. For instance, what is the conductivity of a single curli fibril with gold nanoparticles and what is the critical distance between the fibrils that make them connect? After an extensive literature study, we have decided to simplify this model. Furthermore, when interactions between the curli fibrils have to be taken into account, the model gets too computationally expensive.
+An illustrative view of what our cell looks like during the adding of curli subunits is shown in figure 5.
 </p>
-<br>
-<p>
-<font color="red">[write something about the part where we tried the percolation on this level], low priority</font>
-</p>
-<br>
+<figure>
-<p>
+<img class="theimage" src="https://static.igem.org/mediawiki/2014/e/e3/TU_Delft_2014_Curli_picturestill.png" alt="meh" width="60%" height="60%"/>
-To have a reasonable computational time, we decided to extract our parameters for the colony level modeling from the curli density around the cell. Figure 5 contains a histogram with the amount of curli subunits as a function of the cell radius after 10 hours. <font color="red">[add figure]</font> Note how no curli are found below the actual cell radius. It can be seen from the figure that there is a large peak, followed by a plateau. When this histogram is observed in time, you would notice that at first large curli are being created. Figure 6 shows the length of all curli after 10 hours. <font color="red">figure caption</font>
+<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%"/>
-</p>
-<figure>
-<img src="https://static.igem.org/mediawiki/2014/f/f1/TUDelft_2014_Curli_Length.png" width="100%" height="100%">
 <figcaption>
-Figure 6: por
+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>
-Looking back at figures 1 and 2, the fact that, as can be seen in figure 6, the first curli are much longer that the later ones, can be explained by the fact that there is relatively large curli growth in the beginning, because few CsgB have been produced and therefore, only a few curli fibrils are available for CsgA proteins. After a couple of hours there are more CsgB proteins, thus more curli fibrils, but CsgA protein production does not increase. Therefore, the ratio [CsgA]/[curli fibrils] is much smaller than in the beginning and each curli will grow much slower. A consequence of this is that the ‘newer’ curli fibrils are much shorter. The first histogram of figure 7 shows the amount of curli as function of radius after a growth time of 10 hour. The second histogram shows the curli density as function of radius after a growth time of 10 hour. <font color="red">[add figure]</font>
-</p>
-<a name="conductiveradius"></a>
+<p>
-<h4>Conductive Radius of the Cell</h4>
+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>
-<p>
-We think that a reasonable approximation of the conductivity is the density of the curli around the cell as a function of the radius (see figure 8 in the middle, the blue line). When the density is higher, there are more gold particles, thus higher conductivity. In our simplest model we say that there is a critical density \(\rho_{crit}\) of curli that is needed to have conductivity. The density \(\rho_{curli}\) decreases as function of the radius. The largest radius where \(\rho_{curli} > \rho_{crit}\), we call the conductive radius \(r_{cond}\). With only this simple approximation we can calculate some interesting properties of our system: the time at which we expect percolation to happen and the resistivity of our system. Though this approximation seems to be rather arbitrary, we do have some reasoning for this: <font color="red">[add caption]</font>
-<ul>
 <li>
-First of all, the goal of this parameter is to get information about our system that will be calculated in colony level modeling. We use this parameter in colony level modeling to find connections between cells. To have a continuous path from one electrode to the other electrode, we must have a lot of cells that are connected to each other. In order to know when cells are connected to each other, we have to assume that everything at a certain radius from the cell is conductive; for this radius we use the critical density \(\rho_{crit}\). However, for this to be true the fibrils on one side of the cell must be connected to the fibrils on the other side. The <a href="https://2014.igem.org/Team:TU_Delft-Leiden/Modeling/Techniques#PercolationTheory">Percolation Theory</a> prescribes that this is a sharp transition as a function of the density, so we can choose \(\rho_{crit}\) in such a way that we are very sure that everything at \(\rho_{crit}\) from the cell is conductive.
+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>
-While the precise value of \(\rho_{crit}\) may be unknown and should be measured, we think that we can still get plenty of information about the qualitative behaviour of our system in advance. Figure 8 at the bottom shows the conductive radius \(r_{cond}\) as function of time using \(\rho_{crit}\) as shown in figure 8 in the middle as the red line. Increasing or decreasing \(\rho_{crit}\) would result in a similar \(r_{cond}\) as function of time. Hence, the qualitative behaviour is preserved.
+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:
-</li>
+$$  \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>
-Due to the simplifications that we made in order to be able to model our system, we cannot include interactions or cluster forming between the curli themselves. Using \(\rho_{crit}\), we have an elegant way to filter out modeling errors.
 </li>
 </ul>
 </p>
+<br>
+<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 src="https://static.igem.org/mediawiki/2014/9/95/TUDelft_Curli_prod_r_cond_rho_crit.png" width="100%" height="100%">
+<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 8: por
+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>
-Note that figure 8 at the bottom shows a discrete increase of the conductive radius. This is because \(\rho_{curli}\) is dependent on the radius, and because we run a computer simulation, the radius is therefore discretized. Moreover, if we rerun the script to calculate the conductive radius, we expect some variation in the curve, as \(\rho_{curli}\) is subject to stochasticity. Therefore, we have repeated our simulations on the cell level many times in order to get statistically valid results for the mean and standard deviation of \(r_{cond}\). <br>
+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.
-When we ran our simulation 100 times, we got the results displayed in figure 9. Figure 9 displays the curli density at \(\ t= \ 2 \ hours\) for all cells in the left figure. The orange line represents the average of the simulations. It can be concluded that the intercellular variation is relatively small. This makes sense, since the relative deviation of stochastic processes decreases with the sample size. In the right figure, the mean and standard deviation of the curli density as a function of the radius is shown. <font color="red">insert caption</font>
 </p>
-<figure>
+<img src="https://static.igem.org/mediawiki/2014/d/dc/TUDelft_2014_500x500ChipResults.png" width="60%" height="60%">
-<img src="https://static.igem.org/mediawiki/2014/d/d5/TUDelft_2014_2h_curli_density.png" width="100%" height="100%">
 <figcaption>
-Figure 9: por
+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>
@@ Line 382: / Line 280: @@
 <p>
-Another interesting figure to look at, is the density for different times, shown in figure 10. This figure shows that, corresponding with what we have seen previously, \(\rho_{curli}\) decreases as a function of the radius. Also, it decreases faster as a function of the radius in the first two hours. After two hours, we can see that the curli density increases only for small r, as mainly short curli are added to the system. This agrees with our previous results. <font color="red">insert caption</font>
+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="https://static.igem.org/mediawiki/2014/c/c4/TU_Delft_2014_curli_density.png" width="100%" height="100%">
+<img src="https://static.igem.org/mediawiki/2014/0/02/TUDelft2014_Percolation_InductionDifferences.png" width="60%" height="60%">
 <figcaption>
-Figure 10: por
+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>
-We use the data from figure 10 to make a first approximation for \(\rho_{crit}\). However, the exact value of \(\rho_{crit}\) has to be measured. We set \(\rho_{crit}\) equal to 1% of the maximum \(\rho_{curli}\) observed in figure 9. So, we set \(\rho_{crit} = \ 10^3 \ \frac{number \ of \ curli}{\mu m^{2}}\). The following figure is obtained. <font color="red">insert caption and check if rho_crit value is correct</font>
+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="https://static.igem.org/mediawiki/2014/e/e1/TUDelft_2014_cellradius.png" width="100%" height="100%">
+<img src="https://static.igem.org/mediawiki/2014/c/c5/TUDelft2014_Colony_DifferentSizes.png" width="60%" height="60%">
 <figcaption>
-Figure 11: por
+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>
@@ Line 409: / Line 309: @@
 <p>
-Again, the orange line represents the mean conductive radius. A sharp increase in the conductive radius can be observed for \(t < 1 \ hour\), and after \(t = \ 1 \ hour\) the conductive radius increases slowly. The cellular variation in the second regime is relatively large, as all the lines except the orange line in figure 11 represent single simulations of our system. Different values of \(\rho_{crit}\) result in different characteristic curves for \(r_{cond}\) as shown in figure 12. <font color="red">insert caption</font>
+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/e/ec/TUDelft_2014_Different_Critical_values.png" width="100%" height="100%">
+<img src="https://static.igem.org/mediawiki/2014/f/ff/TUDelft2014_Colony_CellDens.png" width="60%" height="60%">
 <figcaption>
-Figure 12: por
+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>
@@ Line 422: / Line 322: @@
 <p>
-From figure 12, we conclude that low values of \(\rho_{crit}\) result in a sharp increase of \(r_{cond}\) followed by a steady, slow increase of \(r_{cond}\) in time. During the steady, slow increase of \(r_{cond}\) in time, the cellular variation is relatively large. For high values of \(\rho_{crit}\), there is a delayed sharp increase of \(r_{cond}\) and less cellular variation.
+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>
-<br>
-<p>
+<figure>
-<font color="red">make a fit for \(\rho_{curli}\)</font> <br>
+<img src="https://static.igem.org/mediawiki/2014/6/6d/TUDelft2014_Colony_ActivationStrength.png" width="60%" height="60%">
-<font color="red">from this fit determine \(\rho_{crit}\) (transition value between the two different fits)</font> <br>
+<figcaption>
-<font color="red">plot \(r_{cond}\) for this found value of \(\rho_{crit}\)</font> <br>
+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>
-</p>
+</figure>
-<a name="Colony Level"></a>
-<h3>Colony Level Modeling</h3>
-<a name="percolation"></a>
-<h4> Percolation </h4>
-<p>
-<font color="red">still has to be written</font>
-</p>
 <br>
-<a name="resistivity"></a>
-<h4> Resistivity </h4>
 <p>
-To get a quantitative measure for the conductivity between two cells, we will quantify the overlap of two conducting spheres, where we assumed that the conducting spheres represent cells surrounded by curli filaments. To do this, we will subdivide the overlapping region in infinitesimal volumes \(dV\). The infinitesimal conductivity of such an infinitesimal volume is given by:
+Based on the above three figures, there are a couple of recommendations we have for the design of the chip:
-</p>
+<ul>
+<li>
-$$ d \sigma (y) = \ \frac{\rho_1}{r_1} dV \frac{\rho_2}{r_2} dV \tag{}$$
+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>
-<p>
+<li>
-The factor \( 1/r \) is introduced to account for the conductivity of the wires itself. Further away from the cell, the wires need a longer distance to go to the cell. Since we want to know the strongness of the connection between the cells we have to include this factor. For a straight line this inversely proportional to the distance. For a single curli fibril, this relation does not hold. However, we assume that the curli density is high, thus there are many connections between the curli themselves. Then there is a pathway from the origing to \( r \) roughly proportional to the distance from the cell. To find the total conductivity, we integrate on both sides. To account for the fact that both volume elements \(dV\) are the same, we make use of the Dirac-delta function \(\delta_3\) <font color="red">[source]</font>. This gives us the following:
+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.
-</p>
+</li>
+<li>
-$$ \sigma (y) = \int{ \frac{\rho_1(\vec{r_1})\rho_2(\vec{r_2})}{\vec{r_1}\vec{r_2}}\delta_3(\vec{r_2}-f(\vec{r_1}))d^3\vec{r_1}d^3\vec{r_2}} \tag{} $$
+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.
-<p>
+</li>
-The Dirac delta allows us to remove the \(\vec{r_2}\) dependence by expressing these in \(\vec{r_1}\). The still undetermined relation between \(\vec{r_1}\) and \(\vec{r_2}\) is given by \(\vec{r_2} = f(\vec{r_1})\). Applying this removes one of the two volume integrations. Using spherical coordinates, the resulting single volume integration can be written as:
+</ul>
-</p>
-$$ \sigma (y) = \int_{r_0}^{r_{max}} \int_0^{\theta_{max}(r)} \int_0^{2\pi} \rho(r_1)\rho_2(f(r_1))\frac{r_1}{r_2} \sin(\theta_1) d\phi_1 d\theta_1 dr_1 \tag{} $$
-<p>
-Here we have made use of the fact that the density \(\rho\) is only dependent on \(r\) and not on \(\phi\) and \(\theta
-\). The integral over \(phi\) is trivial and gives us a multiplication factor of \(2 \pi\):
-</p>
-$$ \sigma (y) = \ 2 \pi \int_{r_0}^{r_{max}} \int_0^{\theta_{max}(r)} \rho(r_1)\rho_2(f(r_1))\frac{r_1}{r_2} \sin(\theta_1) d\theta_1 dr_1 \tag{} $$
-<p>
-Now that we have reduced our integration to two dimensions, we will work out \(f(\vec{r_1})\). To do this, we introduce the vector from the origin of cell 1 to the origin of cell 2, \(\vec{y}\). This allows us to express \(\vec{r_2}\) in terms of \(\vec{y}\) and \(\vec{r_1}\):
-</p>
-$$ \vec{r_2} = \ \vec{y} - \ \vec{r_1} = \begin{bmatrix}y \\0\\ \end{bmatrix} - \begin{bmatrix} \frac{r_1}{r_2} \cos(\theta_1) \\r_1 \sin(\theta_1)\\ \end{bmatrix} \tag{} $$
-<p>
-Now it is straightforward to express \(r_2\) in terms of \(y\), \(r_1\) and \(\theta_1\):
-</p>
-$$ r_2 = \ |\vec{r_2}| = \ \sqrt{(y - r_1 \cos(\theta_1))^2 + \ r_1^2 \sin^2(\theta_1)} \tag{} $$
-<p>
-Plugging this in yields the following integral:
-</p>
-$$ \sigma (y) =  \ 2 \pi \int_{r_0}^{r_{max}} \int_0^{\theta_{max}(r)}  \frac{\rho(r_1)\rho_2 \left( \sqrt{(y - r_1 \cos(\theta_1))^2 + r_1^2 \sin^2(\theta_1)}\right) r_1 \sin(\theta_1)}{ \sqrt{(y - r_1 \cos(\theta_1))^2 + r_1^2 \sin^2(\theta_1)}}  d\theta_1 dr_1 \tag{} $$
-<p>
-We will now have a closer look at the boundary values for \(r\) and \(\theta\). For \(r\), the maximum value \(r_{max}\) is defined as the cut-off distance, i.e. the radial coordinate where the behavior of the density function changes. The lower boundary \(r_0\) is defined as the \(r_{max}\) from the other cell, so it is equal to \(r_0  = \ y - \ r_{max}\). <font color="red">'[figure to explain this]</font>
-</p>
-<br>
-<p>
-To find the upper limit for \(\theta\), \(\theta_{max}\), we will use the law of cosines. <font color="red">[figure to illustrate this]</font>:
-</p>
-$$ \ r_{max}^2 = \ r_1^2 + y^2  - \ 2 r_1 r_{max} \cos(\theta_{max}) \tag{} $$
-$$ \theta_{max} = \ \arccos \left(\frac{-\ r_{max}^2 + \ r_1^2 + y^2}{2 r_1 r_{max}}\right) \tag{} $$
-<p>
-We will now use the previously <font color="red">[link]</font> found fact that the curli density can be described as:
-</p>
-$$ \rho(r) = \ C_{1}e^{-\frac{r}{C_{2}}}  \tag{} $$
-<p>
-Plugging the boundary values and our expression for \(\rho(r)\), we find the following expression for the conductivity between two cells:
-</p>
-$$  \sigma (y) =  \ 2 \pi C_{3} C_{1} \int_{y - r_{max}}^{r_{max}} \int_0^{\arccos \left(\frac{-\ r_{max}^2 + \ r_1^2 + y^2}{2 r_1 r_{max}}\right)}  e^{-\frac{r_1}{C_{2}}} e^{-\frac{ \sqrt{(y - \ r_1 \cos(\theta_1))^2 + \ r_1^2 \sin^2(\theta_1)}}{C_{2}}} r_1^2 \sin(\theta_1)  d\theta_1 dr_1 \tag{} $$
-<p>
-This integral looks very complicated, but can algebraically be simplified with some clever substitutions. We can rewrite this integral by moving all terms independent of \( \theta \). Furthermore, using that \( sin^2 (\theta_1) + cos^2\theta_1 = 1 \) we get.
-</p>
-$$  \sigma (y) =  \ 2 \pi C_{3} C_{1} \int_{y - r_{max}}^{r_{max}}
-r_1^2 e^{-\frac{r_1}{C_{2}}}
- \int_0^{\arccos \left(\frac{-\ r_{max}^2 + \ r_1^2 + y^2}{2 r_1 r_{max}}\right)}
-  e^{-\frac{ \sqrt{y^2+r_1^2-2yr_1 cos( \theta_1 ) }}{C_{2}}}  \sin(\theta_1)  d\theta_1 dr_1 \tag{} $$
-<p>
-First we must recognize that we can substitute \( x= cos(\theta_1) \) such that \( dx = -\sin(\theta_1) d\theta_1 \). This results in:
-</p>
-$$  \sigma (y) =  - \ 2 \pi C_{3} C_{1} \int_{y - r_{max}}^{r_{max}}
-r_1^2 e^{-\frac{r_1}{C_{2}}}
- \int_0^{\arccos \left(\frac{-\ r_{max}^2 + \ r_1^2 + y^2}{2 r_1 r_{max}}\right)}
-  e^{-\frac{ \sqrt{y^2+r_1^2-2yr_1 x }}{C_{2}}}  dx dr_1 \tag{} $$
-<p>
-In the second integral we recognize something of the form \( \int e^{-\sqrt{a+bx}} dx \) with \( a= \frac{y^2+r_1^2}{C^2_2} \) and \(b=-\frac{2yr_1}{C^2_2} \). Substituting \( h= \sqrt{a+bx} \) with \( dx= \frac{2h}{b} dh \) yields.
-</p>
-$$ \frac{2}{b} \int he^{-h} dh= \frac{-2}{b} \int h d e^{-h} = \frac{-2}{b} (he^{-h} - \int e^{-h} dh )=\frac{-2 e^{-h}(h+1)}{b}$$
-<p>
-Now h, a and b x can be resubtituted:
-</p>
-$$ \frac{-2 e^{-h}(h+1)}{b}=\frac{- e^{-\sqrt{a+bx}}(\sqrt{a+bx}+1)}{b}=\frac{C_2 e^{\frac{\sqrt{y^2+r_1^2-2yr_1cos(\theta_1)}}{C_2}  }(\sqrt{y^2+r_1^2-2yr_1cos(\theta_1)}+C_2)}{yr_1}$$
-<p>
-For the upper boundary, the arccos will cancel with the cos resulting in.
-$$\frac{C_2 e^{-\frac{\sqrt{y^2+r_1^2-yr_1 \left(\frac{-\ r_{max}^2 + \ r_1^2 + y^2}{r_1 r_{max}}\right)}}{C_2}  }
-\left( \sqrt{y^2+r_1^2-yr_1 \left(\frac{-\ r_{max}^2 + \ r_1^2 + y^2}{r_1 r_{max}}\right)}+C_2 \right)}{yr_1}$$
-</p>
-<p>And for the lower boundary this results in:</p>
-$$ \frac{C_2 e^{-\frac{|y-r_1|}{C_2}  }(|y-r_1|+C_2)}{yr_1}$$
-Thus the entire integral now becomes
-$$  \sigma (y) =  \frac{ 2 \pi C_{1} C_2 C_{3}}{y} \int_{y - r_{max}}^{r_{max}}
-r_1 e^{-\frac{|y-r_1|+r_1}{C_2}  }(|y-r_1|+C_2)
--
-r_1 e^{-\frac{r_1+\sqrt{y^2+r_1^2-yr_1 \left(\frac{-\ r_{max}^2 + \ r_1^2 + y^2}{r_1 r_{max}}\right)}}{C_2}  }
-\left( \sqrt{y^2+r_1^2-yr_1 \left(\frac{-\ r_{max}^2 + \ r_1^2 + y^2}{r_1 r_{max}}\right)}+C_2 \right)
- dr_1 \tag{} $$
-<p>
-<font color="red">solving this integral numerically, storing values in table</font> <br>
-<font color="red">giving the conductivity as weights to the adjacency matrix</font> <br>
-<font color="red">from these weights (conductivity) determine if cells are connected or not</font> <br>
-<font color="red">from the resulting adjacency matrix, calculate the resistivity of the system using graph theory</font>
-</p>
-<br>
-<h3> References </h3>
-<p>
-<font color="red">still has to be made</font>
 </p>