(19 intermediate revisions by 4 users not shown) | |||
Line 1: | Line 1: | ||
{{Newcastle}} | {{Newcastle}} | ||
<html> | <html> | ||
+ | <script> | ||
+ | <script type="text/x-mathjax-config"> | ||
+ | MathJax.Hub.Config({ | ||
+ | TeX: { equationNumbers: { autoNumber: "all" } } | ||
+ | }); | ||
+ | </script> | ||
+ | </script> | ||
<style> | <style> | ||
+ | #togglefig | ||
+ | { | ||
+ | mix-blend-mode: multiply; | ||
+ | max-width: 500px; | ||
+ | } | ||
+ | #reflist li | ||
+ | { | ||
+ | font-family: 'Didact Gothic'; | ||
+ | font-size: 16px; | ||
+ | color: #292000; | ||
+ | } | ||
+ | #integrationfig | ||
+ | { | ||
+ | width: 350px; | ||
+ | } | ||
+ | |||
+ | hr { | ||
+ | border-top: 1px solid #292000; | ||
+ | } | ||
+ | |||
#box1{ | #box1{ | ||
float: left; | float: left; | ||
Line 29: | Line 56: | ||
font-size: 16px; | font-size: 16px; | ||
} | } | ||
− | + | h3{ | |
+ | font-size: 26px; | ||
+ | } | ||
p{ | p{ | ||
color: #000; | color: #000; | ||
Line 60: | Line 89: | ||
</style> | </style> | ||
<body> | <body> | ||
− | <div class="column full_size"> | + | <div class="column full_size" style=" |
− | < | + | margin-top: 50px;"> |
+ | <h3>Modelling</h3> | ||
+ | <p>This page documents the dry lab experiments relating to our modelling.</p> | ||
+ | <h2 id="ode-modelling">ODE Modelling</h2> | ||
− | <p>As part of our work, we performed many different types of modelling, from rule-based modelling of our genetic parts to physics analysis of our microfluidics components. On this page, we discuss the ordinary differential equation (ODE) modelling of our components for use in a web-based ‘simulator’. We hope that | + | <p>As part of our work, we performed many different types of modelling, from rule-based modelling of our genetic parts to physics analysis of our microfluidics components. On this page, we discuss the ordinary differential equation (ODE) modelling of our components for use in a web-based <a href="https://2016.igem.org/Team:Newcastle/Software/Hosted/Models">‘simulator’</a>. We hope that our ‘simulator’ will demonstrate how our genetic constructs could fit into an electronic circuit.</p> |
− | < | + | <h2 id="simulator-core">Simulator Core</h2> |
<p>The core of our simulator features a set of ODEs that describe an abstract version of the transcriptional (Equation 1) and translational (Equation 2) processes in a cell.</p> | <p>The core of our simulator features a set of ODEs that describe an abstract version of the transcriptional (Equation 1) and translational (Equation 2) processes in a cell.</p> | ||
Line 80: | Line 112: | ||
k_{mRNA\ degradation} = 0.0058 \\ | k_{mRNA\ degradation} = 0.0058 \\ | ||
k_{translation} = 0.1 | k_{translation} = 0.1 | ||
+ | \tag{3} | ||
</script> <br> | </script> <br> | ||
Since our code is freely available, a modeller with different simulation requirements could easily change these parameters if they were meaningful in their scenario. We also constrain the other parameters in our simulator so that the simulation is stable when the user is given control over these values. For example, it might be logical for the user to select between a weak, medium or strong promoter. In this scenario, the user should not be required to know the exact numerical values of their choice. Although more advanced users may wish to be able to input their own values. For this reason, our simulation uses sane defaults for the transcription rate, <script id="MathJax-Element-1892" type="math/tex">k_{transcription}</script>, which is analogous to promoter strength. <br> | Since our code is freely available, a modeller with different simulation requirements could easily change these parameters if they were meaningful in their scenario. We also constrain the other parameters in our simulator so that the simulation is stable when the user is given control over these values. For example, it might be logical for the user to select between a weak, medium or strong promoter. In this scenario, the user should not be required to know the exact numerical values of their choice. Although more advanced users may wish to be able to input their own values. For this reason, our simulation uses sane defaults for the transcription rate, <script id="MathJax-Element-1892" type="math/tex">k_{transcription}</script>, which is analogous to promoter strength. <br> | ||
Line 86: | Line 119: | ||
medium = 0.03 \\ | medium = 0.03 \\ | ||
weak = 0.001 | weak = 0.001 | ||
+ | \tag{4} | ||
</script> <br> | </script> <br> | ||
− | In our design work, users thought being constrained to these three options limited their ability to fully explore the biological system. For example, it would not allow the user to see the variation in gene product caused by slight noise in the promoter strength value, e.g. switching between | + | In our design work, users thought being constrained to these three options limited their ability to fully explore the biological system. For example, it would not allow the user to see the variation in gene product caused by slight noise in the promoter strength value, e.g. switching between 0.03 and 0.1 quickly. As our system is primarily a way of exploring the constructs we have built we wanted to limit the end user as much as possible. Therefore we set out plans to replace the notion of strong/medium/weak promoters with a more fine-grained scale. One that the user might control with a simple slider for example. To achieve this feat we fit our existing values using regression to a 10 point scale. This gave us an equation that we can use to map a user’s idea of ‘promoter strength’ to a transcription rate. As shown in Figure 1, the resulting equation was <script id="MathJax-Element-1894" type="math/tex">k_{transcription}=0.0008x^{2} + 0.0017x + 0.001</script>. Note that we deliberately ensured that <script id="MathJax-Element-1895" type="math/tex">k_{transcription} \ne 0 \ \forall x \in \{0..10\}</script> to reflect the fact that there is always <em>some</em> background level transcription. In layman's terms, no biological promoter can be completely ‘switched off’.</p> |
<p><figure><img src="https://static.igem.org/mediawiki/2016/e/e7/T--Newcastle--Promoter_Strength_Scale.png" alt="Graph showing trend line and equation fitted to our promoter scale."><figcaption>Figure 1: Trend line and equation fitted to our promoter scale.</figcaption></figure></p> | <p><figure><img src="https://static.igem.org/mediawiki/2016/e/e7/T--Newcastle--Promoter_Strength_Scale.png" alt="Graph showing trend line and equation fitted to our promoter scale."><figcaption>Figure 1: Trend line and equation fitted to our promoter scale.</figcaption></figure></p> | ||
Line 95: | Line 129: | ||
\frac{d(mRNA[x])}{dT} = k_{transcription}[x] - k_{mRNA\ degradation} \cdot mRNA[x] \\ | \frac{d(mRNA[x])}{dT} = k_{transcription}[x] - k_{mRNA\ degradation} \cdot mRNA[x] \\ | ||
\frac{d(Protein[x])}{dT} = (k_{translation} \cdot mRNA[x]) - (k_{protein\ degradation}[x] \cdot Protein[x]) | \frac{d(Protein[x])}{dT} = (k_{translation} \cdot mRNA[x]) - (k_{protein\ degradation}[x] \cdot Protein[x]) | ||
− | </script></p> | + | \tag{5}</script></p> |
Line 101: | Line 135: | ||
<h3 id="justification">Justification</h3> | <h3 id="justification">Justification</h3> | ||
− | <p>We chose to build our simulator around a set of ODEs rather than using a more complex system that produced models in a format like CellML or SBML because it met our design goals. Primarily, this was ease of development and integration into the iGEM wiki system. Our main aim was to allow the simulator to be embedded in the wiki. In turn we felt that it was necessary to support the aims of the iGEM wiki system in preventing ‘bit-rot’, that is not having some aspects of our project hosted somewhere else where their uptime could not be guaranteed. This meant we needed a purely client side implementation and couldn’t ‘call out’ to other computers to perform the computation. When researching how we would construct our system we determined that the most straightforward of modelling our system in JavaScript without using any libraries or complex dependencies was | + | <p>We chose to build our simulator around a set of ODEs rather than using a more complex system that produced models in a format like CellML or SBML because it met our design goals. Primarily, this was ease of development and integration into the iGEM wiki system. Our main aim was to allow the simulator to be embedded in the wiki. In turn we felt that it was necessary to support the aims of the iGEM wiki system in preventing ‘bit-rot’, that is not having some aspects of our project hosted somewhere else where their uptime could not be guaranteed. This meant we needed a purely client side implementation and couldn’t ‘call out’ to other computers to perform the computation. When researching how we would construct our system we determined that the most straightforward of modelling our system in JavaScript without using any libraries or complex dependencies was to hard-code the ODEs.</p> |
<p>Another aspect that makes the ODEs approach useful for us is that the level of abstraction they capture i.e., transcription and translation, corresponds directly with the level at which we design our BioBrick constructs. This makes it easy to translate between formats and allows for us to rapidly add new knowledge to the simulator as we develop our constructs. And vice versa, the simulations at this level translate directly to components in the system which we can control in our constructs. For example, we can directly influence the number of proteins in the cell by increasing the number of coding sequences on a construct. Whereas we can’t influence variables like the amount of free mRNA through our constructs alone.</p> | <p>Another aspect that makes the ODEs approach useful for us is that the level of abstraction they capture i.e., transcription and translation, corresponds directly with the level at which we design our BioBrick constructs. This makes it easy to translate between formats and allows for us to rapidly add new knowledge to the simulator as we develop our constructs. And vice versa, the simulations at this level translate directly to components in the system which we can control in our constructs. For example, we can directly influence the number of proteins in the cell by increasing the number of coding sequences on a construct. Whereas we can’t influence variables like the amount of free mRNA through our constructs alone.</p> | ||
− | <p>There are some trade-offs to this approach. For instance, we lose access to tools that perform local and global sensitivity analysis, parameter estimation, etc., which limits the ‘depth’ of our models. We also lose out on interoperability with future software through interchange standards like SMBL. However, we believe that this doesn’t subtract from the quality of our models because their intended purpose is to encourage the user to think about the implications of interfacing bacteria with electronics rather than providing detailed models. In other words our present systems are not complex enough to warrant modelling at a lower level of abstraction in the simulator.</p> | + | <p>There are some trade-offs to this approach. For instance, we lose access to tools that perform local and global sensitivity analysis, parameter estimation, etc., which limits the ‘depth’ of our models. We also lose out on interoperability with future software through interchange standards like SMBL. However, we believe that this doesn’t subtract from the quality of our models because their intended purpose is to encourage the user to think about the implications of interfacing bacteria with electronics rather than providing detailed models. In other words, our present systems are not complex enough to warrant modelling at a lower level of abstraction in the simulator.</p> |
<p>That is not to say that we didn’t also perform lower level modelling. We also used rule-based modelling in the BioNetGen language using the RuleBender graphical tool. This allowed us to capture our systems in more detail and to ‘sanity check’ the simplified models we use here in our simulator.</p> | <p>That is not to say that we didn’t also perform lower level modelling. We also used rule-based modelling in the BioNetGen language using the RuleBender graphical tool. This allowed us to capture our systems in more detail and to ‘sanity check’ the simplified models we use here in our simulator.</p> | ||
Line 111: | Line 145: | ||
<h3 id="integration">Integration</h3> | <h3 id="integration">Integration</h3> | ||
− | <p>There are a number of ways to integrate the differential equations presented above. The simplest is to iteratively compute the difference and add these to a ‘running total’ over the lifetime of the simulation. That is, for each time step compute <script id="MathJax-Element-1898" type="math/tex">y_{n+1} = y_{n} + hf(t_{n},y_{n})</script> where h is the step size, f is a differential equation such as those above. This method is known as <a href="https://en.wikipedia.org/wiki/Euler_method">Euler’s method</a>. It can be | + | <p>There are a number of ways to integrate the differential equations presented above. The simplest is to iteratively compute the difference and add these to a ‘running total’ over the lifetime of the simulation. That is, for each time step compute <script id="MathJax-Element-1898" type="math/tex">y_{n+1} = y_{n} + hf(t_{n},y_{n})</script> where h is the step size, f is a differential equation such as those above. This method is known as <a href="https://en.wikipedia.org/wiki/Euler_method">Euler’s method</a>. It can be visualised as fitting a series of straight lines to the given differential equation. It is, therefore, more accurate at small step sizes (the global error is proportional to the step size and the local error is proportional to the square of the step size). This is the method we use in our simulator.</p> |
− | <p><figure><img src="https://static.igem.org/mediawiki/2016/e/e0/T--Newcastle--Euler-Graphical.png" alt="A graphical representation of Euler's method">.<figcaption> | + | <p><figure><img id="integrationfig" src="https://static.igem.org/mediawiki/2016/e/e0/T--Newcastle--Euler-Graphical.png" alt="A graphical representation of Euler's method">.<figcaption>Figure 2: Illustration of Euler's method for integration. Source: <a href="https://en.wikipedia.org/wiki/File:Euler_method.svg">Wikipedia</a>.</figcaption></figure></p> |
<h2 id="zinc-uptake-model">Zinc Uptake Model</h2> | <h2 id="zinc-uptake-model">Zinc Uptake Model</h2> | ||
− | <p>Before adding | + | <p>Before adding our constructs to the simulator we first developed the models in <a href="http://uk.mathworks.com/products/matlab/">MATLAB</a>. Figure 2: This was done so that we could quickly get graphical output from our models and thereby assess their validity.</p> |
− | <p>Our <a href="https://2016.igem.org/Team:Newcastle/Parts#bio-varistor">variable resistor</a> works by varying the concentration of free zinc in an electrolyte in order to control the resistance. We wanted our model to reflect the zinc binding and release over time through the synthesis of <a href="http://www.rcsb.org/pdb/explore.do?structureId=1JJD">SmtA</a>, a metallothionein from <em>Cyanobacteria</em>. And to see the corresponding affect on the conductivity. For the electrolyte we chose to model this as zinc sulfate as this is a common electrolyte for which the conductivity data is <a href="http://sites.chem.colostate.edu/diverdi/all_courses/CRC%20reference%20data/electrical%20conductivity%20of%20aqueous%20solutions.pdf" | + | <p>Our <a href="https://2016.igem.org/Team:Newcastle/Parts#bio-varistor">variable resistor</a> works by varying the concentration of free zinc in an electrolyte in order to control the resistance. We wanted our model to reflect the zinc binding and release over time through the synthesis of <a href="http://www.rcsb.org/pdb/explore.do?structureId=1JJD">SmtA</a>, a metallothionein from <em>Cyanobacteria</em>. And to see the corresponding affect on the conductivity. For the electrolyte we chose to model this as zinc sulfate as this is a common electrolyte for which the conductivity data is well known <a href="http://sites.chem.colostate.edu/diverdi/all_courses/CRC%20reference%20data/electrical%20conductivity%20of%20aqueous%20solutions.pdf">(<em>CRC Handbook of Chemistry</em>, 1989)</a>. Because the conductivity is known for a variety of concentrations, we can relate this the number of free ions using the molecular weight of zinc sulfate (zinc sulfate decomposes to <script id="MathJax-Element-17" type="math/tex">Zn^{2+}</script> and <script id="MathJax-Element-18" type="math/tex">SO_{4}^{2-}</script> in solution). We used the available data from the literature to fit an equation relating conductivity to the number of free zinc ions. We determined that this relationship was linear and takes the form <script id="MathJax-Element-19" type="math/tex">conductivity = m \cdot N_{zinc\ molecules} + c</script>. Where <script id="MathJax-Element-20" type="math/tex">m</script> and <script id="MathJax-Element-21" type="math/tex">c</script> are parameters found by our equation fitting (<script id="MathJax-Element-22" type="math/tex">m = 3 \cdot 10^{-23}, c = 2 \cdot 10^{-15}</script>).</p> |
− | <p>At time <script id="MathJax-Element-23" type="math/tex">t_{0}</script> there is some amount of <script id="MathJax-Element-24" type="math/tex">Zn^{2+}</script> ions, in our model this is <script id="MathJax-Element-25" type="math/tex">0.2542M = 1.59613 \cdot 10^{23}</script>. Over time SmtA is expressed and binds to these Zinc molecules, <a href=" | + | <p>At time <script id="MathJax-Element-23" type="math/tex">t_{0}</script> there is some amount of <script id="MathJax-Element-24" type="math/tex">Zn^{2+}</script> ions, in our model this is <script id="MathJax-Element-25" type="math/tex">0.2542M = 1.59613 \cdot 10^{23}</script>. Over time, SmtA is expressed and binds to these Zinc molecules, four ions per protein <a href="#references">(<em>Blindauer et al.</em>, 2002)</a>. Similarly, as the protein is degraded the zinc is released again and can participate in conduction. Our model is a simplified view of the system as it ignores cellular mechanics, like the rate of zinc transport into and out of the cell. Our model treats the system as free SmtA binding to free zinc, with the cell serving only as a producer of protein. This leads us to hypothesise that the cell is unnecessary, on the condition that SmtA remains stable outside of the cell.</p> |
− | <p>From this observation we developed some <a href="https://2016.igem.org/Team:Newcastle/Parts#bio-varistor">additional variants of our variable resistor construct</a> that would allow us to purify the SmtA and examine this question experimentally. We also <a href="https://2016.igem.org/Team:Newcastle/Parts#bio-varistor-protein-only">investigated joining the SmtA protein to a LOV (light-oxygen-voltage) domain</a>. LOV domains are found plant, fungi and bacterial proteins which participate in phototropism or circadian rhythm setting. All LOV proteins have a blue-light sensitive flavin chromophore which means they can be used as light sensitive switches in hybrid proteins. This has previously been explored as | + | <p>From this observation we developed some <a href="https://2016.igem.org/Team:Newcastle/Parts#bio-varistor">additional variants of our variable resistor construct</a> that would allow us to purify the SmtA and examine this question experimentally. We also <a href="https://2016.igem.org/Team:Newcastle/Parts#bio-varistor-protein-only">investigated joining the SmtA protein to a LOV (light-oxygen-voltage) domain</a>. LOV domains are found plant, fungi and bacterial proteins which participate in phototropism or circadian rhythm setting. All LOV proteins have a blue-light sensitive flavin chromophore which means they can be used as light sensitive switches in hybrid proteins. This has previously been explored as LovTAP (<a href="https://2012.igem.org/Team:ETH_Zurich">ETH Zurich iGEM, 2012</a> and <a href="#references"><em>Strickland et al.</em>, 2008</a>). In our case we looked at using a LOV domain to create a light controlled zinc binding mechanism.</p> |
<p>Modelling zinc binding requires a small extension to main ODEs used in our simulator. We came up with another ODE to model the amount of bound zinc.</p> | <p>Modelling zinc binding requires a small extension to main ODEs used in our simulator. We came up with another ODE to model the amount of bound zinc.</p> | ||
− | |||
− | |||
<p><script id="MathJax-Element-26" type="math/tex; mode=display"> | <p><script id="MathJax-Element-26" type="math/tex; mode=display"> | ||
− | \frac{d(Bound\ zinc)}{dT} = (4 \cdot K_{zinc\ binding} \cdot Protein_{free}) - (4 \cdot k_{protein\ degradation} \cdot Protein_{bound}) | + | \frac{d(Bound\ zinc)}{dT} = (4 \cdot K_{zinc\ binding} \cdot Protein_{free}) - (4 \cdot k_{protein\ degradation} \cdot Protein_{bound}) \tag{6} |
</script></p> | </script></p> | ||
<p>The above equation is a variant of the protein production equation, <script id="MathJax-Element-27" type="math/tex">\frac{d(Protein)}{dT}</script>, which captures that every time step, the free SmtA binds with some efficiency to the free zinc ions. This efficiency is modeled as the constant <script id="MathJax-Element-28" type="math/tex">K_{zinc\ binding}</script> which is the percentage of the free proteins that bind each time step. Each protein that binds removes four zinc ions (we later chose to factor our this constant as a variable so that we could explore the impact that a better zinc binding protein would have). At the same time, some of the currently bound proteins are degraded, releasing four zinc ions each. This happens at the same rate as the protein degradation in the standard model.</p> | <p>The above equation is a variant of the protein production equation, <script id="MathJax-Element-27" type="math/tex">\frac{d(Protein)}{dT}</script>, which captures that every time step, the free SmtA binds with some efficiency to the free zinc ions. This efficiency is modeled as the constant <script id="MathJax-Element-28" type="math/tex">K_{zinc\ binding}</script> which is the percentage of the free proteins that bind each time step. Each protein that binds removes four zinc ions (we later chose to factor our this constant as a variable so that we could explore the impact that a better zinc binding protein would have). At the same time, some of the currently bound proteins are degraded, releasing four zinc ions each. This happens at the same rate as the protein degradation in the standard model.</p> | ||
− | <p>You’ll observe that at this point we’ve started to differentiate between bound and unbound proteins, but our model does not currently | + | <p>You’ll observe that at this point we’ve started to differentiate between bound and unbound proteins, but our model does not currently capture this. This inaccuracy prompted us to rethink our model, as by modelling the amount of bound zinc directly we had omitted an important biological process. This was most evident when we ran the model and found it to be unstable over large time steps (<script id="MathJax-Element-29" type="math/tex">\Delta t = 100</script>). While the bound zinc in this model followed an expected trend, closely following the amount of protein produced, it fluctuated between the trend and zero every 500 or so time steps. This was deemed undesirable, so we chose to adapt the model to include all intermediate steps of bound and unbound proteins. we achieved this feat using the same methodology to that of our original model that had zinc binding modelled directly.</p> |
<p><script id="MathJax-Element-30" type="math/tex; mode=display"> | <p><script id="MathJax-Element-30" type="math/tex; mode=display"> | ||
− | \frac{d(Protein_{bound})}{dT} = (K_{zinc\ binding} \cdot (Protein_{free} - Protein_{bound})) - (K_{protein\ degradation} \cdot Protein_{bound}) | + | \frac{d(Protein_{bound})}{dT} = (K_{zinc\ binding} \cdot (Protein_{free} - Protein_{bound})) - (K_{protein\ degradation} \cdot Protein_{bound}) \tag{7} |
</script></p> | </script></p> | ||
<p>The number of bound zinc ions can easily be found by multiplying the number of bound proteins by 4. This new model resulted in the expected behaviour, shown in Figure 3. Note that a larger time step was used to clearly show that the amount of free and bound protein diverge over time.</p> | <p>The number of bound zinc ions can easily be found by multiplying the number of bound proteins by 4. This new model resulted in the expected behaviour, shown in Figure 3. Note that a larger time step was used to clearly show that the amount of free and bound protein diverge over time.</p> | ||
− | <p><figure><img src="https://static.igem.org/mediawiki/2016/5/57/T--Newcastle--Zinc-binding-model-graph.png" alt="Zinc binding graph produced by our model." title=""><figcaption>Graph showing how free and bound zinc concentrations diverge over time in our model.</figcaption></figure></p> | + | <p><figure><img src="https://static.igem.org/mediawiki/2016/5/57/T--Newcastle--Zinc-binding-model-graph.png" alt="Zinc binding graph produced by our model." title=""><figcaption>Figure 3: Graph showing how free and bound zinc concentrations diverge over time in our model.</figcaption></figure></p> |
<p>When we connected this with our conductivity model, however, we noticed that the conductivity was not affected. At this point we realised that to date, we had been modelling only a single cell when we need to be modelling populations of cells to see our desired effects on the correct scale.</p> | <p>When we connected this with our conductivity model, however, we noticed that the conductivity was not affected. At this point we realised that to date, we had been modelling only a single cell when we need to be modelling populations of cells to see our desired effects on the correct scale.</p> | ||
− | <p>For our zinc model this is a simple matter of multiplying the single cell model to reflect the behavior of a whole population. | + | <p>For our zinc model this is a simple matter of multiplying the single cell model to reflect the behavior of a whole population. According to BioNumbers database there are approximately <script id="MathJax-Element-31" type="math/tex">10^{9}</script> cells per ml in an overnight culture (<a href="#references"><em>Milo et al.</em>, 2010</a>). As our chamber holds approximately 2 ml, we will use an estimated population size of of <script id="MathJax-Element-32" type="math/tex">2 \cdot 10^{9}</script> cells.</p> |
+ | |||
+ | <p>Using our population scale model we observe that the concentration of free zinc is not reduced enough to have any impact on the conductivity. We thought this might be a consequence of our model having too low protein expression, ~1400 proteins per cell. This would correspond to a cellular SmtA concentration of approx. 1.5 µM (according to <a href="http://bionumbers.hms.harvard.edu/KeyNumbers.aspx">BioNumbers</a>, 1 µM corresponds to ~1000 proteins in <em>Escherichia coli</em>). We had difficulty finding protein expression quantified in this way as expression levels are usually determined using a reporter protein like GFP and the level given in arbitrary fluorescence units (e.g. <a href="#references"><em>Leveau et al.</em></a>, 2001). The best data we could find was on the <a href="http://bionumbers.hms.harvard.edu//bionumber.aspx?id=112531&ver=2">concentration of CheZ</a> in a typical cell, since this is a natural <em>E. coli</em> protein (with a <a href="http://www.uniprot.org/uniprot/P0A9H9">role in chemotaxis</a>) rather than an engineered expression system the values may not correspond well with our expression system. It is however suitable for modelling purposes. The concentration of CheZ in a typical cell is ~4-12 µM; this equates to between 4000 and 12000 proteins which is higher than the values we currently use in our model.</p> | ||
+ | |||
+ | <p>Assuming that eventually all of the proteins bind the zinc, our population can remove <script id="MathJax-Element-33" type="math/tex">4 \cdot 12000 \cdot 2 \cdot 10^{9} = 9.6 * 10^{13}</script> zinc ions from the solution. This calculation suggests that we should have less zinc in our solution, or a greater cell density. If we wanted to be able to remove all the zinc we would need to start with a <script id="MathJax-Element-34" type="math/tex">{1.594×10^-10} mol / 2 ml = 0.792 \mu M</script> solution. This corresponds to a conductivity of <script id="MathJax-Element-35" type="math/tex">2.88 \cdot 10^{-9}</script> which represents a significant challenge to measure.</p> | ||
+ | |||
+ | <p>Because our model indicates limited success using this approach we used this feedback to explore other approaches to changing electrical conductivity. When we investigated alternative ways to change a mediums resistance we found that pH and conductivity (the reciprocal of resistivity) are intrinsically related. This is because pH is a measure of the number of H+ ions in solution and it is these ions, among others, that conduct current. The relationship between pH and H+ concentration is <script id="MathJax-Element-99" type="math/tex">pH = -log[H+]</script>. Therefore, lower pH solutions have both a higher concentration of H+ ions and higher conductivity.</p> | ||
+ | |||
+ | <p>Altering the resistance in this way has the advantage that even if we were unable to measure resistance we knew we could still integrate this with electronics. We would be able to do this using readily available <a href="https://www.sparkfun.com/products/10972">pH sensors</a> which we knew about due to our discussions with <a href="2016.igem.org/Team:Newcastle/Collaborations#iwbda">Luis Ortiz at IWBDA 2016</a> who was working with this technique.</p> | ||
+ | |||
+ | <p>We adapted the basic genetic toggle switch design of <a href="#references"><em>Gardner et al</em>. (2000)</a> as shown in figure 4 by adding <span class="gene-name">lacZ</span> as the reporter (figure 5). When <span class="gene-name">lacZ</span> is expressed the E. coli can metabolise lactose in their growth media, this results in acid production which lowers the pH of the growth media and thereby also lowers the resistance.</p> | ||
− | < | + | <figure><img id="togglefig" src="https://static.igem.org/mediawiki/2016/a/a5/T--Newcastle--ToggleSwitchFigure1-Gardneretal.png" /><figcaption>Figure 4: The basic toggle switch of <a href="#references"><em>Gardner et al</em>. (2000)</a>. Repressor 1 inhibits transcription from Promoter 1 and is induced by Inducer 1. Repressor 2 inhibits transcription from Promoter 2 and is induced by Inducer 2.</figcaption></figure> |
− | < | + | <figure><img src="https://static.igem.org/mediawiki/2016/3/33/T--Newcastle--phToggleConstruct.png" /><figcaption>Figure 5: SBOL Visual Diagram of the Design of our pH Toggle Switch Device.</figcaption></figure> |
− | <p> | + | <p>To toggle the expression of <span class="gene-name">lacZ</span> we chose a TetR/LacI scheme as this had been shown to be effective in previous iGEM projects like <a href="https://2014.igem.org/Team:Wageningen_UR/project/kill-switch">Wageningen 2014</a>. In this scheme the addition of aTC to the medium inhibits the off signal (tetR expression) causing the on signal (<span class="gene-name">lacI</span>) and reporter to be expressed, this is the 'on' state. The addition of IPTG on the other hand causes the repression of the on signal (<span class="gene-name">lacI</span>), allowing expression of the off signal (<span class="gene-name">tetR</span>) which inhibits the further expression of the reporter. This is the off state. This two state system is shown in figure 5. In the future we wold like to further augment this design with a mechanism for reducing the pH in the off state, for example through <a href="https://2006.igem.org/wiki/index.php/Urease">Urease</a> production. However this remains a challenge as the Urease protein is very complex (5 kb of DNA encoding 7 genes).</p> |
<h2 id="biological-capacitor-model">Biological Capacitor Model</h2> | <h2 id="biological-capacitor-model">Biological Capacitor Model</h2> | ||
− | <p>As with all our models the biological capacitor is built around the core equations for transcription and translation. For our <a href="https://2016.igem.org/Team:Newcastle/Parts">biological capacitor</a> the only extension upon our existing model base is that instead of keeping track of one protein we have to keep track of three: lambda repressor, 434 repressor and a reporter protein. In our construct, we chose GFP as our reporter gene. However, in the final construct, this protein would ideally be the protein that causes the capacitor to charge so that they could be chained together. The following set of ODEs mathematically describe the system:</p> | + | <p>As with all our models the biological capacitor is built around the core equations for transcription and translation. For our <a href="https://2016.igem.org/Team:Newcastle/Parts">biological capacitor</a> the only extension upon our existing model base is that instead of keeping track of one protein we have to keep track of three: lambda repressor, 434 repressor, and a reporter protein. In our construct, we chose GFP as our reporter gene. However, in the final construct, this protein would ideally be the protein that causes the capacitor to charge so that they could be chained together. The following set of ODEs mathematically describe the system:</p> |
<p><script id="MathJax-Element-1936" type="math/tex; mode=display"> | <p><script id="MathJax-Element-1936" type="math/tex; mode=display"> | ||
Line 166: | Line 208: | ||
\frac{d(Protein_{Lambda})}{dT} = K_{translation} \cdot mRNA - K_{deg\ protein} \cdot Protein_{Lambda}\\ | \frac{d(Protein_{Lambda})}{dT} = K_{translation} \cdot mRNA - K_{deg\ protein} \cdot Protein_{Lambda}\\ | ||
\frac{d(Protein_{Output})}{dT} = K_{translation} \cdot mRNA - K_{deg\ protein} \cdot Protein_{Output} | \frac{d(Protein_{Output})}{dT} = K_{translation} \cdot mRNA - K_{deg\ protein} \cdot Protein_{Output} | ||
+ | \tag{8} | ||
</script></p> | </script></p> | ||
Line 174: | Line 217: | ||
<p>As with our previous models, there are transcription and translation rate constants that directly control the amount of protein synthesised. By using the same strong promoter values as in our other models, we get the standard behaviour. Namely, that the number of proteins in the cell gradually increases up to some carrying capacity at which point they degrade at the same rate as they are produced and the amount of protein in the cell levels off.</p> | <p>As with our previous models, there are transcription and translation rate constants that directly control the amount of protein synthesised. By using the same strong promoter values as in our other models, we get the standard behaviour. Namely, that the number of proteins in the cell gradually increases up to some carrying capacity at which point they degrade at the same rate as they are produced and the amount of protein in the cell levels off.</p> | ||
− | <p><figure><img src="https://static.igem.org/mediawiki/2016/7/70/T--Newcastle--Simulator-Capacitor-Std-Behaviour.png" alt="Graph showing behaviour of biological capacitor if left to charge indefinitely." title=""><figcaption>Figure | + | <p><figure><img src="https://static.igem.org/mediawiki/2016/7/70/T--Newcastle--Simulator-Capacitor-Std-Behaviour.png" alt="Graph showing behaviour of biological capacitor if left to charge indefinitely." title=""><figcaption>Figure 6: Graph showing the expected 'levelling-off' behaviour present in our capacitor model.</figcaption></figure></p> |
− | <p>A simulation that simply accumulates proteins is not of significant biological interest | + | <p>A simulation that simply accumulates proteins is not of significant biological interest. Thus we needed a way to model the addition of our repressor L-arabinose. We observed that fundamentally, inducers and repressors alter the promoter strength. Since there will always be leaky constitutive production, the addition of an inducer simple increases the strength, and the addition of a repressor decreases the strength of the promoter.</p> |
− | <p>In our model because 434 and lambda repressor are under the control of the same promoter they effectively share a <em>transscripton rate</em> which corresponds to the strength of that promoter. Consequently to simulate the addition of a repressor in our model we need only vary the value of this transcriptional constant. Figure | + | <p>In our model because 434 and lambda repressor are under the control of the same promoter they effectively share a <em>transscripton rate</em> which corresponds to the strength of that promoter. Consequently to simulate the addition of a repressor in our model we need only vary the value of this transcriptional constant. Figure 6 shows the behaviour of the system when the promoter is weakened at <script type="math/tex">t=4000</script>. We changed the promoter from <em>strong</em> to <em>weak</em> using the values outlined in the ‘simulator core’ section.</p> |
− | <p><figure><img src="https://static.igem.org/mediawiki/2016/0/01/T--Newcastle--Simulator-Capacitor-Rept4000.png" alt="Graph showing simulation of capacitor system behaviour when repressor is added at time t=4000." title=""><figcaption>Figure | + | <p><figure><img src="https://static.igem.org/mediawiki/2016/0/01/T--Newcastle--Simulator-Capacitor-Rept4000.png" alt="Graph showing simulation of capacitor system behaviour when repressor is added at time t=4000." title=""><figcaption>Figure 7: Behaviour of the system when the promoter is weakened at t=4000.</figcaption></figure></p> |
− | <p>As can be seen in Figure | + | <p>As can be seen in Figure 7, the 434 repressor quickly decays to the point where there is more lambda repressor than 434 repressor. At this stage, we would expect our output protein to be produced as 434 repressor would out compete lambda repressor and prevent transcription of the output protein.</p> |
− | <p>We tried different approaches to modelling this change. Our first method used a stochastic approach. This approach modelled the competition between lambda repressor and 434 repressor to induce/repress production (i.e. transcription and translation) of our output protein. In this variant of the model we used a ‘lottery’ system to decide the transcriptional rate of the output protein. For each time step, if the 434 repressor ‘won’, then the promoter is ‘weak’; otherwise the promoter is ‘strong’ | + | <p>We tried different approaches to modelling this change. Our first method used a stochastic approach. This approach modelled the competition between lambda repressor and 434 repressor to induce/repress production (i.e. transcription and translation) of our output protein. In this variant of the model we used a ‘lottery’ system to decide the transcriptional rate of the output protein. For each time step, if the 434 repressor ‘won’, then the promoter is ‘weak’; otherwise the promoter is ‘strong’. The way the lottery system works is that each protein is given a number of ‘tickets’ in direct proportion to the total amount of protein that exists. A ticket is then chosen at random and if it belongs to 434 repressor, then it ‘wins’, otherwise lambda wins. This can easily be simulated using the proportion of each protein to allocate tickets on the line 0-1 and then using a random number from 0 to 1 to determine the winner. This results in the following behaviour in the system.</p> |
− | <p><figure><img src="https://static.igem.org/mediawiki/2016/d/d3/T--Newcastle--Simulator-capacitor-stochastic.png" alt="Graph showing simulated behaviour of biological capacitor using stochastic approach." title=""><figcaption>Figure | + | <p><figure><img src="https://static.igem.org/mediawiki/2016/d/d3/T--Newcastle--Simulator-capacitor-stochastic.png" alt="Graph showing simulated behaviour of biological capacitor using stochastic approach." title=""><figcaption>Figure 8: Graph showing simulated behaviour of biological capacitor using stochastic approach.</figcaption></figure></p> |
<p>As can be seen in Figure 6, the stochastic model results in leaky behaviour. This type of behaviour was not the intended behaviour of our construct. The output signal follows the same trend as the lambda and 434 repressor, but with a large amount of noise. In our construct, we would expect that the expression of the output signal would remain low until the amount of 434 repressor was greater than the amount of lambda repressor at which point there would be a spike in the output signal which would decay as the repressors do. We think that this is similar to having the promoter for the output be ‘weak’, while having it be ‘strong’ when there is more 434 than lambda repressor. It is easy to adapt the model to this behaviour, as this is the same technique we use to model the addition of L-arabinose.</p> | <p>As can be seen in Figure 6, the stochastic model results in leaky behaviour. This type of behaviour was not the intended behaviour of our construct. The output signal follows the same trend as the lambda and 434 repressor, but with a large amount of noise. In our construct, we would expect that the expression of the output signal would remain low until the amount of 434 repressor was greater than the amount of lambda repressor at which point there would be a spike in the output signal which would decay as the repressors do. We think that this is similar to having the promoter for the output be ‘weak’, while having it be ‘strong’ when there is more 434 than lambda repressor. It is easy to adapt the model to this behaviour, as this is the same technique we use to model the addition of L-arabinose.</p> | ||
− | <p><figure><img src="https://static.igem.org/mediawiki/2016/e/eb/T--Newcastle--Simulator-cap-correct.png" alt="Graph showing the correct simulated behaviour of our biological capacitor construct." title=""><figcaption>Figure | + | <p><figure><img src="https://static.igem.org/mediawiki/2016/e/eb/T--Newcastle--Simulator-cap-correct.png" alt="Graph showing the correct simulated behaviour of our biological capacitor construct." title=""><figcaption>Figure 9: Graph showing the correct simulated behaviour of our biological capacitor construct.</figcaption></figure></p> |
<p>Figure 7 shows the desired behaviour of our system, which confirms that our construct needs to work as described above to exhibit the correct behaviour. If on the other hand, the interaction between the lambda and 434 repressor behaves as in the stochastic simulation then our construct will not work as expected. Although we expect our construct to behave as in this latter simulation we will need to verify this empirically.</p> | <p>Figure 7 shows the desired behaviour of our system, which confirms that our construct needs to work as described above to exhibit the correct behaviour. If on the other hand, the interaction between the lambda and 434 repressor behaves as in the stochastic simulation then our construct will not work as expected. Although we expect our construct to behave as in this latter simulation we will need to verify this empirically.</p> | ||
Line 196: | Line 239: | ||
<h3 id="circuit-analysis">Circuit Analysis</h3> | <h3 id="circuit-analysis">Circuit Analysis</h3> | ||
− | <p>Because our components are to be integrated into electrical circuits, we needed a way of modelling how their resistances would affect other circuit components as well as a means of determining the current though, and voltage across, each of the components as this effects the heat | + | <p>Because our components are to be integrated into electrical circuits, we needed a way of modelling how their resistances would affect other circuit components as well as a means of determining the current though, and voltage across, each of the components as this effects the heat dissipated in the component as a result of the Joule effect.</p> |
− | <p>After some research we settled on using <a href="https://en.wikipedia.org/wiki/Superposition_theorem">‘analysis by superposition’</a> as our circuit analysis technique. This allows us to have | + | <p>After some research we settled on using <a href="https://en.wikipedia.org/wiki/Superposition_theorem">‘analysis by superposition’</a> as our circuit analysis technique. This allows us to have multiple voltage sources in our circuit which is important as we want to be able to mix microbial fuel cell sources and standard DC sources like batteries and see the effect. Circuit analysis by superposition considers each voltage source at a time (‘shorting’ out the others) and then summing the resulting voltage values to get the final voltage drop for all components.</p> |
− | <p>Before we could perform our analysis we needed to decide on our representation of our circuits. For | + | <p>Before we could perform our analysis we needed to decide on our representation of our circuits. For our components, we have chosen to consider our systems as a series of resistors and sources, grounded at either end. This creates a linear circuit like the figure below.</p> |
− | <p><figure><img src="https://static.igem.org/mediawiki/2016/2/27/T--Newcastle--CircuitModel.png" alt="Linear circuit model" title=""><figcaption>Figure | + | <p><figure><img src="https://static.igem.org/mediawiki/2016/2/27/T--Newcastle--CircuitModel.png" alt="Linear circuit model" title=""><figcaption>Figure 10: Illustration of the circuit level abstraction we are using in our model.</figcaption></figure></p> |
− | <p>Computationally we build our circuits as linked lists for easy traversal | + | <p>Computationally we build our circuits as linked lists for easy traversal and simple editing operations. Additionally, the node structure of the linked list allows the system to generalise to parallel circuits, which would be represented by a graph data structure if at a later date we choose to implement this functionality.</p> |
<p>The algorithm for analysing such circuits is simple, and is illustrated in the worked example below.</p> | <p>The algorithm for analysing such circuits is simple, and is illustrated in the worked example below.</p> | ||
Line 214: | Line 257: | ||
<p>Consider the circuit shown here.</p> | <p>Consider the circuit shown here.</p> | ||
− | <p><figure><img src="https://static.igem.org/mediawiki/2016/f/f6/T--Newcastle--WorkedCircuitModel.png" alt="Circuit diagram used in the worked example." title=""><figcaption>Figure | + | <p><figure><img src="https://static.igem.org/mediawiki/2016/f/f6/T--Newcastle--WorkedCircuitModel.png" alt="Circuit diagram used in the worked example." title=""><figcaption>Figure 11: The circuit used in the worked example.</figcaption></figure></p> |
<p>Step 1: Short each voltage source in turn and compute the voltage across the resistors.</p> | <p>Step 1: Short each voltage source in turn and compute the voltage across the resistors.</p> | ||
Line 234: | Line 277: | ||
<p>We checked this analysis by modelling the same circuit in <a href="http://www.falstad.com/circuit/circuitjs.html">an existing circuit simulator</a> as shown here, where we got the same results.</p> | <p>We checked this analysis by modelling the same circuit in <a href="http://www.falstad.com/circuit/circuitjs.html">an existing circuit simulator</a> as shown here, where we got the same results.</p> | ||
− | <p><figure><img src="https://static.igem.org/mediawiki/2016/b/b3/T--Newcastle--Existing-Sim-Circuit.png" alt="Diagram of example circuit modelled using existing software." title=""><figcaption>Figure | + | <p><figure><img src="https://static.igem.org/mediawiki/2016/b/b3/T--Newcastle--Existing-Sim-Circuit.png" alt="Diagram of example circuit modelled using existing software." title=""><figcaption>Figure 12: Our worked example circuit (Figure 11) modelled using existing software, <a href="http://www.falstad.com/circuit/circuitjs.html">CircuitJS</a> by <a href="http://lushprojects.com/circuitjs/">Iain Sharp</a>.</figcaption></figure></p> |
<p>You can also see the same circuit in our very own simulator implementing this technique.</p> | <p>You can also see the same circuit in our very own simulator implementing this technique.</p> | ||
− | <p><figure><img src="https://static.igem.org/mediawiki/2016/6/66/T--Newcastle--Our-Sim-Circuit.png" alt="Diagram of example circuit modelled using our software." title=""><figcaption>Figure | + | <p><figure><img src="https://static.igem.org/mediawiki/2016/6/66/T--Newcastle--Our-Sim-Circuit.png" alt="Diagram of example circuit modelled using our software." title=""><figcaption>Figure 13: Our worked example circuit (Figure 11) modelled using our worked software.</figcaption></p> |
<h4 id="temperature-change">Temperature Change</h4> | <h4 id="temperature-change">Temperature Change</h4> | ||
− | <p>As some of our <a href="https://2016.igem.org/Team:Newcastle/Parts#bio-bulb">parts are based on exploiting the the heat-shock response of <em>E. coli</em></a> it is important that we be able to determine the temperature change over time for our components. Using what we learnt in our physics modelling using <a href="https://www.comsol.com/">COMSOL</a> we knew we could approximate the temperature rise from the | + | <p>As some of our <a href="https://2016.igem.org/Team:Newcastle/Parts#bio-bulb">parts are based on exploiting the the heat-shock response of <em>E. coli</em></a> it is important that we be able to determine the temperature change over time for our components. Using what we learnt in our physics modelling using <a href="https://www.comsol.com/">COMSOL</a> we knew we could approximate the temperature rise from the Joule heating effect using the power dissipated in the component and the thermal conductivity of the media (which we determined with the <a href="https://2016.igem.org/Team:Newcastle/Collaborations">help of the Exter iGEM team</a>).</p> |
<p>For our models we use <script id="MathJax-Element-2250" type="math/tex">P = I^{2}R</script> and <script id="MathJax-Element-2251" type="math/tex">\Delta T = \frac{P \cdot L}{2 \cdot K\cdot A}</script>. Where K is the thermal conductivity (for LB broth this is <script id="MathJax-Element-2252" type="math/tex">\approx 0.58</script>), L is the length and A is the area of the container.</p> | <p>For our models we use <script id="MathJax-Element-2250" type="math/tex">P = I^{2}R</script> and <script id="MathJax-Element-2251" type="math/tex">\Delta T = \frac{P \cdot L}{2 \cdot K\cdot A}</script>. Where K is the thermal conductivity (for LB broth this is <script id="MathJax-Element-2252" type="math/tex">\approx 0.58</script>), L is the length and A is the area of the container.</p> | ||
Line 249: | Line 292: | ||
<p>We also looked at using <script id="MathJax-Element-2253" type="math/tex">Q = m \cdot c \cdot \Delta T</script> where <script id="MathJax-Element-2254" type="math/tex">Q = Power \cdot Time</script> and <script id="MathJax-Element-2255" type="math/tex">c</script> is the specfic heat capacity, which we know for LB is similar to water (<script id="MathJax-Element-2256" type="math/tex">=4200</script>). Both of these equations give similar order of magnitude approximations. Unfortunately they are not complex enough to capture aspects liek the cooling effect of the surrounding air on the components which limites the accuracy of our model somewhat and can lead to incorrect data in the simulator.</p> | <p>We also looked at using <script id="MathJax-Element-2253" type="math/tex">Q = m \cdot c \cdot \Delta T</script> where <script id="MathJax-Element-2254" type="math/tex">Q = Power \cdot Time</script> and <script id="MathJax-Element-2255" type="math/tex">c</script> is the specfic heat capacity, which we know for LB is similar to water (<script id="MathJax-Element-2256" type="math/tex">=4200</script>). Both of these equations give similar order of magnitude approximations. Unfortunately they are not complex enough to capture aspects liek the cooling effect of the surrounding air on the components which limites the accuracy of our model somewhat and can lead to incorrect data in the simulator.</p> | ||
− | <p><figure><img src="https://static.igem.org/mediawiki/2016/f/fd/T--Newcastle--COMSOL-Chamber.png" alt="Multiphysics temperature model of our microfluidics chambers." title=""><figcaption>Figure | + | <p><figure><img src="https://static.igem.org/mediawiki/2016/f/fd/T--Newcastle--COMSOL-Chamber.png" alt="Multiphysics temperature model of our microfluidics chambers." title=""><figcaption>Figure 14: Multiphysics temperature model of our microfluidics chambers.</figcaption></figure></p> |
− | <p>Given more time we would have tried to include some of these aspects in our simulator or found a way to include data from external programs like COMSOL which we were able to use to perform this more | + | <p>Given more time we would have tried to include some of these aspects in our simulator or found a way to include data from external programs like COMSOL which we were able to use to perform this more sophisticated analysis.</p> |
<h3 id="the-ui">The UI</h3> | <h3 id="the-ui">The UI</h3> | ||
Line 257: | Line 300: | ||
<p>We chose to build a UI for our models using HTML so that other people could explore the ways our technology would work in an easy way. Ideally, through our iGEM wiki. This is an important part of being a foundational advance team, showing how our technologies could be used.</p> | <p>We chose to build a UI for our models using HTML so that other people could explore the ways our technology would work in an easy way. Ideally, through our iGEM wiki. This is an important part of being a foundational advance team, showing how our technologies could be used.</p> | ||
− | <p>To speed up development our UI was constructed around the existing JavaScript models, which we ported from the above MATLAB models and tied to an integrator which is always running in the background. We used a number of libraries to aid our development, including <a href="https://jqueryui.com/draggable/">jQuery UI</a> to allow for drag and drop circuit construction, the <a href="https://picturepan2.github.io/spectre/">Spectre.css</a> CSS framework for basic page layout and <a href="http://www.chartjs.org/">Chart.js</a> to show the output from our models.</p> | + | <p>To speed up development our UI was constructed around the existing JavaScript models, which we ported from the above MATLAB models and tied to an integrator which is always running in the background. We used a number of libraries to aid our development, including <a href="https://jqueryui.com/draggable/">jQuery UI</a> to allow for drag and drop circuit construction, the <a href="https://picturepan2.github.io/spectre/">Spectre.css</a> CSS framework for basic page layout and, <a href="http://www.chartjs.org/">Chart.js</a> to show the output from our models.</p> |
− | <p><figure><img src="https://static.igem.org/mediawiki/2016/a/a9/T--Newcastle--Our-Sim-Circuit2.png" alt="Screenshot of bacterial lightbulb and power in our simulator." title=""><figcaption>Figure | + | <p><figure><img src="https://static.igem.org/mediawiki/2016/a/a9/T--Newcastle--Our-Sim-Circuit2.png" alt="Screenshot of bacterial lightbulb and power in our simulator." title=""><figcaption>Figure 15: Screenshot of bacterial lightbulb and power source in our simulator.</figcaption></p> |
<p>You can <a href="https://2016.igem.org/Team:Newcastle/Software/Hosted/Models">try out more circuits</a> in our editor to see the effect of using a bacterial resistor among other things.</p> | <p>You can <a href="https://2016.igem.org/Team:Newcastle/Software/Hosted/Models">try out more circuits</a> in our editor to see the effect of using a bacterial resistor among other things.</p> | ||
Line 267: | Line 310: | ||
<h4 id="references">References</h4> | <h4 id="references">References</h4> | ||
− | < | + | <ol id="reflist"> |
− | Leveau, J.H. and Lindow, S.E., 2001. Predictive and interpretive simulation of green fluorescent protein expression in reporter bacteria. Journal of bacteriology, 183(23), pp.6752-6762. < | + | <li>Blindauer, C. A., Harrison, M. D., Robinson, A. K., Parkinson, J. A., Bowness, P. W., Sadler, P. J. and Robinson, N. J. (2002), <a href="http://onlinelibrary.wiley.com/doi/10.1046/j.1365-2958.2002.03109.x/full">Multiple bacteria encode metallothioneins and SmtA-like zinc fingers</a>. Molecular Microbiology, 45: 1421–1432. doi:10.1046/j.1365-2958.2002.03109.x</li> |
− | Milo, R., Jorgensen, P., Moran, U., Weber, G. and Springer, M., 2010. BioNumbers—the database of key numbers in molecular and cell biology. Nucleic acids research, 38(suppl 1), pp.D750-D753. < | + | <li>Gardner, T.S., Cantor, C.R. and Collins, J.J., 2000. <a href="http://www.nature.com/nature/journal/v403/n6767/full/403339a0.html">Construction of a genetic toggle switch in Escherichia coli</a>. Nature, 403(6767), pp.339-342.</li> |
− | Strickland, D., Moffat, K. and Sosnick, T.R., 2008. Light-activated DNA binding in a designed allosteric protein. Proceedings of the National Academy of Sciences, 105(31), pp.10709-10714.</ | + | <li>Leveau, J.H. and Lindow, S.E., 2001. <a href="http://jb.asm.org/content/183/23/6752.full">Predictive and interpretive simulation of green fluorescent protein expression in reporter bacteria</a>. Journal of bacteriology, 183(23), pp.6752-6762.</li> |
+ | <li>Milo, R., Jorgensen, P., Moran, U., Weber, G. and Springer, M., 2010. <a href="http://kirschner.med.harvard.edu/files/bionumbers/fundamentalBioNumbersHandout.pdf">BioNumbers—the database of key numbers in molecular and cell biology</a>. Nucleic acids research, 38(suppl 1), pp.D750-D753.</li> | ||
+ | <li>Strickland, D., Moffat, K. and Sosnick, T.R., 2008. <a href="http://www.pnas.org/content/105/31/10709.figures-only">Light-activated DNA binding in a designed allosteric protein</a>. Proceedings of the National Academy of Sciences, 105(31), pp.10709-10714.</li> | ||
+ | </ol> | ||
</div> | </div> | ||
</body> | </body> | ||
</html> | </html> |
Latest revision as of 18:38, 19 October 2016
Modelling
This page documents the dry lab experiments relating to our modelling.
ODE Modelling
As part of our work, we performed many different types of modelling, from rule-based modelling of our genetic parts to physics analysis of our microfluidics components. On this page, we discuss the ordinary differential equation (ODE) modelling of our components for use in a web-based ‘simulator’. We hope that our ‘simulator’ will demonstrate how our genetic constructs could fit into an electronic circuit.
Simulator Core
The core of our simulator features a set of ODEs that describe an abstract version of the transcriptional (Equation 1) and translational (Equation 2) processes in a cell.
The rate at which these processes occur is controlled by a number of paramters, and . These constants describe the rate of transcription rate, degradation of mRNA, translation rate and protein degradation rate respectively. In our modelling we focus only on variation in the and parameters; we can fluctuate these values various promoter strengths and protein degradation tags. The other parameters are unrelated to any of the elements of our system that we control. Consequently, we have assigned these fixed values in our simulator.
Since our code is freely available, a modeller with different simulation requirements could easily change these parameters if they were meaningful in their scenario. We also constrain the other parameters in our simulator so that the simulation is stable when the user is given control over these values. For example, it might be logical for the user to select between a weak, medium or strong promoter. In this scenario, the user should not be required to know the exact numerical values of their choice. Although more advanced users may wish to be able to input their own values. For this reason, our simulation uses sane defaults for the transcription rate, , which is analogous to promoter strength.
In our design work, users thought being constrained to these three options limited their ability to fully explore the biological system. For example, it would not allow the user to see the variation in gene product caused by slight noise in the promoter strength value, e.g. switching between 0.03 and 0.1 quickly. As our system is primarily a way of exploring the constructs we have built we wanted to limit the end user as much as possible. Therefore we set out plans to replace the notion of strong/medium/weak promoters with a more fine-grained scale. One that the user might control with a simple slider for example. To achieve this feat we fit our existing values using regression to a 10 point scale. This gave us an equation that we can use to map a user’s idea of ‘promoter strength’ to a transcription rate. As shown in Figure 1, the resulting equation was . Note that we deliberately ensured that to reflect the fact that there is always some background level transcription. In layman's terms, no biological promoter can be completely ‘switched off’.
The remainder of our simulator is built around these simple equations. The rest of this page will discuss how the other components were simulated using this technique and then integrated into the simulator. Our constructs were found to be easily adaptable to this form of model. However, it is not applicable in all situations. Our equations capture a simple ‘PRET’ (Promoter-Ribosome-Expression-Terminator) construct. As such it can be generalized to any number of protein expressions by having a tuple for each gene product in the model and separate counts of mRNA and protein.
Justification
We chose to build our simulator around a set of ODEs rather than using a more complex system that produced models in a format like CellML or SBML because it met our design goals. Primarily, this was ease of development and integration into the iGEM wiki system. Our main aim was to allow the simulator to be embedded in the wiki. In turn we felt that it was necessary to support the aims of the iGEM wiki system in preventing ‘bit-rot’, that is not having some aspects of our project hosted somewhere else where their uptime could not be guaranteed. This meant we needed a purely client side implementation and couldn’t ‘call out’ to other computers to perform the computation. When researching how we would construct our system we determined that the most straightforward of modelling our system in JavaScript without using any libraries or complex dependencies was to hard-code the ODEs.
Another aspect that makes the ODEs approach useful for us is that the level of abstraction they capture i.e., transcription and translation, corresponds directly with the level at which we design our BioBrick constructs. This makes it easy to translate between formats and allows for us to rapidly add new knowledge to the simulator as we develop our constructs. And vice versa, the simulations at this level translate directly to components in the system which we can control in our constructs. For example, we can directly influence the number of proteins in the cell by increasing the number of coding sequences on a construct. Whereas we can’t influence variables like the amount of free mRNA through our constructs alone.
There are some trade-offs to this approach. For instance, we lose access to tools that perform local and global sensitivity analysis, parameter estimation, etc., which limits the ‘depth’ of our models. We also lose out on interoperability with future software through interchange standards like SMBL. However, we believe that this doesn’t subtract from the quality of our models because their intended purpose is to encourage the user to think about the implications of interfacing bacteria with electronics rather than providing detailed models. In other words, our present systems are not complex enough to warrant modelling at a lower level of abstraction in the simulator.
That is not to say that we didn’t also perform lower level modelling. We also used rule-based modelling in the BioNetGen language using the RuleBender graphical tool. This allowed us to capture our systems in more detail and to ‘sanity check’ the simplified models we use here in our simulator.
Integration
There are a number of ways to integrate the differential equations presented above. The simplest is to iteratively compute the difference and add these to a ‘running total’ over the lifetime of the simulation. That is, for each time step compute where h is the step size, f is a differential equation such as those above. This method is known as Euler’s method. It can be visualised as fitting a series of straight lines to the given differential equation. It is, therefore, more accurate at small step sizes (the global error is proportional to the step size and the local error is proportional to the square of the step size). This is the method we use in our simulator.
.
Zinc Uptake Model
Before adding our constructs to the simulator we first developed the models in MATLAB. Figure 2: This was done so that we could quickly get graphical output from our models and thereby assess their validity.
Our variable resistor works by varying the concentration of free zinc in an electrolyte in order to control the resistance. We wanted our model to reflect the zinc binding and release over time through the synthesis of SmtA, a metallothionein from Cyanobacteria. And to see the corresponding affect on the conductivity. For the electrolyte we chose to model this as zinc sulfate as this is a common electrolyte for which the conductivity data is well known (CRC Handbook of Chemistry, 1989). Because the conductivity is known for a variety of concentrations, we can relate this the number of free ions using the molecular weight of zinc sulfate (zinc sulfate decomposes to and in solution). We used the available data from the literature to fit an equation relating conductivity to the number of free zinc ions. We determined that this relationship was linear and takes the form . Where and are parameters found by our equation fitting ().
At time there is some amount of ions, in our model this is . Over time, SmtA is expressed and binds to these Zinc molecules, four ions per protein (Blindauer et al., 2002). Similarly, as the protein is degraded the zinc is released again and can participate in conduction. Our model is a simplified view of the system as it ignores cellular mechanics, like the rate of zinc transport into and out of the cell. Our model treats the system as free SmtA binding to free zinc, with the cell serving only as a producer of protein. This leads us to hypothesise that the cell is unnecessary, on the condition that SmtA remains stable outside of the cell.
From this observation we developed some additional variants of our variable resistor construct that would allow us to purify the SmtA and examine this question experimentally. We also investigated joining the SmtA protein to a LOV (light-oxygen-voltage) domain. LOV domains are found plant, fungi and bacterial proteins which participate in phototropism or circadian rhythm setting. All LOV proteins have a blue-light sensitive flavin chromophore which means they can be used as light sensitive switches in hybrid proteins. This has previously been explored as LovTAP (ETH Zurich iGEM, 2012 and Strickland et al., 2008). In our case we looked at using a LOV domain to create a light controlled zinc binding mechanism.
Modelling zinc binding requires a small extension to main ODEs used in our simulator. We came up with another ODE to model the amount of bound zinc.
The above equation is a variant of the protein production equation, , which captures that every time step, the free SmtA binds with some efficiency to the free zinc ions. This efficiency is modeled as the constant which is the percentage of the free proteins that bind each time step. Each protein that binds removes four zinc ions (we later chose to factor our this constant as a variable so that we could explore the impact that a better zinc binding protein would have). At the same time, some of the currently bound proteins are degraded, releasing four zinc ions each. This happens at the same rate as the protein degradation in the standard model.
You’ll observe that at this point we’ve started to differentiate between bound and unbound proteins, but our model does not currently capture this. This inaccuracy prompted us to rethink our model, as by modelling the amount of bound zinc directly we had omitted an important biological process. This was most evident when we ran the model and found it to be unstable over large time steps (). While the bound zinc in this model followed an expected trend, closely following the amount of protein produced, it fluctuated between the trend and zero every 500 or so time steps. This was deemed undesirable, so we chose to adapt the model to include all intermediate steps of bound and unbound proteins. we achieved this feat using the same methodology to that of our original model that had zinc binding modelled directly.
The number of bound zinc ions can easily be found by multiplying the number of bound proteins by 4. This new model resulted in the expected behaviour, shown in Figure 3. Note that a larger time step was used to clearly show that the amount of free and bound protein diverge over time.
When we connected this with our conductivity model, however, we noticed that the conductivity was not affected. At this point we realised that to date, we had been modelling only a single cell when we need to be modelling populations of cells to see our desired effects on the correct scale.
For our zinc model this is a simple matter of multiplying the single cell model to reflect the behavior of a whole population. According to BioNumbers database there are approximately cells per ml in an overnight culture (Milo et al., 2010). As our chamber holds approximately 2 ml, we will use an estimated population size of of cells.
Using our population scale model we observe that the concentration of free zinc is not reduced enough to have any impact on the conductivity. We thought this might be a consequence of our model having too low protein expression, ~1400 proteins per cell. This would correspond to a cellular SmtA concentration of approx. 1.5 µM (according to BioNumbers, 1 µM corresponds to ~1000 proteins in Escherichia coli). We had difficulty finding protein expression quantified in this way as expression levels are usually determined using a reporter protein like GFP and the level given in arbitrary fluorescence units (e.g. Leveau et al., 2001). The best data we could find was on the concentration of CheZ in a typical cell, since this is a natural E. coli protein (with a role in chemotaxis) rather than an engineered expression system the values may not correspond well with our expression system. It is however suitable for modelling purposes. The concentration of CheZ in a typical cell is ~4-12 µM; this equates to between 4000 and 12000 proteins which is higher than the values we currently use in our model.
Assuming that eventually all of the proteins bind the zinc, our population can remove zinc ions from the solution. This calculation suggests that we should have less zinc in our solution, or a greater cell density. If we wanted to be able to remove all the zinc we would need to start with a solution. This corresponds to a conductivity of which represents a significant challenge to measure.
Because our model indicates limited success using this approach we used this feedback to explore other approaches to changing electrical conductivity. When we investigated alternative ways to change a mediums resistance we found that pH and conductivity (the reciprocal of resistivity) are intrinsically related. This is because pH is a measure of the number of H+ ions in solution and it is these ions, among others, that conduct current. The relationship between pH and H+ concentration is . Therefore, lower pH solutions have both a higher concentration of H+ ions and higher conductivity.
Altering the resistance in this way has the advantage that even if we were unable to measure resistance we knew we could still integrate this with electronics. We would be able to do this using readily available pH sensors which we knew about due to our discussions with Luis Ortiz at IWBDA 2016 who was working with this technique.
We adapted the basic genetic toggle switch design of Gardner et al. (2000) as shown in figure 4 by adding lacZ as the reporter (figure 5). When lacZ is expressed the E. coli can metabolise lactose in their growth media, this results in acid production which lowers the pH of the growth media and thereby also lowers the resistance.
![](https://static.igem.org/mediawiki/2016/a/a5/T--Newcastle--ToggleSwitchFigure1-Gardneretal.png)
![](https://static.igem.org/mediawiki/2016/3/33/T--Newcastle--phToggleConstruct.png)
To toggle the expression of lacZ we chose a TetR/LacI scheme as this had been shown to be effective in previous iGEM projects like Wageningen 2014. In this scheme the addition of aTC to the medium inhibits the off signal (tetR expression) causing the on signal (lacI) and reporter to be expressed, this is the 'on' state. The addition of IPTG on the other hand causes the repression of the on signal (lacI), allowing expression of the off signal (tetR) which inhibits the further expression of the reporter. This is the off state. This two state system is shown in figure 5. In the future we wold like to further augment this design with a mechanism for reducing the pH in the off state, for example through Urease production. However this remains a challenge as the Urease protein is very complex (5 kb of DNA encoding 7 genes).
Biological Capacitor Model
As with all our models the biological capacitor is built around the core equations for transcription and translation. For our biological capacitor the only extension upon our existing model base is that instead of keeping track of one protein we have to keep track of three: lambda repressor, 434 repressor, and a reporter protein. In our construct, we chose GFP as our reporter gene. However, in the final construct, this protein would ideally be the protein that causes the capacitor to charge so that they could be chained together. The following set of ODEs mathematically describe the system:
Since our construct contains twice as many coding sequences for 434 repressor as lambda repressor its concentration increases twice as fast which is reflected in the above equations. By virtue of it growing twice as fast, it also degrades twice as fast because of the equations above. This is the desired behaviour since in our construct we use a ssRA degradation tag to have the 434 repressor degrade faster. In reality it may not be the case that the degradation tag results in a 2x speed up in degradation but this can be accomodated in the model. We would do this by having seperate translational constants such that and where captures the relative degradation ratio of lambda to 434 repressor and could be determined experimentally.
Our biological capacitor was constructed to mimic the behaviour of an electrical capacitor. Rather than accumulating charge, the capacitor accumulates proteins that are used to drive an output when the input signal stops being applied. For demonstrative purposes of showing we can partially mimic this behaviour, we assume that the signal is applied constitutively and then switched off by the addition of a chemical. In our construct this is L-arabinose, but from a modelling perspective, all that needs to change is the promoter strength.
As with our previous models, there are transcription and translation rate constants that directly control the amount of protein synthesised. By using the same strong promoter values as in our other models, we get the standard behaviour. Namely, that the number of proteins in the cell gradually increases up to some carrying capacity at which point they degrade at the same rate as they are produced and the amount of protein in the cell levels off.
A simulation that simply accumulates proteins is not of significant biological interest. Thus we needed a way to model the addition of our repressor L-arabinose. We observed that fundamentally, inducers and repressors alter the promoter strength. Since there will always be leaky constitutive production, the addition of an inducer simple increases the strength, and the addition of a repressor decreases the strength of the promoter.
In our model because 434 and lambda repressor are under the control of the same promoter they effectively share a transscripton rate which corresponds to the strength of that promoter. Consequently to simulate the addition of a repressor in our model we need only vary the value of this transcriptional constant. Figure 6 shows the behaviour of the system when the promoter is weakened at . We changed the promoter from strong to weak using the values outlined in the ‘simulator core’ section.
As can be seen in Figure 7, the 434 repressor quickly decays to the point where there is more lambda repressor than 434 repressor. At this stage, we would expect our output protein to be produced as 434 repressor would out compete lambda repressor and prevent transcription of the output protein.
We tried different approaches to modelling this change. Our first method used a stochastic approach. This approach modelled the competition between lambda repressor and 434 repressor to induce/repress production (i.e. transcription and translation) of our output protein. In this variant of the model we used a ‘lottery’ system to decide the transcriptional rate of the output protein. For each time step, if the 434 repressor ‘won’, then the promoter is ‘weak’; otherwise the promoter is ‘strong’. The way the lottery system works is that each protein is given a number of ‘tickets’ in direct proportion to the total amount of protein that exists. A ticket is then chosen at random and if it belongs to 434 repressor, then it ‘wins’, otherwise lambda wins. This can easily be simulated using the proportion of each protein to allocate tickets on the line 0-1 and then using a random number from 0 to 1 to determine the winner. This results in the following behaviour in the system.
As can be seen in Figure 6, the stochastic model results in leaky behaviour. This type of behaviour was not the intended behaviour of our construct. The output signal follows the same trend as the lambda and 434 repressor, but with a large amount of noise. In our construct, we would expect that the expression of the output signal would remain low until the amount of 434 repressor was greater than the amount of lambda repressor at which point there would be a spike in the output signal which would decay as the repressors do. We think that this is similar to having the promoter for the output be ‘weak’, while having it be ‘strong’ when there is more 434 than lambda repressor. It is easy to adapt the model to this behaviour, as this is the same technique we use to model the addition of L-arabinose.
Figure 7 shows the desired behaviour of our system, which confirms that our construct needs to work as described above to exhibit the correct behaviour. If on the other hand, the interaction between the lambda and 434 repressor behaves as in the stochastic simulation then our construct will not work as expected. Although we expect our construct to behave as in this latter simulation we will need to verify this empirically.
Circuit Analysis
Because our components are to be integrated into electrical circuits, we needed a way of modelling how their resistances would affect other circuit components as well as a means of determining the current though, and voltage across, each of the components as this effects the heat dissipated in the component as a result of the Joule effect.
After some research we settled on using ‘analysis by superposition’ as our circuit analysis technique. This allows us to have multiple voltage sources in our circuit which is important as we want to be able to mix microbial fuel cell sources and standard DC sources like batteries and see the effect. Circuit analysis by superposition considers each voltage source at a time (‘shorting’ out the others) and then summing the resulting voltage values to get the final voltage drop for all components.
Before we could perform our analysis we needed to decide on our representation of our circuits. For our components, we have chosen to consider our systems as a series of resistors and sources, grounded at either end. This creates a linear circuit like the figure below.
Computationally we build our circuits as linked lists for easy traversal and simple editing operations. Additionally, the node structure of the linked list allows the system to generalise to parallel circuits, which would be represented by a graph data structure if at a later date we choose to implement this functionality.
The algorithm for analysing such circuits is simple, and is illustrated in the worked example below.
Circuit Analysis Example
Consider the circuit shown here.
Step 1: Short each voltage source in turn and compute the voltage across the resistors.
Shorting the 1.5V source leaves only the 9V source, the voltage across each resistor can then be found using the voltage divider formula: .
Therefore for R1 we have, . And for R2 we have, .
We can check our answers make sense by seeing that they sum to , ✓.
Shorting the 9V source leaves the 1.5V source. The voltage across R1 is . Note that the voltage is positive even though the source is after the resistor because it is ‘dragging up’ the ground. Across R2 the resistance is .
Step 2: Sum the total voltage across each resistor.
For R1 the total voltage drop is and for R2 it is . Again we can check that our total voltage drop matches the voltage supplied in the circuit, ✓.
Step 3: Finally we can derive any of the other data we need using the voltage drop through each component along with other equations. For example we can compute the current using Ohm’s law, . That is, the current through R1 is and through R2 is . As you would expect the current is the same all the way through the circuit.
We checked this analysis by modelling the same circuit in an existing circuit simulator as shown here, where we got the same results.
You can also see the same circuit in our very own simulator implementing this technique.
Temperature Change
As some of our parts are based on exploiting the the heat-shock response of E. coli it is important that we be able to determine the temperature change over time for our components. Using what we learnt in our physics modelling using COMSOL we knew we could approximate the temperature rise from the Joule heating effect using the power dissipated in the component and the thermal conductivity of the media (which we determined with the help of the Exter iGEM team).
For our models we use and . Where K is the thermal conductivity (for LB broth this is ), L is the length and A is the area of the container.
We also looked at using where and is the specfic heat capacity, which we know for LB is similar to water (). Both of these equations give similar order of magnitude approximations. Unfortunately they are not complex enough to capture aspects liek the cooling effect of the surrounding air on the components which limites the accuracy of our model somewhat and can lead to incorrect data in the simulator.
Given more time we would have tried to include some of these aspects in our simulator or found a way to include data from external programs like COMSOL which we were able to use to perform this more sophisticated analysis.
The UI
We chose to build a UI for our models using HTML so that other people could explore the ways our technology would work in an easy way. Ideally, through our iGEM wiki. This is an important part of being a foundational advance team, showing how our technologies could be used.
To speed up development our UI was constructed around the existing JavaScript models, which we ported from the above MATLAB models and tied to an integrator which is always running in the background. We used a number of libraries to aid our development, including jQuery UI to allow for drag and drop circuit construction, the Spectre.css CSS framework for basic page layout and, Chart.js to show the output from our models.
You can try out more circuits in our editor to see the effect of using a bacterial resistor among other things.
References
- Blindauer, C. A., Harrison, M. D., Robinson, A. K., Parkinson, J. A., Bowness, P. W., Sadler, P. J. and Robinson, N. J. (2002), Multiple bacteria encode metallothioneins and SmtA-like zinc fingers. Molecular Microbiology, 45: 1421–1432. doi:10.1046/j.1365-2958.2002.03109.x
- Gardner, T.S., Cantor, C.R. and Collins, J.J., 2000. Construction of a genetic toggle switch in Escherichia coli. Nature, 403(6767), pp.339-342.
- Leveau, J.H. and Lindow, S.E., 2001. Predictive and interpretive simulation of green fluorescent protein expression in reporter bacteria. Journal of bacteriology, 183(23), pp.6752-6762.
- Milo, R., Jorgensen, P., Moran, U., Weber, G. and Springer, M., 2010. BioNumbers—the database of key numbers in molecular and cell biology. Nucleic acids research, 38(suppl 1), pp.D750-D753.
- Strickland, D., Moffat, K. and Sosnick, T.R., 2008. Light-activated DNA binding in a designed allosteric protein. Proceedings of the National Academy of Sciences, 105(31), pp.10709-10714.