Difference between revisions of "Team:Oxford/Modeling"

 
(48 intermediate revisions by 5 users not shown)
Line 2: Line 2:
  
 
<html>
 
<html>
−
     <div class="container-fluid page-heading">
+
<style>
−
         <h3>Modeling</h3>
+
/*main colour*/
 +
    .navbar-default .navbar-brand,
 +
    .contents-sidebar .nav>.active>a,
 +
    .contents-sidebar .nav>.active>a,
 +
    .contents-sidebar .nav>li>a:hover,
 +
    .contents-sidebar .nav>li>a:focus,
 +
    h2,
 +
    h1,
 +
    .main-container .section a {
 +
        color: #1B46DD;
 +
    }
 +
    .contents-sidebar .nav>.active>a,
 +
    .contents-sidebar .nav>li>a:hover,
 +
    .contents-sidebar .nav>li>a:focus {
 +
        border-left: 2px solid #1B46DD;
 +
    }
 +
    .image.lightbox, #notebook-key-button {
 +
        background-color: #1B46DD;
 +
    }
 +
 
 +
    /*complimentary colour*/
 +
    .navbar-default .navbar-brand:hover,
 +
    .definition:hover, .definition:focus,
 +
    ol li::before,
 +
    .slim ul li:before,
 +
    .table>thead>tr>th,
 +
    .algorithm ol li::before,
 +
    .quote,
 +
    .quote h3,
 +
    a,
 +
    a:visited,
 +
    a:hover
 +
    {
 +
        color:#32C97F;
 +
    }
 +
    .definition {
 +
        border-bottom: 1px dotted #32C97F;
 +
    }
 +
    #notebook-key-button.active {
 +
        background-color: #32C97F;
 +
    }
 +
    .popover.right>.arrow::after{
 +
        border-right-color: #32C97F;
 +
    }
 +
    .popover.bottom>.arrow::after {
 +
        border-bottom-color: #32C97F;
 +
    }
 +
    .popover.left>.arrow::after {
 +
        border-left-color: #32C97F;
 +
    }
 +
</style>
 +
     <div class="container-fluid page-heading" style="background-image:url(https://static.igem.org/mediawiki/2015/a/ab/Ox_modelling-background.jpg)">
 +
         <h3>Modelling</h3>
 
     </div>
 
     </div>
 
     <div class="container-fluid">
 
     <div class="container-fluid">
 
         <div class="row">
 
         <div class="row">
 
             <div class="col-md-9">
 
             <div class="col-md-9">
−
                 <div class="section" id="introduction">
+
                 <div class="slim">
−
                     <div class="slim">
+
                     <div class="section" id="introduction">
 
                         <h2>Introduction</h2>
 
                         <h2>Introduction</h2>
 
                         <p>
 
                         <p>
−
                             Modeling allows us to observe what happens to a system when you change the parameters that describe it. By quantifying the system’s response to these variables, we can aid the wet lab team in their design of their biological system. By repeating the modeling, data-fitting and refining schedule we approach a sufficient description of our required design.
+
                             Mathematical <a class="definition" title="model" data-content="A simplified or idealised description of a system or process, usually mathematical, that can be used to predict how it will behave.">modelling</a> plays a crucial role in Synthetic Biology by acting as a link between the conception and the physical realisation of a biological circuit. Our modelling team has evaluated the effectiveness of initial designs, and has provided insight into how the system can (or must) be improved.
 
                         </p>
 
                         </p>
 
                         <p>
 
                         <p>
−
                             Our modeling is made up of three pieces: gene expression, diffusion and biofilm modeling. By modeling gene expression, we can inform the lab team which parameters they can choose freely - and which ones they must pick carefully - to maximise the output of our product. The gene copy number is one example of these parameters. By modeling diffusion, we can incorporate into our previous model the effect of our protein being separated from the biofilm it is trying to interact with. If this is a large effect, the lab team must change the design of their system. Finally, by modeling the formation of a biofilm and its interaction with our product, we can decide upon the required concentrations of product we need to clear the urinary tract infection. In all cases we use both deterministic and stochastic models to investigate our system.
+
                             Our team experimentally validated that <em>Escherichia coli</em> can secrete enzymes which break down the biofilms associated with urinary infections. However, it is difficult to directly measure whether our enzymes are produced in a sufficient quantity to be a more effective treatment than antibiotics. We measured gene expression and diffusion of widely-used chemicals, and then used our model to estimate the number of <em>E. coli</em> cells that would make our project a more effective treatment than antibiotics. We expect to have to improve our system to make it realistic.
 +
                        </p>
 +
                        <p>
 +
                            To help readers of all kinds and specialisations understand this page we have produced guides for all the modelling techniques used in this section. They are available in our Modelling Tutorial page and will be linked to when appropriate.
 
                         </p>
 
                         </p>
 
                     </div>
 
                     </div>
−
                </div>
+
                    <div class="section-spacer"></div>
−
                <div class="section-spacer"></div>
+
                    <div class="section" id="characterising-our-cells">
−
                <div class="section" id="gene-expression">
+
                         <h2>Gene expression rates</h2>
−
                    <div class="slim">
+
                        <p>
−
                         <h2>Gene Expression</h2>
+
                            In this section we look at our cells in isolation in order to assess their functionality and answer important questions such as “how long does it take to produce a certain concentration of product?” The end result - the final concentration of useful enzyme that is produced in the cell - is required for our diffusion model.
−
                         <div id="gene-expression-deterministic">
+
                        </p>
−
                             <h3>Deterministic</h3>
+
                         <div id="characterising-our-cells-arab">
 +
                             <h3>Arabinose-induced expression</h3>
 
                             <p>
 
                             <p>
−
                                 A deterministic model evolves according to ordinary or partial differential equations. We model the following system:
+
                                 We have decided to use an <a class="definition" title="arabinose" data-content="A sugar which is commonly used to induce gene expression.">arabinose</a>-<a class="definition" title="induced" data-content="When we say expression is induced, we mean that it only happens when the inducer is present.">induced</a> <a class="definition" title="promoter" data-content="The section of DNA that the RNA polymerase enzyme binds to before it starts making the RNA strand - it is needed to start transcription, so it sits in the DNA before a gene. Regulates whether a gene is &quot;on&quot; or &quot;off&quot; and to what extent. They can be made to be sensitive to certain conditions so that if a bacterium senses a change in environment it can up or down regulate the expression of a certain gene (so there will be more or less of the protein encoded by that gene, produced in the cell).">promoter</a> for the <a class="definition" title="gene expression" data-content="The production of a protein within the cell.">expression</a> of a number of our <a class="definition" title="protein" data-content="An essential part of all living organisms. They are long and fold up into complicated structures and are made up of amino acids.">proteins</a>. This promoter can be modelled as the following chemical system:
 
                             </p>
 
                             </p>
−
                                \[\overset{\alpha_{1}}{\rightarrow}Arab\quad\overset{\alpha_{5}}{\rightarrow}AraC\]
 
  
−
                                \[Arab+AraC\mathrel{\mathop{\rightleftharpoons}^{\mathrm{\alpha_{2}}}_{\mathrm{\alpha_{3}}}}Arab:AraC\]
+
                            \[(Arab:AraC)\overset{K}{\rightarrow}mRNA\overset{\alpha}{\rightarrow}P\]
−
 
+
                            \[mRNA\overset{\gamma_{1}}{\rightarrow}\phi\quad P\overset{\gamma_{2}}{\rightarrow}\phi\]
−
                                \[\overset{K}{\rightarrow}mRNA\overset{\alpha_{4}}{\rightarrow}P\]
+
−
 
+
−
                                \[Arab\overset{\gamma_{1}}{\rightarrow}\phi\quad AraC\overset{\gamma_{2}}{\rightarrow}\phi\quad mRNA\overset{\gamma_{3}}{\rightarrow}\phi\quad P\overset{\gamma_{4}}{\rightarrow}\phi\]
+
 
                             <p>
 
                             <p>
−
                                 Using the following deterministic equations [<a href="#Ref4">4</a>]:
+
                                 Our promoter, pBAD, binds to <a class="definition" title="AraC" data-content="A compound which stops our E. coli from producing enzymes">AraC</a> and this represses <a class="definition" title="transcription" data-content="The process of converting DNA to mRNA.">transcription</a> of <a class="definition" title="mRNA" data-content="Messenger RNA (mRNA) carries the information of the gene we wish to express to the ribosome. The protein is built at the ribosome.">mRNA</a>. Arabinose will bind to AraC and form the Arab:AraC compound, allowing transcription to occur.
 +
                            </p>
 +
                            <p>
 +
                                For this system we will assume that AraC is always in large concentration and that its binding to arabinose happens on a faster time scale to transcription. Therefore, we do not need to consider the individual concentrations of arabinose and AraC, instead we just need to include the concentration of the complex (Arab:AraC). The rate \(K\) is not just a simple constant and is given as the <a class="definition" title="Hill function" data-content="The Hill function looks like a smooth step on a stairwell. Depending on the Hill coefficient n, the step can be made to be sharper or smoother.">Hill function</a> in the equations below.
 +
                            </p>
 +
                            <p>
 +
                                Using <a href="https://2015.igem.org/Team:Oxford/Modeling/Tutorial#gene-expression-networks-michaelis">Michaelis-Mentin kinetics</a>, we arrive at the equations:
 
                             </p>
 
                             </p>
−
                                \[\dfrac{d\left[Arab\right]}{dt}=\alpha_{1}+\alpha_{2}\left[Arab\right]\left[AraC\right]-\alpha_{3}\left[Arab:AraC\right]-\gamma_{1}\left[Arab\right]\]
 
  
−
                                 \[\dfrac{d\left[Arab:AraC\right]}{dt}=\alpha_{3}\left[Arab:AraC\right]-\alpha_{2}\left[Arab\right]\left[AraC\right]\]
+
                                 \[\dfrac{d[mRNA]}{dt}=K_{max}\dfrac{[Arab:AraC]^{n}}{K_{half}^{n}+[Arab:AraC]^{n}}-\gamma_{1}[mRNA]\]
  
−
                                \[\dfrac{d\left[AraC\right]}{dt}=\alpha_{5}+\alpha_{2}\left[Arab\right]\left[AraC\right]-\alpha_{3}\left[Arab:AraC\right]-\gamma_{2}\left[AraC\right]\]
+
                                 \[\dfrac{d\left[P\right]}{dt}=\alpha\left[mRNA\right]-\gamma_{2}\left[P\right]\]
−
 
+
−
                                \[\dfrac{d[mRNA]}{dt}=K_{max}\dfrac{[Arab:AraC]^{n}}{K_{half}^{n}+[Arab:AraC]^{n}}-\gamma_{3}[mRNA]\]
+
−
 
+
−
                                 \[\dfrac{d\left[P\right]}{dt}=\alpha_{4}\left[mRNA\right]-\gamma_{4}\left[P\right]\]
+
 
                             <p>
 
                             <p>
−
                                 Where we define the symbols
+
                                 Where \([Arab]\), \([AraC]\), \([Arab:AraC]\), \([mRNA]\) and \([P]\) represent the concentrations of arabinose, AraC, Arab:AraC, mRNA and our product protein respectively. We define the remaining symbols in the table below.
−
                                 <table border="1">
+
                                 <table class="table table-striped">
−
                                     <tr>
+
                                     <thead>
 
                                         <th>Symbol</th>
 
                                         <th>Symbol</th>
 
                                         <th>Definition</th>
 
                                         <th>Definition</th>
 
                                         <th>Initial Value/Literature Value</th>
 
                                         <th>Initial Value/Literature Value</th>
 
                                         <th>Fitted</th>
 
                                         <th>Fitted</th>
−
                                     </tr>
+
                                     </thead>
 
                                     <tr>
 
                                     <tr>
−
                                         <td>\([Arab]\)</td>
+
                                         <td>\(\alpha\)</td>
−
                                        <td>The concentration of Arabinose</td>
+
                                         <td><a class="definition" title="translation" data-content="The process of converting mRNA to a protein. Occurs at the ribosome.">Translation</a> rate</td>
−
                                        <td>\(1\times10^{-4}M\)</td>
+
                                         <td>\(15ntd\: s^{-1}\)/length of sequence [<a href="#references">6</a>]</td>
−
                                         <td>-</td>
+
                                         <td>\(6.6ntd\: s^{-1}\)/length of sequence</td>
−
                                    </tr>
+
−
                                    <tr>
+
−
                                        <td>\([AraC]\)</td>
+
−
                                        <td>The concentration of AraC</td>
+
−
                                        <td>\(1\times10^{-5}M\)</td>
+
−
                                        <td>-</td>
+
−
                                    </tr>
+
−
                                    <tr>
+
−
                                        <td>\([Arab:AraC]\)</td>
+
−
                                        <td>The concentration of associated Arabinose and AraC</td>
+
−
                                        <td>\(0\)</td>
+
−
                                        <td>-</td>
+
−
                                    </tr>
+
−
                                    <tr>
+
−
                                        <td>\([mRNA]\)</td>
+
−
                                        <td>The concentration of mRNA</td>
+
−
                                        <td>\(0\)</td>
+
−
                                        <td>-</td>
+
−
                                    </tr>
+
−
                                    <tr>
+
−
                                        <td>\([P]\)</td>
+
−
                                        <td>The concentration of our product</td>
+
−
                                        <td>\(0\)</td>
+
−
                                        <td>-</td>
+
−
                                    </tr>
+
−
                                    <tr>
+
−
                                        <td>\(\alpha_{1}\)</td>
+
−
                                        <td>Basal production of Arabinose</td>
+
−
                                        <td>???</td>
+
−
                                        <td>?</td>
+
−
                                    </tr>
+
−
                                    <tr>
+
−
                                        <td>\(\alpha_{2}\)</td>
+
−
                                        <td>Association constant</td>
+
−
                                        <td>\(2.8\times10^{7}s^{-1}\) [<a href="#Ref7">7</a>]</td>
+
−
                                        <td>?</td>
+
−
                                    </tr>
+
−
                                    <tr>
+
−
                                        <td>\(\alpha_{3}\)</td>
+
−
                                        <td>Dissociation constant</td>
+
−
                                        <td>\(0.022s^{-1}\) [<a href="#Ref7">7</a>]</td>
+
−
                                        <td>?</td>
+
−
                                    </tr>
+
−
                                    <tr>
+
−
                                        <td>\(\alpha_{4}\)</td>
+
−
                                        <td>Translation rate</td>
+
−
                                         <td>\(15ntd\: s^{-1}\)/length of sequence [<a href="#Ref6">6</a>]</td>
+
−
                                        <td>?</td>
+
−
                                    </tr>
+
−
                                    <tr>
+
−
                                         <td>\(\alpha_{5}\)</td>
+
−
                                        <td>Basal production of AraC</td>
+
−
                                        <td>???</td>
+
−
                                        <td>?</td>
+
 
                                     </tr>
 
                                     </tr>
 
                                     <tr>
 
                                     <tr>
 
                                         <td>\(\gamma_{1}\)</td>
 
                                         <td>\(\gamma_{1}\)</td>
−
                                         <td> Degradation rate of Arabinose </td>
+
                                         <td><a class="definition" title="combining the rates" data-content="Reference 5 links to a cell division time of 22.5 minutes. Reference 10 links to a mRNA half life of 6.8 minutes. The combined rate is therefore given by ln(2)/22.5/60+ln(2)/6.8/60.">Combined</a> <a class="definition" title="degradation rate" data-content="Proteins and mRNA are unstable and will decay into other products. The degradation rate tells us how quickly this process happens.">degradation</a> and <a class="definition" title="dilution rate" data-content="Bacterial cells divide up to every half hour. Each time this happens, the amount of mRNA or protein present in the cell halves.">dilution</a> rate of mRNA</td>
−
                                         <td>\(5.13\times10^{-4}s^{-1}\) [<a href="#Ref5">5</a>]</td>
+
                                         <td>\(2.2\times10^{-3}s^{-1}\) [<a href="#references">5</a>, <a href="#references">10</a>]</td>
−
                                         <td>?</td>
+
                                         <td>\(1.1\times10^{-2}s^{-1}\)</td>
 
                                     </tr>
 
                                     </tr>
 
                                     <tr>
 
                                     <tr>
 
                                         <td>\(\gamma_{2}\)</td>
 
                                         <td>\(\gamma_{2}\)</td>
−
                                         <td>Degradation rate of AraC</td>
+
                                         <td><a class="definition" title="combining the rates" data-content="Reference 5 links to a cell division time of 22.5 minutes. Reference 11 links to a GFP half life of 24 hours.">Combined</a> degradation and dilution rate of <a class="definition" title="GFP" data-content="Green fluorescent protein. It is relatively easy to measure the amount of green light that this protein produces.">GFP</a></td>
−
                                        <td>\(5.13\times10^{-4}s^{-1}\) [<a href="#Ref5">5</a>]</td>
+
                                         <td>\(5.2\times10^{-4}s^{-1}\) [<a href="#references">5</a>, <a href="#references">11</a>]</td>
−
                                        <td>?</td>
+
                                         <td>\(1.1\times10^{-2}s^{-1}\)</td>
−
                                    </tr>
+
−
                                    <tr>
+
−
                                        <td>\(\gamma_{3}\)</td>
+
−
                                        <td>Degradation rate of mRNA</td>
+
−
                                         <td>\(5.13\times10^{-4}s^{-1}\) [<a href="#Ref5">5</a>]</td>
+
−
                                        <td>?</td>
+
−
                                    </tr>
+
−
                                    <tr>
+
−
                                        <td>\(\gamma_{4}\)</td>
+
−
                                        <td>Degradation rate of product</td>
+
−
                                         <td>\(5.13\times10^{-4}s^{-1}\) [<a href="#Ref5">5</a>]</td>
+
−
                                        <td>?</td>
+
 
                                     </tr>
 
                                     </tr>
 
                                     <tr>
 
                                     <tr>
 
                                         <td>\(K_{max}\)</td>
 
                                         <td>\(K_{max}\)</td>
−
                                         <td>Maximal transcription rate</td>
+
                                         <td>Maximal <a class="definition" title="transcription" data-content="The synthesis of mRNA from DNA.">transcription</a> rate</td>
−
                                         <td>\(50ntd\: s^{-1}\)/length of sequence [<a href="#Ref6">6</a>]</td>
+
                                         <td>\(50ntd\: s^{-1}\)/length of sequence [<a href="#references">6</a>]</td>
−
                                         <td>?</td>
+
                                         <td>\(47ntd\: s^{-1}\)/length of sequence</td>
 
                                     </tr>
 
                                     </tr>
 
                                     <tr>
 
                                     <tr>
 
                                         <td>\(K_{half}\)</td>
 
                                         <td>\(K_{half}\)</td>
 
                                         <td>Half-maximal transcription rate</td>
 
                                         <td>Half-maximal transcription rate</td>
−
                                         <td>\(160\mu M\) [<a href="#Ref8">8</a>]</td>
+
                                         <td>\(160\mu M\) [<a href="#references">7</a>]</td>
−
                                         <td>?</td>
+
                                         <td>\(100\mu M\)</td>
 
                                     </tr>
 
                                     </tr>
 
                                     <tr>
 
                                     <tr>
 
                                         <td>\(n\)</td>
 
                                         <td>\(n\)</td>
−
                                         <td>Hill coefficient</td>
+
                                         <td><a class="definition" title="Hill coefficient" data-content="A measure of the smoothness of the graph of gene expression.">Hill coefficient</a></td>
−
                                         <td>\(2.65\) [<a href="#Ref3">3</a>]</td>
+
                                         <td>\(2.65\) [<a href="#references">8</a>]</td>
−
                                         <td>?</td>
+
                                         <td>\(2.73\)</td>
 
                                     </tr>
 
                                     </tr>
 
                                 </table>
 
                                 </table>
 
                             </p>
 
                             </p>
 
                             <p>
 
                             <p>
−
                                 To evaluate these parameters, we will fit the data created by the wet lab team. We allow a tolerance of a factor of ten either side of our initial guess for each parameter. The fitted parameters will be given in the fourth column and will be found using the MATLAB function fmincon. To solve the equations numerically we use the MATLAB equation solver ode15s. We plot three solutions for various parameters below.
+
                                 This table contains literature values for the parameters, found from a number of sources. We then measured <a class="definition" title="GFP" data-content="Green fluorescent protein. It is relatively easy to measure the amount of green light that this protein produces.">GFP</a> expression in <em>E. coli</em> to extract experimental values. Here is a plot showing the fit of this model to our experimental data.
 
                             </p>
 
                             </p>
 
                             <div class="image image-full">
 
                             <div class="image image-full">
−
                                 <img src="https://static.igem.org/mediawiki/2015/7/74/OxiGEM_Modeling_Arabinose_01.png"/>
+
                                 <img src="https://static.igem.org/mediawiki/2015/d/de/OxiGEM_Gene_Fitter.png" alt="Fitting our gene expression data to the theoretical model" />
 +
                                <p>
 +
                                    Results showing GFP concentration as a function of time, matched to our deterministic model. Errors are given to one standard deviation and an arbitrary scaling factor is included as a fitted parameter.
 +
                                </p>
 
                             </div>
 
                             </div>
−
                        </div>
 
−
                        <div id="gene-expression-stochastic">
 
−
                            <h3>Stochastic</h3>
 
 
                             <p>
 
                             <p>
−
                                 A stochastic model is anything which includes some random element, and is useful for modeling small systems or those with some uncertainty. We produced a video tutorial on how to model gene expression networks stochastically using Gillespie's algorithm.[<a href="#Ref2">2</a>, <a href="#Ref4">4</a>].
+
                                 We can now calculate the limiting concentrations that our products will be expressed. The dilution rate and Hill coefficient of our cells is the same for GFP and our proteins, but the transcription and translation rates are dependent on the sequence length of the protein. Here is a table showing the relevant proteins and sequence lengths:
 
                             </p>
 
                             </p>
−
                             <div class="algorithm">
+
                             <table class="table table-striped">
−
                                 <h3>The Gillespie Algorithm</h3>
+
                                 <thead>
−
                                 <ol>
+
                                    <tr>
−
                                     <li>Find random numbers \(r_{1}\) and \(r_{2}\) which are uniformly distributed between 0 and 1.</li>
+
                                        <th>
−
 
+
                                            Product
−
                                     <li>Find propensities \(\alpha_{i}\) of \(N\) possible reactions. Find \(\alpha_{0}=\sum\limits_{i=1}^N \alpha_{i}\).</li>
+
                                        </th>
−
 
+
                                        <th>
−
                                     <li>Compute time \(\tau=\dfrac{1}{\alpha_{0}}\ln\left(\dfrac{1}{r_{1}}\right)\).</li>
+
                                            Sequence Length (/bp)
−
 
+
                                        </th>
−
                                     <li>Find the reaction \(j\) such that \(\dfrac{1}{\alpha_{0}}\sum\limits_{i=1}^{j-1} \alpha_{i}\leq r_{2}<\dfrac{1}{\alpha_{0}}\sum\limits_{i=1}^{j} \alpha_{i}\).</li>
+
                                    </tr>
−
 
+
                                </thead>
−
                                     <li>Implement reaction \(j\) and increment time \(t\) to \(t+\tau\). </li>
+
                                 <tr>
−
 
+
                                     <td>
−
                                     <li>Repeat steps 1-5 until a sufficient time or computation budget has been exhausted.</li>
+
                                        pBAD HisB DNase DsbA
−
                                 </ol>
+
                                    </td>
−
 
+
                                    <td>
 +
                                        621
 +
                                     </td>
 +
                                </tr>
 +
                                <tr>
 +
                                     <td>
 +
                                        pBAD HisB MccS
 +
                                    </td>
 +
                                    <td>
 +
                                        414
 +
                                     </td>
 +
                                </tr>
 +
                                <tr>
 +
                                    <td>
 +
                                        pBAD HisB Art-175 DsbA
 +
                                    </td>
 +
                                    <td>
 +
                                        987
 +
                                    </td>
 +
                                </tr>
 +
                                <tr>
 +
                                     <td>
 +
                                        pBAD HisB Art-175 YebF
 +
                                    </td>
 +
                                    <td>
 +
                                        1284
 +
                                     </td>
 +
                                </tr>
 +
                                <tr>
 +
                                    <td>
 +
                                        pBAD HisB Art-E
 +
                                    </td>
 +
                                    <td>
 +
                                        632
 +
                                    </td>
 +
                                </tr>
 +
                                <tr>
 +
                                    <td>
 +
                                        pBAD HisB Art-175 Fla
 +
                                    </td>
 +
                                    <td>
 +
                                        1095
 +
                                    </td>
 +
                                </tr>
 +
                                <tr>
 +
                                    <td>
 +
                                        pBAD HisB Art-175
 +
                                    </td>
 +
                                    <td>
 +
                                        936
 +
                                    </td>
 +
                                </tr>
 +
                                <tr>
 +
                                    <td>
 +
                                        pBAD HisB DNase
 +
                                    </td>
 +
                                    <td>
 +
                                        570
 +
                                    </td>
 +
                                </tr>
 +
                            </table>
 +
                            <p>
 +
                                We now can run our model of the system by solving the set of equations using the MATLAB <a class="definition" title="ODE" data-content="Ordinary differential equation. Describes the small-scale changes of variables.">ordinary differential equation</a> solver ode15s. Below is a plot of the concentration of product against time for each protein expressed with this inducer-promoter pair where the expression is induced by a step function:
 +
                            </p>
 +
                            <div class="image image-full">
 +
                                 <img src="https://static.igem.org/mediawiki/2015/f/f6/Ox_arab_induced_proteins.png"/>
 +
                                <p>
 +
                                    Model data for each of the enzymes we plan to release, using the parameters we found from our experimental data. We found our limiting concentrations were of order \(\mu M\).
 +
                                </p>
 
                             </div>
 
                             </div>
 
                             <p>
 
                             <p>
−
                                 Turning the reaction rates in our model into propensities is simple - for example, the propensity for the association reaction
+
                                 The advantage of this method is that we have not had to directly measure expression data for all of our enzymes, which is a difficult process. We conclude that we <em>should</em> obtain enzyme expression of order \(\mu M\) within 350 minutes. However, the scaling factor we introduced in our fitting function is no substitute for a calibration curve to match GFP fluorescence with GFP concentration. For this reason, we conservatively estimate that our proteins are expressed at \(nM\) concentration.
−
                                \(Arab+AraC\rightarrow Arab:AraC\) is \(\alpha_{3}\left[Arab\right]\left[AraC\right]\).
+
 
                             </p>
 
                             </p>
 
                         </div>
 
                         </div>
 
                     </div>
 
                     </div>
 
                 </div>
 
                 </div>
−
                 <div class="section-spacer"></div>
+
                 <div class="image-massive">
−
                 <div class="section" id="diffusion">
+
                    <img src="https://static.igem.org/mediawiki/2015/5/5d/Ox_HenryLarge.jpeg" alt="Our modeller Henry" />
−
                     <div class="slim">
+
                </div>
−
                         <h2>Diffusion</h2>
+
                 <div class="slim">
−
                    </div>
+
                     <div class="section" id="delivery">
−
                    <div id="diffusion-deterministic">
+
                         <h2>Delivery</h2>
−
                         <div class="slim">
+
                        <p>
−
                             <h3>Deterministic</h3>
+
                            With the information about the rates of production and concentrations of our products, we can look at how the products behave once they leave the cell. Our enzymes are first secreted from the cells, and then out of our containment beads to the biofilms they target. We can provide an estimate of the time scale that our project is working on and assess any need for optimisation of enzyme efficiency.
 +
                        </p>
 +
                         <div id="delivery-dispersin">
 +
                             <h3>Dispersin B</h3>
 
                             <p>
 
                             <p>
−
                                 The standard from for the diffusion equation,
+
                                 Dispersin B is one of the anti-<a class="definition" title="biofilm" data-content="A layer of polysaccharides and extracellular DNA that render antibiotics ineffective. They stick to surfaces, which is why you find them on urinary catheters.">biofilm</a> agents we are using in our project and will be the focus of this delivery section. As such we will assume that conclusions reached apply to all of our enzymes.
 
                             </p>
 
                             </p>
−
                                \[\dfrac{d[P]}{dt}=D\dfrac{d^{2}[P]}{dx^{2}}\]
 
 
                             <p>
 
                             <p>
−
                                 where \([P]\) is the concentration of our product, \(t\) is time and \(x\) is distance and \(D\) is the diffusion constant, can be solved analytically in any number of dimensions:
+
                                 A concentration of Dispersin B of 60μg/ml is required to destroy a biofilm that has already formed on a surface [<a href="#references">1</a>]. This equates to a concentration of 1.5μM. This is higher than the steady-state gene expression concentration we can expect from our cells, meaning that our system cannot rely solely on diffusion to transport our enzymes to the biofilm. We will therefore model these diffusion systems assuming that our cells are expressing at a 2μM concentration and later we will look at optimising the gene expression to this level.
−
                            </p>
+
−
                                \[[P]\left(x,t\right)=\dfrac{1}{\left(4\pi Dt\right)^{\frac{N}{2}}}\exp\left(\dfrac{-\sum\limits_{i=1}^N x_{i}^{2}}{4Dt}\right)\]
+
−
                            <p>
+
−
                                where \(x_{i}\) is a generalised position co-ordinate and \(N\) is the dimensionality of the system. This solution only applies to a source of [P] at the position \(x=0\). We can also solve for a grid of 1 or 2D points using finite difference models given in the <a href="#modeling-tutorial">tutorial section</a>.
+
 
                             </p>
 
                             </p>
 
                         </div>
 
                         </div>
−
                    </div>
+
                        <div id="delivery-beads">
−
                    <div id="diffusion-stochastic">
+
                             <h3>Beads</h3>
−
                        <div class="slim">
+
                             <div id="delivery-beads-diffusion">
−
                             <h3>Stochastic</h3>
+
                                 <h4>Diffusion</h4>
−
                             <p>
+
                                <p>
−
                                 We have the option of using three stochastic models for diffusion - random walk, random velocity and compartmental Gillespie.
+
                                    The bead delivery system consists of our cells being contained in <a class="definition" title="alginate" data-content="An inexpensive, non-toxic chemical which can be used to trap cells in place, while letting proteins move freely.">alginate</a> spheres. Water is passed through a container filled with the beads allowing our enzymes to diffuse from the alginate to where they are required. More details about the design of the system can be found <a href="https://2015.igem.org/Team:Oxford/Design#beads">here</a>. Below is a diagram of our cells secreting proteins out of the beads that they are contained in.
−
                            </p>
+
                                </p>
−
                            <h4>Random Walk</h4>
+
                                <div class="image image-full">
−
                            <p>
+
                                    <img  src="https://static.igem.org/mediawiki/2015/6/6a/OxiGEM_Beads_Description_Diagram.png" alt="Diagram of secretion through our AlgiBeads" />
−
                                We initially take a collection of particles and, for each one, compute the number
+
                                    <p>
−
                            </p>
+
                                        Diagram of secretion through our containment beads.
−
                                \[dr=\sqrt{2Ddt}\xi_{r}\]
+
                                    </p>
−
                            <p>
+
                                </div>
−
                                where \(dt\) is the size of our timestep, \(dr\) is the absolute distance moved in the timestep and \(\xi_{r}\) is a normally distributed random number with mean 0 and variance 1. To calculate the distance that particles have actually moved depends on what dimensionality we are operating in. We work in polar or spherical poar co-ordinates in 2D and 3D respectively, and so sample angles \(\phi\) and \(\theta\) and uniformly distributed in the ranges [0,\(2\pi\)] and [0,\(\pi\)].
+
                                <p>
−
                            </p>
+
                                    To determine the <a class="definition" title="convection mass transfer coefficient" data-content="A measure of how quickly a particle can move from one place to another.">convection mass transfer coefficient</a> of Dispersin B from our gel spheres we looked at the diffusion data obtained from <a href="https://2015.igem.org/Team:Oxford/Beads#03-09-2015">this experiment</a> involving the diffusion of <a class="definition" title="crystal violet" data-content="A staining chemical used to mark cell structures and make them more visible under light microscopy. It is also the chemical used in Gram’s staining method for classifying bacteria into those that are Gram positive (the cell walls are stained) and Gram negative (the cell walls are not stained).">crystal violet</a> from our beads. By analysing the system we can produce a theoretical formula for the concentration of crystal violet in the bulk water as a function of time:
−
                            <h4>Random Velocity</h4>
+
                                </p>
−
                            <p>
+
−
                                We initially take a collection of particles at position \(\mathbf{r_{0}}\). We give them random velocities uniformly distributed in the range \(\left[-2\mathbf{v},2\mathbf{v}\right]\) (where \(v\) is determined by thermal fluctuations) and decide that the direction of these velocities should change direction (or not) with some probability after every timestep. We summarise this process in the algorithm below.
+
−
                            </p>
+
−
                            <div class="algorithm">
+
−
                                <h3>Random Velocity Algorithm</h3>
+
−
                                <ol>
+
−
                                    <li>For each particle, generate a random speed \(s\) in a random direction.</li>
+
  
−
                                     <li>Compute the position of the particle after a timestep \(\Delta t\).</li>
+
                                \[c_{f}=\dfrac{c_{bo}}{1+\frac{V_{f}}{V_{b}}}\left(1-\exp\left(\dfrac{-K_{m}A_{b}\left(1+\frac{V_{f}}{V_{b}}\right)t}{V_{f}}\right)\right)\]
 +
                                <table class="table table-striped">
 +
                                     <thead>
 +
                                        <th>Symbol</th>
 +
                                        <th>Definition</th>
 +
                                        <th>Value</th>
 +
                                        <th>Units</th>
 +
                                    </thead>
 +
                                    <tr>
 +
                                        <td>\(A_{b}\)</td>
 +
                                        <td>Total surface area of the beads</td>
 +
                                        <td>\(0.024\)</td>
 +
                                        <td>\(m^{2}\)</td>
 +
                                    </tr>
 +
                                    <tr>
 +
                                        <td>\(V_{b}\)</td>
 +
                                        <td>Total volume of beads</td>
 +
                                        <td>\(1.4\times10^{-5}\)</td>
 +
                                        <td>\(m^{3}\)</td>
 +
                                    </tr>
 +
                                    <tr>
 +
                                        <td>\(c_{bo}\)</td>
 +
                                        <td>Initial dye concentration in beads</td>
 +
                                        <td>\(0.025\)</td>
 +
                                        <td>\(M\)</td>
 +
                                    </tr>
 +
                                    <tr>
 +
                                        <td>\(V_{f}\)</td>
 +
                                        <td>Volume of fluid surrounding the beads</td>
 +
                                        <td>\(V_{f}=V_{fo}-\dfrac{1\times10^{-6}}{10}t\)</td>
 +
                                        <td>\(m^{3}\)</td>
 +
                                    </tr>
 +
                                    <tr>
 +
                                        <td>\(V_{fo}\)</td>
 +
                                        <td>Initial volume of fluid surrounding the beads</td>
 +
                                        <td>\(1\times10^{-4}\)</td>
 +
                                        <td>\(m^{3}\)</td>
 +
                                    </tr>
 +
                                    <tr>
 +
                                        <td>\(t\)</td>
 +
                                        <td>Time</td>
 +
                                        <td>\(-\)</td>
 +
                                        <td>\(min\)</td>
 +
                                    </tr>
 +
                                    <tr>
 +
                                        <td>\(c_{f}\)</td>
 +
                                        <td>Concentration of fluid surrounding beads</td>
 +
                                        <td>\(-\)</td>
 +
                                        <td>\(M\)</td>
 +
                                    </tr>
 +
                                    <tr>
 +
                                        <td>\(K_{m}\)</td>
 +
                                        <td><a class="definition" title="convection mass transfer coefficient" data-content="A measure of how quickly a particle can move from one place to another.">Convection mass transfer coefficient</a></td>
 +
                                        <td>To be fitted</td>
 +
                                        <td>\(mmin^{-1}\)</td>
 +
                                    </tr>
 +
                                </table>
 +
                                <p>
 +
                                    The volume of fluid is also a function of time in order to account for the removal of 1ml of water every 10 minutes. The area and volume of the beads is equal to that of 660 spheres with diameter 3.39mm.
 +
                                </p>
 +
                                <p>
 +
                                    However, the number of beads is an estimate. Because of this, in order to fit the curve to the experimental data we must scale the experimental data by an unknown factor. Therefore we pre-multiply our equation with an arbitrary scaling factor which, along with the convection diffusion coefficient - \(K_{m}\) - is determined by our fitting function.
 +
                                </p>
 +
                                <p>
 +
                                    By fitting the model to our data we returned the value of \(K_{m} = 1.7265\times 10^{-5} mmin^{-1}\).
 +
                                </p>
 +
                                <div class="image image-full">
 +
                                    <img src="https://static.igem.org/mediawiki/2015/b/b7/Ox_fitteddiffusion.png"/>
 +
                                    <p>
 +
                                        We measure the concentration of crystal violet dye as it diffuses out of our containment beads. Errors are given to one standard deviation and data is fitted to a deterministic model to find the convection mass transfer co-efficient. From this we can determine we would require \(100m^{3}\) of beads to reach the required concentration of our own enzymes.
 +
                                    </p>
 +
                                </div>
 +
                                <p>
 +
                                    <a class="definition" title="Dispersin B" data-content="One of the anti-biofilm agents we are using in our project.">Dispersin B</a> is a significantly larger molecule than crystal violet so this diffusion coefficient will not be close to that for Dispersin B. To account for this, we estimated the convective mass transfer coefficient \(K_{m}\) for Dispersin B using that obtained for crystal violet by assuming that \(K_{m}\) is proportional to the diffusion constant in water D.
 +
                                </p>
  
−
                                    <li>Generate a random number \(r\) uniformly distributed between 0 and 1.</li>
+
                                \[\left(K_{m}\right)_{DispersinB} = \dfrac{D_{DispersinB}}{D_{crystal violet}}\left(K_{m}\right)_{crystal violet} \]
  
−
                                     <li>If \(r<\dfrac{s^{2}}{2D}\Delta t\), calculate a new random direction of velocity for the particle.</li>
+
                                <table class="table table-striped">
 +
                                     <thead>
 +
                                        <th>Symbol</th>
 +
                                        <th>Definition</th>
 +
                                        <th>Value</th>
 +
                                        <th>Units</th>
 +
                                    </thead>
 +
                                    <tr>
 +
                                        <td>
 +
                                            \(D_{crystal violet}\)
 +
                                        </td>
 +
                                        <td>
 +
                                            Mass diffusivity of crystal violet in water
 +
                                        </td>
 +
                                        <td>
 +
                                            \(2.87\times10^{9}\)[<a href="#references">2</a>]
 +
                                        </td>
 +
                                        <td>
 +
                                            \(\mu m^{2}s^{-1}\)
 +
                                        </td>
 +
                                    </tr>
 +
                                    <tr>
 +
                                        <td>
 +
                                            \(D_{Dispersin B}\)
 +
                                        </td>
 +
                                        <td>
 +
                                            Mass diffusivity of Dispersin B in water
 +
                                        </td>
 +
                                        <td>
 +
                                            \(100\) [<a href="#references">3</a>]
 +
                                        </td>
 +
                                        <td>
 +
                                            \(\mu m^{2}s^{-1}\)
 +
                                        </td>
 +
                                    </tr>
 +
                                </table>
 +
                                <p>
 +
                                    By substituting in these values we arrive at \(\left(K_{m}\right)_{DispersinB} = 6.03\times10^{-13} mmin^{-1}\)
 +
                                </p>
 +
                            </div>
 +
                            <div id="delivery-beads-mass-exchange">
 +
                                <h4>Mass Exchange</h4>
 +
                                <p>
 +
                                    This result allows us to theorise a mass exchange system. As a first estimate we will assume that the flow through the beads is sufficiently slow to use the convection diffusion coefficient found above. It is also assumed that the gene expression happens on a faster time scale than the diffusion from the beads to the water, enabling us to assume the concentration of enzyme in the beads remains constant. This is supported by our gene expression models. We can now visualize how the concentrations of the fluid will vary with distance along the mass exchanger:
 +
                                </p>
 +
                                <div class="image image-full">
 +
                                    <img src="https://static.igem.org/mediawiki/2015/e/e0/Ox_mass_exchange_daig.jpg"/>
 +
                                    <p>
 +
                                        Visualisation of the concentrations of the fluid and the beads along our mass exchanger
 +
                                    </p>
 +
                                </div>
 +
                                <p>
 +
                                    The overall system can now be described with the equation:
 +
                                </p>
  
−
                                     <li>Repeat steps 2-5 for each particle until a certain time or computational budget is exhausted.</li>
+
                                \[J = K_mA\dfrac{c_{fo}-c_{fi}}{\ln\left(\dfrac{c_{B}-c_{fi}}{c_{B}-c_{fo}}\right)}\]
−
                                 </ol>
+
                                <p>
 +
                                     Therefore
 +
                                </p>
 +
 
 +
                                \[A = J\dfrac{\ln\left(\dfrac{c_{B}-c_{fi}}{c_{B}-c_{fo}}\right)}{K_{m}\left(c_{fo}-c_{fi}\right)}\]
 +
                                <p>
 +
                                    Where \(J=Q\left(c_{fo}-c_{fi}\right)\) and \(Q\) is the volume flow rate of water. We have chosen a flow rate range of 10-100ml/min as this is accepted as a safe artificial bladder fill rate [<a href="#references">4</a>]. This range results in the following number of beads required to reach the desired concentration:
 +
                                </p>
 +
                                <div class="image image-full">
 +
                                    <img src="https://static.igem.org/mediawiki/2015/1/18/Ox_numbervsflow.png"/>
 +
                                    <p>
 +
                                        Relationship between the number of bacteria-containment beads required to reach a particular flow rate of our enzymes. These are the flow rates we require for practical use.
 +
                                    </p>
 +
                                </div>
 +
                                <p>
 +
                                    Therefore a volume of between \(20.3-203m^3\) of beads is required, assuming a <a class="definition" title="packing efficiency" data-content="Our containment-beads are spherical. Therefore, in a collection of beads, some space will be taken up by the gaps between the beads.">packing efficiency</a> of 64% [<a href="#references">9</a>].
 +
                                </p>
 +
                                <p>
 +
                                    This estimation relied upon the flow of fluid around the beads being sufficiently slow such that it may be approximated as stationary, so that mass transfer occurs as natural convection. However, because of the (likely) large volume of beads compared to the cross-sectional area of the catheter, this flow of fluid may have a non-negligable velocity.
 +
                                 </p>
 
                             </div>
 
                             </div>
−
                            </p>
 
−
                            <h4>Compartmental diffusion; Gillespie</h4>
 
−
                            <p>
 
−
                                We can divide our region into a system of compartments
 
−
                            </p>
 
−
                                \[A_{1} \mathrel{\mathop{\rightleftharpoons}^{\mathrm{k}}_{\mathrm{k}}} A_{2} \mathrel{\mathop{\rightleftharpoons}^{\mathrm{k}}_{\mathrm{k}}} ... \mathrel{\mathop{\rightleftharpoons}^{\mathrm{k}}_{\mathrm{k}}} A_{N-1} \mathrel{\mathop{\rightleftharpoons}^{\mathrm{k}}_{\mathrm{k}}} A_{N}\]
 
−
                            <p>
 
−
                                such that the only difference in their propensities is related to the number of particles contained in each compartment. This will produce a correct model of diffusion if the proportionality constant \(k\) is chosen as
 
−
                            </p>
 
−
                                \[k=\dfrac{D}{(dx)^{2}}\]
 
−
                            <p>
 
−
                                where \(dx\) is the physical size of each compartment [<a href="#Ref2">2</a>]. We now run the Gillespie algorithm but with \(2(N-1)\) reactions to consider.
 
−
                            </p>
 
−
                            <h4>Comparison of methods</h4>
 
−
                            <p>
 
−
                                All of the methods provide similar levels of accuracy (depending on the size of the compartments used in the Gillespie model) and so we choose based on the simplicity of implementation. The random walk method is generally the fastest and easiest to code, so we choose this one.
 
−
                            </p>
 
−
                            <h4>The challenge</h4>
 
−
                            <p>
 
−
                                Since we are modelling diffusion in the catheter, the topology of the system is a ring, neglecting edge effects of the cylinder our bacteria will be placed in. To do this we use a random walk model in two dimensions. We confine our particles to a ring of a fixed radius and initially coat the ring with particles.
 
−
                            </p>
 
 
                         </div>
 
                         </div>
 
                     </div>
 
                     </div>
−
                </div>
+
                    <div class="section-spacer">
−
                <div class="section-spacer"></div>
+
−
                <div class="section" id="modeling-tutorial">
+
−
                    <div class="slim">
+
−
                        <h2>Modeling Tutorial</h2>
+
 
                     </div>
 
                     </div>
−
                     <div id="modeling-tutorial-deterministic">
+
                     <div class="section" id="conclusion">
−
                        <div class="slim">
+
                        <h2>Conclusion</h2>
−
                            <h3>Deterministic</h3>
+
                        <p>
−
                            <p>
+
                             Using gene expression and diffusion models, we estimated that we would need around \(100m^3\) of beads to deliver enough enzymes to clear a urinary-infection-associated biofilm. Our treatment will only be more effective than antibiotics if we make our enzymes many orders of magnitude more efficient. This led us to consider an alternative design.
−
                                Our goal is to model a system of reactions such as the reversible reaction
+
                         </p>
−
                             </p>
+
−
                                \[A \mathrel{\mathop{\rightleftharpoons}^{\mathrm{k_{+}}}_{\mathrm{k_{-}}}} B\]
+
−
                            <p>
+
−
                                Where \(A\) and \(B\) are our two species and \(k_{+}\) and \(k_{-}\) will be determined shortly. To do this, we make the intuitive assumption that the rate at which something reacts is proportional to how much of it we've got. That is, the rate of change of \(B\) with time is
+
−
                            </p>
+
−
                                \[\dfrac{dB}{dt}=k_{+}A\]
+
−
                            <p>
+
−
                                where \(k_{+}\) is simply a proportionality constant. This assumption is called the law of mass action. This of course is not the complete picture, as we forgot to include the fact that \(B\) can react back to form \(A\). So the complete system of equations to solve is
+
−
                            </p>
+
−
                            \[\dfrac{dA}{dt}=k_{-}B-k_{+}A\]
+
−
                            \[\dfrac{dB}{dt}=k_{+}A-k_{-}B\]
+
−
                            <p>
+
−
                                At this point, you can type these into your favourite computing program and solve them.
+
−
                            </p>
+
−
                        </div>
+
−
                    </div>
+
−
                    <div id="modeling-tutorial-stochastic">
+
−
                        <div class="slim">
+
−
                            <h3>Stochastic</h3>
+
−
                            <p>
+
−
                                The simplest way to add an element of 'randomness' into an equation is to add a random term into the differential equation. This turns it into a stochastic differential equation, such as:
+
−
                            </p>
+
−
                                \[\dfrac{dx}{dt}=c+r\]
+
−
                            <div class="image image-right">
+
−
                                <img src="https://static.igem.org/mediawiki/2015/1/1f/OxiGEM_Stochastic_Tutorial_01.jpg"/>
+
−
                            </div>
+
−
                            <p>
+
−
                                which determines the evolution of a particle travelling at a speed \(c\), with a random number \(r\) added into the mix just for fun. We choose \(r\) to be normally distributed with a mean value of 0 and a variance of, say, 1. We evaluate the position of the particle at multiple timesteps (with some small time \(dt\) separating them) and plot our solution in the Figure below.
+
−
                            </p>
+
−
                            <p>
+
−
                                We can do better than this when dealing with systems such as our example in the previous section. We use the Gillespie algorithm to find:
+
−
                            </p>
+
−
                            <ol>
+
−
                                <li>The <strong>time</strong> taken between reactions</li>
+
−
                                <li>The <strong>reaction</strong> that took place in that time.</li>
+
−
                            </ol>
+
−
                            <p>
+
−
                                as neither of these could be determined by our initial method [<a href="#Ref2">2</a>, <a href="#Ref4">4</a>].
+
−
                            </p>
+
−
                            <div class="image image-left">
+
−
                                <img src="https://static.igem.org/mediawiki/2015/7/72/OxiGEM_Stochastic_Tutorial_02.png"/>
+
−
                            </div>
+
−
                            <p>
+
−
                                <strong>1)</strong> To find the time taken, we make an assumption about the form of the probability distribution that the time between reactions holds. It makes sense to think of it as a falling exponential, such that reactions more often happen quickly rather than slowly, and that the more probable a reaction is, the shorter the time between subsequent reactions. This naturally leads to the form
+
−
                            </p>
+
−
                                \[P(t)=e^{-\alpha_{0}t}\]
+
−
                            <p>
+
−
                                where \(P(t)\) is the probability of any reaction occuring in time \(t\), and \(\alpha_{0}\) tells us something about the likelihood of the reaction. We plot this in the Figure. In fact we call \(\alpha_{0}\) the sum of the propensities \(\alpha_{i}\) (where the dummy variable \(i\) runs from 1 to the number of possible reactions, \(N\), that could occur). In the case of our example system, the propensity for the reaction \(A\) to \(B\) would be \(k_{+}A\) and the propensity for the reaction \(B\) to \(A\) is \(k_{-}B\) so in this case \(\alpha_{0}=k_{+}A+k_{-}B\).
+
−
                            </p>
+
−
                            <p>
+
−
                                To decide the time it took for a reaction to occur, we sample a random number \(r_{1}\) in the range [0,1] and stick that number on the y axis of the graph. We read off from our graph a corresponding time and we're done - we now know when the next reaction took place.
+
−
                            </p>
+
−
                            <p>
+
−
                                <strong>2)</strong> To decide which reaction took place, we need the probabilities of each reaction occuring being proportional to the ratio of their propensities. We split up the domain of 0 to 1 into the normalised ratios of the propensities of our problem, recalling the propensities we calculated. We could have, in this case, selected the backwards reaction B to A. We can generalise this to the form given in the Gillespie Algorithm above.
+
−
                            </p>
+
−
                            <!--CONTENT MISSING-->
+
−
                            <div class="image image-full">
+
−
                                <img src="https://static.igem.org/mediawiki/2015/d/da/OxiGEM_Stochastic_Tutorial_03.jpg"/>
+
−
                            </div>
+
−
                         </div>
+
−
                    </div>
+
−
                    <div id="modeling-tutorial-finite-difference">
+
−
                        <div class="slim">
+
−
                            <h3>Finite Difference Models</h3>
+
−
                            <p>
+
−
                                For finite difference models, we set up a grid of positions that we want to solve our equations at, separated by some fixed difference \(dx\) and at time steps separated by time \(dt\). Subscripts of our concentration P represent the spacial co-ordinates of the grid-point, while superscript represents the location in time. We define a constant \(r=\dfrac{dt}{(dx)^{2}}\) which must be less than 0.5 for the scheme to be convergent.
+
−
                            </p>
+
−
                            <p>
+
−
                                The 1D formula to iterate over for the finite difference model of diffusion is
+
−
                            </p>
+
−
                                \[P_{j}^{n+1}=(1-2r)P_{j}^{n}+r(P_{j-1}^{n}+P_{j+1}^{n})\]
+
−
                            <p>
+
−
                                The 2D formula to iterate over is
+
−
                            </p>
+
−
                                \[P_{i,j}^{n+1}=(1-4r)P_{i,j}^{n}+r(P_{i,j-1}^{n}+P_{i,j+1}^{n}+P_{i-1,j}^{n}+P_{i+1,j}^{n})\]
+
−
                            <p>
+
−
                                Similarly, for 3D:
+
−
                            </p>
+
−
                                \[P_{i,j,k}^{n+1}=(1-6r)P_{i,j,k}^{n}+r(P_{i,j-1,k}^{n}+P_{i,j+1,k}^{n}+P_{i-1,j,k}^{n}+P_{i+1,j,k}^{n}+P_{i,j,k-1}^{n}+P_{i,j,k+1}^{n})\]
+
−
                            <p>
+
−
                                we (arbitrarily) choose our boundary conditions to be continuous at the boundary.
+
−
                            </p>
+
−
                        </div>
+
 
                     </div>
 
                     </div>
 
                 </div>
 
                 </div>
Line 371: Line 454:
 
                     <h2>References</h2>
 
                     <h2>References</h2>
 
                     <ol class="references">
 
                     <ol class="references">
−
                         <li>Synthetic Biology - A Primer, Paul S. Freemont, Richard I Kenne et al., 2012 Imperial College Press<a name="Ref1"></a></li>
+
                         <li>Jeffrey B. Kaplan; Dispersin B polypeptides and uses thereof. Patent PI 8580551, Nov 12, 2013</li>
  
−
                         <li>R. Erban, S. J. Chapman, and P. Maini. A practical guide to stochastic simulations of reaction- diffusion processes, 2007. Available <a href="http://arxiv.org/abs/0704.1908">here</a><a name="Ref2"></a></li>
+
                         <li><a href="http://physicalpharmacy2013.blogspot.co.uk/2013/05/practical-4.html">http://physicalpharmacy2013.blogspot.co.uk/2013/05/practical-4.html</a></li>
  
−
                         <li>Salto R, Delgado A, Michán C, Marqués S, Ramos JL. Modulation of the function of the signal receptor domain of XylR, a member of a family of prokaryotic enhancer-like positive regulators. J Bacteriol. 1998 Feb180(3):600-4. p.601 right column<a name="Ref3"></a></li>
+
                         <li>"Physical Biology of the Cell", Rob Phillips, Jane Kondev and Julie Theriot (2009). Page 110</li>
  
−
                         <li>Brian Ingalls,Mathematical Modeling in Systems Biology: An Introduction, 2012 University of Waterloo<a name="Ref4"></a></li>
+
                         <li>Kim S-Y, Ko SH, Shin MJ, et al. Phasic Changes in Bladder Compliance During Filling Cystometry of the Neurogenic Bladder. Annals of Rehabilitation Medicine. 2014;38(3):342-346. doi:10.5535/arm.2014.38.3.342.</li>
  
−
                         <li>Liang ST, Ehrenberg M, Dennis P, Bremer H. Decay of rplN and lacZ mRNA in Escherichia coli. J Mol Biol. 1999 May 14 288(4):521-38. p.524 right column bottom paragraph<a name="Ref5"></a></li>
+
                         <li>Liang ST, Ehrenberg M, Dennis P, Bremer H. Decay of rplN and lacZ mRNA in Escherichia coli. J Mol Biol. 1999 May 14 288(4):521-38. p.524 right column bottom paragraph</li>
  
−
                         <li>Proshkin S, Rahmouni AR, Mironov A, Nudler E. Cooperation between translating ribosomes and RNA polymerase in transcription elongation. Science. 2010 Apr 23 328(5977):504-8. p.505 table 1<a name="Ref6"></a></li>
+
                         <li>Proshkin S, Rahmouni AR, Mironov A, Nudler E. Cooperation between translating ribosomes and RNA polymerase in transcription elongation. Science. 2010 Apr 23 328(5977):504-8. p.505 table 1</li>
  
−
                         <li>Nelson HC, Sauer RT. Lambda repressor mutations that increase the affinity and specificity of operator binding.Cell. 1985 Sep42(2):549-58. p.552 table 3<a name="Ref7"></a></li>
+
                         <li>Sourjik V, Berg HC. Functional interactions between receptors in bacterial chemotaxis. Nature. 2004 Mar 25 428(6981):437-41.p.439 left column top paragraph</li>
  
−
                         <li>Sourjik V, Berg HC. Functional interactions between receptors in bacterial chemotaxis. Nature. 2004 Mar 25 428(6981):437-41.p.439 left column top paragraph<a name="Ref8"></a></li>
+
                         <li>Salto R, Delgado A, Michán C, Marqués S, Ramos JL. Modulation of the function of the signal receptor domain of XylR, a member of a family of prokaryotic enhancer-like positive regulators. J Bacteriol. 1998 Feb180(3):600-4. p.601 right column</li>
 +
 
 +
                        <li>S. Torquato, T. M. Truskett, and P. G. Debenedetti Is Random Close Packing of Spheres Well Defined? 2000 Phys. Rev. Lett. 84, 2064</li>
 +
 
 +
                        <li>Selinger DW, Saxena RM, Cheung KJ, Church GM, Rosenow C. Global RNA half-life analysis in Escherichia coli reveals positional patterns of transcript degradation. Genome Res. 2003 Feb13(2):216-23. p.217 left column 2nd paragraph</li>
 +
 
 +
                        <li>Andersen JB, Sternberg C, Poulsen LK, Bjorn SP, Givskov M, Molin S. New Unstable Variants of Green Fluorescent Protein for Studies of Transient Gene Expression in Bacteria. Appl Environ Microbiol. 1998 Jun64(6):2240-6.</li>
 
                     </ol>
 
                     </ol>
 
                 </div>
 
                 </div>
Line 395: Line 484:
 
                     </li>
 
                     </li>
 
                     <li>
 
                     <li>
−
                         <a href="#gene-expression">Gene Expression</a>
+
                         <a href="#characterising-our-cells">Gene expression rates</a>
−
                        <ul class="nav nav-stacked">
+
−
                            <li><a href="#gene-expression-deterministic">Deterministic</a></li>
+
−
                            <li><a href="#gene-expression-stochastic">Stochastic</a></li>
+
−
                        </ul>
+
−
                    </li>
+
−
                    <li>
+
−
                        <a href="#diffusion">Diffusion</a>
+
 
                         <ul class="nav nav-stacked">
 
                         <ul class="nav nav-stacked">
−
                             <li><a href="#diffusion-deterministic">Deterministic</a></li>
+
                             <li><a href="#characterising-our-cells-arab">Arabinose-induced promoter</a></li>
−
                            <li><a href="#diffusion-stochastic">Stochastic</a></li>
+
 
                         </ul>
 
                         </ul>
 
                     </li>
 
                     </li>
 
                     <li>
 
                     <li>
−
                         <a href="#modeling-tutorial">Modeling Tutorial</a>
+
                         <a href="#delivery">Delivery</a>
 
                         <ul class="nav nav-stacked">
 
                         <ul class="nav nav-stacked">
−
                             <li><a href="#modeling-tutorial-deterministic">Deterministic</a></li>
+
                             <li><a href="#delivery-dispersin">Dispersin B</a></li>
−
                             <li><a href="#modeling-tutorial-stochastic">Stochastic</a></li>
+
                             <li>
−
                            <li><a href="#modeling-tutorial-finite-difference">Finite Difference Models</a></li>
+
                                <a href="#delivery-beads">Beads</a>
 +
                                <ul class="nav nav-stacked">
 +
                                    <li><a href="#delivery-beads-diffusion">Diffusion</a></li>
 +
                                    <li><a href="#delivery-beads-mass-exchange">Mass Exchanger</a></li>
 +
                                </ul>
 +
                            </li>
 
                         </ul>
 
                         </ul>
 +
                        <a href="#conclusion">Conclusion</a>
 
                     </li>
 
                     </li>
 
                     <li>
 
                     <li>

Latest revision as of 03:00, 19 September 2015

Modelling

Introduction

Mathematical modelling plays a crucial role in Synthetic Biology by acting as a link between the conception and the physical realisation of a biological circuit. Our modelling team has evaluated the effectiveness of initial designs, and has provided insight into how the system can (or must) be improved.

Our team experimentally validated that Escherichia coli can secrete enzymes which break down the biofilms associated with urinary infections. However, it is difficult to directly measure whether our enzymes are produced in a sufficient quantity to be a more effective treatment than antibiotics. We measured gene expression and diffusion of widely-used chemicals, and then used our model to estimate the number of E. coli cells that would make our project a more effective treatment than antibiotics. We expect to have to improve our system to make it realistic.

To help readers of all kinds and specialisations understand this page we have produced guides for all the modelling techniques used in this section. They are available in our Modelling Tutorial page and will be linked to when appropriate.

Gene expression rates

In this section we look at our cells in isolation in order to assess their functionality and answer important questions such as “how long does it take to produce a certain concentration of product?” The end result - the final concentration of useful enzyme that is produced in the cell - is required for our diffusion model.

Arabinose-induced expression

We have decided to use an arabinose-induced promoter for the expression of a number of our proteins. This promoter can be modelled as the following chemical system:

\[(Arab:AraC)\overset{K}{\rightarrow}mRNA\overset{\alpha}{\rightarrow}P\] \[mRNA\overset{\gamma_{1}}{\rightarrow}\phi\quad P\overset{\gamma_{2}}{\rightarrow}\phi\]

Our promoter, pBAD, binds to AraC and this represses transcription of mRNA. Arabinose will bind to AraC and form the Arab:AraC compound, allowing transcription to occur.

For this system we will assume that AraC is always in large concentration and that its binding to arabinose happens on a faster time scale to transcription. Therefore, we do not need to consider the individual concentrations of arabinose and AraC, instead we just need to include the concentration of the complex (Arab:AraC). The rate \(K\) is not just a simple constant and is given as the Hill function in the equations below.

Using Michaelis-Mentin kinetics, we arrive at the equations:

\[\dfrac{d[mRNA]}{dt}=K_{max}\dfrac{[Arab:AraC]^{n}}{K_{half}^{n}+[Arab:AraC]^{n}}-\gamma_{1}[mRNA]\] \[\dfrac{d\left[P\right]}{dt}=\alpha\left[mRNA\right]-\gamma_{2}\left[P\right]\]

Where \([Arab]\), \([AraC]\), \([Arab:AraC]\), \([mRNA]\) and \([P]\) represent the concentrations of arabinose, AraC, Arab:AraC, mRNA and our product protein respectively. We define the remaining symbols in the table below.

Symbol Definition Initial Value/Literature Value Fitted
\(\alpha\) Translation rate \(15ntd\: s^{-1}\)/length of sequence [6] \(6.6ntd\: s^{-1}\)/length of sequence
\(\gamma_{1}\) Combined degradation and dilution rate of mRNA \(2.2\times10^{-3}s^{-1}\) [5, 10] \(1.1\times10^{-2}s^{-1}\)
\(\gamma_{2}\) Combined degradation and dilution rate of GFP \(5.2\times10^{-4}s^{-1}\) [5, 11] \(1.1\times10^{-2}s^{-1}\)
\(K_{max}\) Maximal transcription rate \(50ntd\: s^{-1}\)/length of sequence [6] \(47ntd\: s^{-1}\)/length of sequence
\(K_{half}\) Half-maximal transcription rate \(160\mu M\) [7] \(100\mu M\)
\(n\) Hill coefficient \(2.65\) [8] \(2.73\)

This table contains literature values for the parameters, found from a number of sources. We then measured GFP expression in E. coli to extract experimental values. Here is a plot showing the fit of this model to our experimental data.

Fitting our gene expression data to the theoretical model

Results showing GFP concentration as a function of time, matched to our deterministic model. Errors are given to one standard deviation and an arbitrary scaling factor is included as a fitted parameter.

We can now calculate the limiting concentrations that our products will be expressed. The dilution rate and Hill coefficient of our cells is the same for GFP and our proteins, but the transcription and translation rates are dependent on the sequence length of the protein. Here is a table showing the relevant proteins and sequence lengths:

Product Sequence Length (/bp)
pBAD HisB DNase DsbA 621
pBAD HisB MccS 414
pBAD HisB Art-175 DsbA 987
pBAD HisB Art-175 YebF 1284
pBAD HisB Art-E 632
pBAD HisB Art-175 Fla 1095
pBAD HisB Art-175 936
pBAD HisB DNase 570

We now can run our model of the system by solving the set of equations using the MATLAB ordinary differential equation solver ode15s. Below is a plot of the concentration of product against time for each protein expressed with this inducer-promoter pair where the expression is induced by a step function:

Model data for each of the enzymes we plan to release, using the parameters we found from our experimental data. We found our limiting concentrations were of order \(\mu M\).

The advantage of this method is that we have not had to directly measure expression data for all of our enzymes, which is a difficult process. We conclude that we should obtain enzyme expression of order \(\mu M\) within 350 minutes. However, the scaling factor we introduced in our fitting function is no substitute for a calibration curve to match GFP fluorescence with GFP concentration. For this reason, we conservatively estimate that our proteins are expressed at \(nM\) concentration.

Our modeller Henry

Delivery

With the information about the rates of production and concentrations of our products, we can look at how the products behave once they leave the cell. Our enzymes are first secreted from the cells, and then out of our containment beads to the biofilms they target. We can provide an estimate of the time scale that our project is working on and assess any need for optimisation of enzyme efficiency.

Dispersin B

Dispersin B is one of the anti-biofilm agents we are using in our project and will be the focus of this delivery section. As such we will assume that conclusions reached apply to all of our enzymes.

A concentration of Dispersin B of 60μg/ml is required to destroy a biofilm that has already formed on a surface [1]. This equates to a concentration of 1.5μM. This is higher than the steady-state gene expression concentration we can expect from our cells, meaning that our system cannot rely solely on diffusion to transport our enzymes to the biofilm. We will therefore model these diffusion systems assuming that our cells are expressing at a 2μM concentration and later we will look at optimising the gene expression to this level.

Beads

Diffusion

The bead delivery system consists of our cells being contained in alginate spheres. Water is passed through a container filled with the beads allowing our enzymes to diffuse from the alginate to where they are required. More details about the design of the system can be found here. Below is a diagram of our cells secreting proteins out of the beads that they are contained in.

Diagram of secretion through our AlgiBeads

Diagram of secretion through our containment beads.

To determine the convection mass transfer coefficient of Dispersin B from our gel spheres we looked at the diffusion data obtained from this experiment involving the diffusion of crystal violet from our beads. By analysing the system we can produce a theoretical formula for the concentration of crystal violet in the bulk water as a function of time:

\[c_{f}=\dfrac{c_{bo}}{1+\frac{V_{f}}{V_{b}}}\left(1-\exp\left(\dfrac{-K_{m}A_{b}\left(1+\frac{V_{f}}{V_{b}}\right)t}{V_{f}}\right)\right)\]
Symbol Definition Value Units
\(A_{b}\) Total surface area of the beads \(0.024\) \(m^{2}\)
\(V_{b}\) Total volume of beads \(1.4\times10^{-5}\) \(m^{3}\)
\(c_{bo}\) Initial dye concentration in beads \(0.025\) \(M\)
\(V_{f}\) Volume of fluid surrounding the beads \(V_{f}=V_{fo}-\dfrac{1\times10^{-6}}{10}t\) \(m^{3}\)
\(V_{fo}\) Initial volume of fluid surrounding the beads \(1\times10^{-4}\) \(m^{3}\)
\(t\) Time \(-\) \(min\)
\(c_{f}\) Concentration of fluid surrounding beads \(-\) \(M\)
\(K_{m}\) Convection mass transfer coefficient To be fitted \(mmin^{-1}\)

The volume of fluid is also a function of time in order to account for the removal of 1ml of water every 10 minutes. The area and volume of the beads is equal to that of 660 spheres with diameter 3.39mm.

However, the number of beads is an estimate. Because of this, in order to fit the curve to the experimental data we must scale the experimental data by an unknown factor. Therefore we pre-multiply our equation with an arbitrary scaling factor which, along with the convection diffusion coefficient - \(K_{m}\) - is determined by our fitting function.

By fitting the model to our data we returned the value of \(K_{m} = 1.7265\times 10^{-5} mmin^{-1}\).

We measure the concentration of crystal violet dye as it diffuses out of our containment beads. Errors are given to one standard deviation and data is fitted to a deterministic model to find the convection mass transfer co-efficient. From this we can determine we would require \(100m^{3}\) of beads to reach the required concentration of our own enzymes.

Dispersin B is a significantly larger molecule than crystal violet so this diffusion coefficient will not be close to that for Dispersin B. To account for this, we estimated the convective mass transfer coefficient \(K_{m}\) for Dispersin B using that obtained for crystal violet by assuming that \(K_{m}\) is proportional to the diffusion constant in water D.

\[\left(K_{m}\right)_{DispersinB} = \dfrac{D_{DispersinB}}{D_{crystal violet}}\left(K_{m}\right)_{crystal violet} \]
Symbol Definition Value Units
\(D_{crystal violet}\) Mass diffusivity of crystal violet in water \(2.87\times10^{9}\)[2] \(\mu m^{2}s^{-1}\)
\(D_{Dispersin B}\) Mass diffusivity of Dispersin B in water \(100\) [3] \(\mu m^{2}s^{-1}\)

By substituting in these values we arrive at \(\left(K_{m}\right)_{DispersinB} = 6.03\times10^{-13} mmin^{-1}\)

Mass Exchange

This result allows us to theorise a mass exchange system. As a first estimate we will assume that the flow through the beads is sufficiently slow to use the convection diffusion coefficient found above. It is also assumed that the gene expression happens on a faster time scale than the diffusion from the beads to the water, enabling us to assume the concentration of enzyme in the beads remains constant. This is supported by our gene expression models. We can now visualize how the concentrations of the fluid will vary with distance along the mass exchanger:

Visualisation of the concentrations of the fluid and the beads along our mass exchanger

The overall system can now be described with the equation:

\[J = K_mA\dfrac{c_{fo}-c_{fi}}{\ln\left(\dfrac{c_{B}-c_{fi}}{c_{B}-c_{fo}}\right)}\]

Therefore

\[A = J\dfrac{\ln\left(\dfrac{c_{B}-c_{fi}}{c_{B}-c_{fo}}\right)}{K_{m}\left(c_{fo}-c_{fi}\right)}\]

Where \(J=Q\left(c_{fo}-c_{fi}\right)\) and \(Q\) is the volume flow rate of water. We have chosen a flow rate range of 10-100ml/min as this is accepted as a safe artificial bladder fill rate [4]. This range results in the following number of beads required to reach the desired concentration:

Relationship between the number of bacteria-containment beads required to reach a particular flow rate of our enzymes. These are the flow rates we require for practical use.

Therefore a volume of between \(20.3-203m^3\) of beads is required, assuming a packing efficiency of 64% [9].

This estimation relied upon the flow of fluid around the beads being sufficiently slow such that it may be approximated as stationary, so that mass transfer occurs as natural convection. However, because of the (likely) large volume of beads compared to the cross-sectional area of the catheter, this flow of fluid may have a non-negligable velocity.

Conclusion

Using gene expression and diffusion models, we estimated that we would need around \(100m^3\) of beads to deliver enough enzymes to clear a urinary-infection-associated biofilm. Our treatment will only be more effective than antibiotics if we make our enzymes many orders of magnitude more efficient. This led us to consider an alternative design.

References

  1. Jeffrey B. Kaplan; Dispersin B polypeptides and uses thereof. Patent PI 8580551, Nov 12, 2013
  2. http://physicalpharmacy2013.blogspot.co.uk/2013/05/practical-4.html
  3. "Physical Biology of the Cell", Rob Phillips, Jane Kondev and Julie Theriot (2009). Page 110
  4. Kim S-Y, Ko SH, Shin MJ, et al. Phasic Changes in Bladder Compliance During Filling Cystometry of the Neurogenic Bladder. Annals of Rehabilitation Medicine. 2014;38(3):342-346. doi:10.5535/arm.2014.38.3.342.
  5. Liang ST, Ehrenberg M, Dennis P, Bremer H. Decay of rplN and lacZ mRNA in Escherichia coli. J Mol Biol. 1999 May 14 288(4):521-38. p.524 right column bottom paragraph
  6. Proshkin S, Rahmouni AR, Mironov A, Nudler E. Cooperation between translating ribosomes and RNA polymerase in transcription elongation. Science. 2010 Apr 23 328(5977):504-8. p.505 table 1
  7. Sourjik V, Berg HC. Functional interactions between receptors in bacterial chemotaxis. Nature. 2004 Mar 25 428(6981):437-41.p.439 left column top paragraph
  8. Salto R, Delgado A, Michán C, Marqués S, Ramos JL. Modulation of the function of the signal receptor domain of XylR, a member of a family of prokaryotic enhancer-like positive regulators. J Bacteriol. 1998 Feb180(3):600-4. p.601 right column
  9. S. Torquato, T. M. Truskett, and P. G. Debenedetti Is Random Close Packing of Spheres Well Defined? 2000 Phys. Rev. Lett. 84, 2064
  10. Selinger DW, Saxena RM, Cheung KJ, Church GM, Rosenow C. Global RNA half-life analysis in Escherichia coli reveals positional patterns of transcript degradation. Genome Res. 2003 Feb13(2):216-23. p.217 left column 2nd paragraph
  11. Andersen JB, Sternberg C, Poulsen LK, Bjorn SP, Givskov M, Molin S. New Unstable Variants of Green Fluorescent Protein for Studies of Transient Gene Expression in Bacteria. Appl Environ Microbiol. 1998 Jun64(6):2240-6.