Difference between revisions of "Team:Oxford/Modeling/Tutorial"

Line 2: Line 2:
  
 
<html>
 
<html>
−
     <div class="container-fluid page-heading" style="background-image: url(https://static.igem.org/mediawiki/2015/6/6d/Ox_modelling_will_poster.png)">
+
     <div class="container-fluid page-heading">
−
         <h3>Modeling Tutorials</h3>
+
         <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>
−
                             On this page we aim to explain the concepts and methodology behind the <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.">model</a>ling we have done in our project. We hope this will make our modeling section more accessible for all members of the iGEM community.
+
                             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>
 +
                            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-networks">
+
                         <h2>Gene expression rates</h2>
−
                    <div class="slim">
+
                         <p>
−
                         <h2>Gene Expression Networks</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-networks-mass-action">
+
                        </p>
−
                            <h3>Mass Action</h3>
+
                         <div id="characterising-our-cells-arab">
−
                            <p>
+
                             <h3>Arabinose-induced expression</h3>
−
                                Our goal is to model a system of reactions such as the reversible reaction
+
−
                             </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 id="gene-expression-networks-michaelis">
+
−
                             <h3>Michaelis-Menten Kinetics</h3>
+
 
                             <p>
 
                             <p>
−
                                 It is possible to formulate a rate law that describes <a class="definition" title="enzyme" data-content="A molecule which speeds up a chemical reaction - a biological catalyst. The reaction does not involve these molecules.">enzyme</a>-catalysed reactions, this is known as Michaelis-Menten kinetics. In order to derive this rate law we will look at a generalised system of chemical events detailed below:
+
                                 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>
  
−
                             \[S + E\mathop{\rightleftharpoons}C_{1}\mathop{\rightleftharpoons}C_{2}\mathop{\rightleftharpoons}P + E\]
+
                             \[(Arab:AraC)\overset{K}{\rightarrow}mRNA\overset{\alpha}{\rightarrow}P\]
−
 
+
                             \[mRNA\overset{\gamma_{1}}{\rightarrow}\phi\quad P\overset{\gamma_{2}}{\rightarrow}\phi\]
−
                             <p><em>
+
−
                                In this system \(S\) represents our substrate, \(E\) our enzyme, \(C_{1}\) our enzyme-substrate complex, \(C_{2}\) our enzyme-product complex, and \(P\) our product.
+
−
                            </em></p>
+
 
                             <p>
 
                             <p>
−
                                 Initially we will make two simplifications. Firstly, we will combine \(C_{1}\) and C2 into a single complex, \(C\), as we will assume that the time-scale of the conversion \(C_{1}\mathop{\rightleftharpoons}C_{2}\) is much faster than that of the association and dissociation events. Secondly, we will assume that the product never binds with the free enzyme. These two assumptions lead to the simplified network:
+
                                 Our promoter, pBAD, is known as a double repressor. <a class="definition" title="AraC" data-content="A compound which stops our E. coli from producing enzymes">AraC</a> binds to pBAD which 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>. By introducing arabinose into the system, AraC will bind to arabinose and form the Arab:AraC compound. Transcription can then occur.
 
                             </p>
 
                             </p>
−
 
−
                            \[S + E\mathrel{\mathop{\rightleftharpoons}^{k_{1}}_{k_{-1}}}C\mathop{\rightarrow}^{k_{2}}P + E\]
 
−
 
 
                             <p>
 
                             <p>
−
                                 Using the laws of mass action (detailed above) we arrive at the following differential equation model:
+
                                 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>
−
 
−
                            \[\dfrac{d}{dt}s(t)=-k_{1}s(t)e(t) + k_{-1}c(t)\]
 
−
 
−
                            \[\dfrac{d}{dt}e(t)=k_{-1}c(t) - k_{1}s(t)e(t) + k_{2}c(t)\]
 
−
 
−
                            \[\dfrac{d}{dt}c(t)=-k_{-1}c(t) + k_{1}s(t)e(t) - k_{2}c(t)\]
 
−
 
−
                            \[\dfrac{d}{dt}p(t) = k_{2}c(t)\]
 
−
 
−
                            <p><em>
 
−
                                Concentrations are denoted as the lowercase letter eg the concentration of \(S\) is given by \(s\).
 
−
                            </em></p>
 
−
 
 
                             <p>
 
                             <p>
−
                                 You may have spotted that the enzyme is not consumed in this series of reactions. Therefore the total amount of enzyme remains constant. We reflect this in the expression \(e_{T}=e+c\). We can now use this expression to eliminate \(e(t)\) from our model, leaving:
+
                                 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}{dt}s(t)=-k_{1}s(t)(e_{T}-c(t)) + k_{-1}c(t)\]
+
                                \[\dfrac{d[mRNA]}{dt}=K_{max}\dfrac{[Arab:AraC]^{n}}{K_{half}^{n}+[Arab:AraC]^{n}}-\gamma_{1}[mRNA]\]
−
 
+
−
                            \[\dfrac{d}{dt}c(t)=-k_{-1}c(t) + k_{1}s(t)(e_{T}-c(t)) - k_{2}c(t)\]
+
−
 
+
−
                            \[\dfrac{d}{dt}p(t) = k_{2}c(t)\]
+
  
 +
                                \[\dfrac{d\left[P\right]}{dt}=\alpha\left[mRNA\right]-\gamma_{2}\left[P\right]\]
 
                             <p>
 
                             <p>
−
                                 We have one further simplification to make. This is the rapid equilibrium approximation, by which we assume that equilibrium between \(s+e\) and \(c\) on a much faster time scale than the reaction of \(c\) to \(p\). With this approximation we can write the following equation:
+
                                 Where we define the symbols as:
 +
                                <table class="table table-striped">
 +
                                    <thead>
 +
                                        <th>Symbol</th>
 +
                                        <th>Definition</th>
 +
                                        <th>Initial Value/Literature Value</th>
 +
                                        <th>Fitted</th>
 +
                                    </thead>
 +
                                    <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\)</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>\(15ntd\: s^{-1}\)/length of sequence [<a href="#references">6</a>]</td>
 +
                                        <td>\(6.6ntd\: s^{-1}\)/length of sequence</td>
 +
                                    </tr>
 +
                                    <tr>
 +
                                        <td>\(\gamma_{1}\)</td>
 +
                                        <td><a class="definition" title="degradation rate" data-content="A measure of how quickly something is decreasing. For example, the amount of mRNA in a cell decreases as that cell divides and multiplies.">Degradation rate</a> of mRNA</td>
 +
                                        <td>\(5.13\times10^{-4}s^{-1}\) [<a href="#references">5</a>]</td>
 +
                                        <td>\(1.1\times10^{-2}s^{-1}\)</td>
 +
                                    </tr>
 +
                                    <tr>
 +
                                        <td>\(\gamma_{2}\)</td>
 +
                                        <td>Degradation rate of product</td>
 +
                                        <td>\(5.13\times10^{-4}s^{-1}\) [<a href="#references">5</a>]</td>
 +
                                        <td>\(1.1\times10^{-2}s^{-1}\)</td>
 +
                                    </tr>
 +
                                    <tr>
 +
                                        <td>\(K_{max}\)</td>
 +
                                        <td>Maximal transcription rate</td>
 +
                                        <td>\(50ntd\: s^{-1}\)/length of sequence [<a href="#references">6</a>]</td>
 +
                                        <td>\(47ntd\: s^{-1}\)/length of sequence</td>
 +
                                    </tr>
 +
                                    <tr>
 +
                                        <td>\(K_{half}\)</td>
 +
                                        <td>Half-maximal transcription rate</td>
 +
                                        <td>\(160\mu M\) [<a href="#references">7</a>]</td>
 +
                                        <td>\(100\mu M\)</td>
 +
                                    </tr>
 +
                                    <tr>
 +
                                        <td>\(n\)</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="#references">8</a>]</td>
 +
                                        <td>\(2.73\)</td>
 +
                                    </tr>
 +
                                </table>
 
                             </p>
 
                             </p>
−
 
−
                            \[0=-k_{-1}c(t) + k_{1}s(t)(e_{T}-c(t)) - k_{2}c(t)\]
 
−
 
 
                             <p>
 
                             <p>
−
                                 Leading to:
+
                                 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."> expression in <em>E. coli</em> to extract experimental values. Here are the details of our fitting function.
 
                             </p>
 
                             </p>
−
 
+
                             <div class="image image-full">
−
                             \[c^{*}(t)= \dfrac{k_{1}s(t)e_{T}}{k_{-1} + k_{2} + k_{1}s(t)}\]
+
                                 <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>
−
                            <p>
+
                                    Results showing GFP concentration as a function of time, matched to our <a class="definition" title="deterministic model" data-content="A deterministic model predicts a single outcome from a given set of circumstances.">deterministic model</a>. Errors are given to one standard deviation and an arbitrary scaling factor is included as a fitted parameter.
−
                                 Now we can see that \(c\) is no longer an independent variable, but instead tracks the other variables. This can be referred to as \(c\) being in quasi-steady state. By including this new expression into our model we are left with:
+
                                </p>
−
                            </p>
+
−
 
+
−
                            \[\dfrac{d}{dt}s(t) = -\dfrac{k_{2}k_{1}s(t)e_{T}}{k_{-1} + k_{2} + k_{1}s(t)}\]
+
−
 
+
−
                            \[\dfrac{d}{dt}p(t) = \dfrac{k_{2}k_{1}s(t)e_{T}}{k_{-1} + k_{2} + k_{1}s(t)}\]
+
−
 
+
−
                            <p>
+
−
                                With this we have an expression describing our enzyme catalysed system in the form of a single reaction. The rate of this reaction is known as a Michaelis-Menten rate law. By defining \(K_{max}=k_{2}e_{T}\) and \(K_{half}=\dfrac{k_{-1}+k_{2}}{k_{1}}\) we can express in the more familiar form:
+
−
                            </p>
+
−
 
+
−
                            \[rate \: of \: S \mathop{\rightarrow} P = K_{max}\dfrac{s}{K_{half} + s}\]
+
−
 
+
−
                            <p>
+
−
                                It is worth noting that there are a number of ways to derive this form, each involving different approximations. These methods may lead to the constants \(K_{max}\) and \(K_{half}\) being defined differently.
+
−
                            </p>
+
−
 
+
−
                        </div>
+
−
                        <div id="gene-expression-networks-stochastic">
+
−
                            <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 <a class="definition" title="stochastic" data-content="A stochastic model predicts a set of possible outcomes weighted by their likelihoods, or probabilities. It allows for an element of randomness.">stochastic</a> 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>
 
                             </div>
 
                             <p>
 
                             <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.
+
                                 We can now calculate the limiting concentrations that our products will be expressed at. There are mutliple products being expressed, each of different sequence lengths. Here is a table showing the relevant proteins and sequence lengths:
−
                            </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>
 
                             </p>
−
                             <ol>
+
                             <table class="table table-striped">
−
                                 <li>The <strong>time</strong> taken between reactions</li>
+
                                 <thead>
−
                                 <li>The <strong>reaction</strong> that took place in that time.</li>
+
                                    <tr>
−
                             </ol>
+
                                        <th>
 +
                                            Product
 +
                                        </th>
 +
                                        <th>
 +
                                            Sequence Length (/bp)
 +
                                        </th>
 +
                                    </tr>
 +
                                 </thead>
 +
                                <tr>
 +
                                    <td>
 +
                                        pBAD HisB DNase DsbA
 +
                                    </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>
 
                             <p>
−
                                 as neither of these could be determined by our initial method [<a href="#Ref2">2</a>, <a href="#Ref4">4</a>].
+
                                 We now can run our model of the system by solving the set of equations using the MATLAB 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:
 
                             </p>
 
                             </p>
 
                             <div class="image image-full">
 
                             <div class="image image-full">
−
                                 <img src="https://static.igem.org/mediawiki/2015/7/72/OxiGEM_Stochastic_Tutorial_02.png"/>
+
                                 <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>
−
                                 <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
+
                                 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.
 
                             </p>
 
                             </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>
 
−
                            <div class="image image-full">
 
−
                                <img src="https://static.igem.org/mediawiki/2015/d/da/OxiGEM_Stochastic_Tutorial_03.jpg"/>
 
−
                            </div>
 
 
                         </div>
 
                         </div>
 
                     </div>
 
                     </div>
−
                    <video class="image-massive" poster="https://static.igem.org/mediawiki/2015/6/6d/Ox_modelling_will_poster.png" controls>
 
−
                        <source src="https://static.igem.org/mediawiki/2015/b/b2/OxiGEM_Modelling_With_Will.mp4" type="video/mp4">
 
−
                        Your browser does not support the video tag.
 
−
                    </video>
 
 
                 </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 class="section" id="delivery">
 +
                         <h2>Delivery</h2>
 
                         <p>
 
                         <p>
−
                             In this section we will derive some of the equations used in the <a href="https://2015.igem.org/Team:Oxford/Modeling#delivery">Delivery Section</a>.
+
                             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>
 
                         </p>
−
                         <div id="diffusion-beads">
+
                         <div id="delivery-dispersin">
−
                             <h3>Diffusion out of our Beads</h3>
+
                             <h3>Dispersin B</h3>
 
                             <p>
 
                             <p>
−
                                 In order to find the convection mass transfer coefficient, \(K_m\) of our beads in water we needed derive a theoretical curve to fit our experimental data to. To do this, we made a number of approximations. Our beads are assumed to be spherical and of the diameter. We also used a lumped capacitance model, in which the concentration of substance is assumed to be uniform within the bead. We will also make use of the equivalence of heat and mass transfer.
+
                                 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.
 
                             </p>
 
                             </p>
 
                             <p>
 
                             <p>
−
                                 Our derivation starts from Newton's law of cooling
+
                                 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>
 +
                        </div>
 +
                        <div id="delivery-beads">
 +
                            <h3>Beads</h3>
 +
                            <div id="delivery-beads-diffusion">
 +
                                <h4>Diffusion</h4>
 +
                                <p>
 +
                                    The bead delivery system consists of our cells being contained in alginate spheres. Water is passed through the 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.
 +
                                </p>
 +
                                <p>
 +
                                    To determine the convection mass transfer coefficient 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 form for the concentration of crystal violet in the bulk water as a function of time:
 +
                                </p>
  
−
                            \[Q = hA\Delta T\]
+
                                \[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)\]
−
                            <p>
+
                                <table class="table table-striped">
−
                                Where \(Q\) is the flow of energy, \(h\) is the convection coefficient, and \(\Delta T\) is the difference in temperature accross the boundary. Or rather its mass equivalent
+
                                    <thead>
−
                            </p>
+
                                        <th>Symbol</th>
−
                            \[\dot{m}=K_{m}A_{b}\left(c_{b}-c_{f}\right)\]
+
                                        <th>Definition</th>
−
                            <p>
+
                                        <th>Value</th>
−
                                Where the mass flow rate, \(\dot{m}\), across a boundary is driven by the concentration difference across it. In our case the concentration difference is between the concentration of particles, \(c_f\) in the fluid and the beads, \(c_b\). Using the relationship
+
                                        <th>Units</th>
−
                            </p>
+
                                    </thead>
−
                            \[\dot{m}=V_{f}\dot{c_{f}}=-V_{b}\dot{c_{b}}\]
+
                                    <tr>
−
                            <p>
+
                                        <td>\(A_{b}\)</td>
−
                                With \(V_f\) and \(c_f\) as the concentration and volume of the fluid surrounding the beads, as well as \(A\), \(V_b\) and \(c_b\) representing the total area, volume and concentration of the beads themselves. We arrive at
+
                                        <td>Total surface area of the beads</td>
−
                            </p>
+
                                        <td>\(0.0238\)</td>
 +
                                        <td>\(m^{2}\)</td>
 +
                                    </tr>
 +
                                    <tr>
 +
                                        <td>\(V_{b}\)</td>
 +
                                        <td>Total volume of beads</td>
 +
                                        <td>\(1.3463\times10^{-5}\)</td>
 +
                                        <td>\(m^{3}\)</td>
 +
                                    </tr>
 +
                                    <tr>
 +
                                        <td>\(c_{bo}\)</td>
 +
                                        <td>Initial concentration in beads</td>
 +
                                        <td>\(0.02451107\)</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>Convection mass diffusion coefficient</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>
 +
                                    Our fitting script, detailed here, 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 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>
 +
                                    Dispersin B is a significantly larger molecule than crystal violet so this diffusion coefficient will not be close to that for Dispersin B. To correct this we need to make use of similarity. More specifically we take the Sherwood Numbers of the systems to be equal therefore:
 +
                                </p>
  
−
                            \[\dot{c_{b}}=-\dfrac{V_{f}}{V_{b}}\dot{c_{f}}\quad c_{b}=c_{bo}-\dfrac{V_{f}}{V_{b}}c_{f}\]
+
                                \[\left(\dfrac{K_{m}R}{D}\right)_{crystal violet} = \left(\dfrac{K_{m}R}{D}\right)_{Dispersin B}\]
  
−
                            <p>
+
                                <table class="table table-striped">
−
                                In this \(c_{b0}\) is the initial concentration of the beads. By substituting this into the mass equivalent of Newton's law of cooling we arrive at the expression
+
                                    <thead>
−
                            </p>
+
                                        <th>Symbol</th>
−
                            \[V_{f}\dot{c_{f}}=K_{m}A_{b}\left(c_{bo}-\left(1+\dfrac{V_{f}}{V_{b}}\right)c_{f}\right)\]
+
                                        <th>Definition</th>
−
                            <p>
+
                                        <th>Value</th>
−
                                This expression can be rearranged and integrated to with respect to time to give
+
                                        <th>Units</th>
−
                            </p>
+
                                    </thead>
−
                            \[c_{f}=\dfrac{c_{bo}}{\left(1+\frac{V_{f}}{V_{b}}\right)}-\dfrac{V_{f}}{K_{m}A_{b}\left(1+\frac{V_{f}}{V_{b}}\right)}\exp\left(-\dfrac{K_{m}A_{b}\left(1+\frac{V_{f}}{V_{b}}\right)\left(t+k\right)}{V_{f}}\right)\]
+
                                    <tr>
−
                            <p>
+
                                        <td>
−
                                where \(k\) constant that comes from integration.
+
                                            \(D_{crystal violet}\)
−
                            </p>
+
                                        </td>
−
                            <p>
+
                                        <td>
−
                                At \(t=0\) \(c_f=0\), leaving us with
+
                                            Mass diffusivity of crystal violet in water
−
                            </p>
+
                                        </td>
−
                            \[\dfrac{c_{bo}}{\left(1+\frac{V_{f}}{V_{b}}\right)}=\dfrac{V_{f}}{K_{m}A_{b}\left(1+\frac{V_{f}}{V_{b}}\right)}\exp\left(-\dfrac{K_{m}A_{b}\left(1+\frac{V_{f}}{V_{b}}\right)k}{V_{f}}\right)\]
+
                                        <td>
−
                            <p>
+
                                            \(2.8652\times10^{9}\)[<a href="#references">2</a>]
−
                                resulting in
+
                                        </td>
−
                            </p>
+
                                        <td>
−
                            \[k=-\dfrac{V_{f}}{K_{m}A_{b}\left(1+\frac{V_{f}}{V_{b}}\right)}\ln\left(\dfrac{K_{m}A_{b}c_{bo}}{V_{f}}\right)\]
+
                                            \(\mu m^{2}s^{-1}\)
−
                            <p>
+
                                        </td>
−
                                and therefore
+
                                    </tr>
−
                            </p>
+
                                    <tr>
−
                            \[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)}{V_{f}}t\right)\right)\]
+
                                        <td>
−
                             <p>
+
                                            \(D_{Dispersin B}\)
−
                                This shows how the concentration of fluid varies over time when there is a constant volume of fluid surrounding the beads. However, in the experiment, 1ml of water was removed every 10 minutes. We account for this by taking \(V_f=V_{fo}-\dfrac{1\times10^{-6}}{10}\). Where volumes are given in \(m^3\), time is in minutes and \(V_fo\) is the initial volume of the fluid surrounding the beads.
+
                                        </td>
−
                            </p>
+
                                        <td>
−
                            <p>
+
                                            Mass diffusivity of Dispersin B in water
−
                                Now we have derived the final expression, have a look at how this result was fitted to the experimental data to determine the value of \(K_m\) <a href="https://2015.igem.org/Team:Oxford/Modeling#delivery-beads-diffusion">here</a>.
+
                                        </td>
−
                            </p>
+
                                        <td>
−
                        </div>
+
                                            \(100\)[<a href="#references">3</a>]
−
                        <div id="diffusion-mass-exchanger">
+
                                        </td>
−
                            <h3>Mass Exchanger</h3>
+
                                        <td>
−
                            <p>
+
                                            \(\mu m^{2}s^{-1}\)
−
                                In our analysis of the use of beads in a mass exchange system we used the expression
+
                                        </td>
−
                            </p>
+
                                    </tr>
−
                            \[J = K_mA\dfrac{c_{fo}-c_{fi}}{\ln\left(\dfrac{c_{B}-c_{fi}}{c_{B}-c_{fo}}\right)}\]
+
                                    <tr>
−
                            <p>
+
                                        <td>
−
                                but where did this come from?
+
                                            \(R\)
−
                            </p>
+
                                        </td>
−
                            <p>
+
                                        <td>
−
                                For this derivation we need to take a look at the mass flow out of the beads across a differential area, \(dA\). To help visulise the system below is a diagram showing how the conentration of proteins varies through the mass exchanger.
+
                                            Radius of bead
−
                            </p>
+
                                        </td>
−
                            <div class="image image-full">
+
                                        <td>
−
                                <img src="https://static.igem.org/mediawiki/2015/e/e0/Ox_mass_exchange_daig.jpg"/>
+
                                            \(1.695\)
 +
                                        </td>
 +
                                        <td>
 +
                                            \(mm\)
 +
                                        </td>
 +
                                    </tr>
 +
                                </table>
 +
                                <p>
 +
                                    By rearranging this 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>
 +
 
 +
                                \[J = K_mA\dfrac{c_{fo}-c_{fi}}{\ln\left(\dfrac{c_{B}-c_{fi}}{c_{B}-c_{fo}}\right)}\]
 +
                                <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 packing efficiency 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>
 
−
                                The mass flow rate across a differential area, \(dJ\) can be expressed as
 
−
                            </p>
 
−
                            \[dJ=K_{m}\left(c_{B}-c_{f}\right)dA\]
 
−
                            <p>
 
−
                                This expression comes from the mass equivalent of Newton's law of cooling. We may use the relationship
 
−
                            </p>
 
−
                            \[dJ=Qdc_{f}\]
 
−
                            <p>
 
−
                                Where \(Q\) is the volume flow rate of water through the mass exchanger. By combining these 2 relationships we arrive at
 
−
                            </p>
 
−
                            \[d\Delta c=-\dfrac{dJ}{Q}=\dfrac{-K_{m}\left(\Delta c\right)dA}{Q}\]
 
−
                            \[d\Delta c=\dfrac{-K_{m}\left(\Delta c\right)dA}{Q}\]
 
−
                            <p>
 
−
                                In this expression \(\Delta c=c_b-c_f\) at a distance \(x\) from the start of the mass exchanger. Rearranging and integrating leads to
 
−
                            </p>
 
−
                            \[\ln\left(\dfrac{\Delta c_{i}}{\Delta c_{o}}\right)=\dfrac{K_{m}}{Q}A\]
 
−
                            <p>
 
−
                                Where \(\Delta c_{i}\) and \(\Delta c_{o}\) are the concentration differences at the inlet and outlet of the mass exchanger respectively. Thus \(\Delta c_{i} = c_B - c_{fi}\) and \(\Delta c_{o} = c_B - c_{fo}\).
 
−
                            </p>
 
−
                            <p>
 
−
                                The total mass flow, \(J\), can also be expressed as \(J=Q\left(c_{fo}-c_{fi}\right)\). Substituting this into the previous expression gives us
 
−
                            </p>
 
−
                            \[\ln\left(\dfrac{\Delta c_{i}}{\Delta c_{o}}\right)=\dfrac{K_{m}A\left(c_{fo}-c_{fi}\right)}{J}\]
 
−
                            <p>
 
−
                                Which can be rearranged into the form
 
−
                            </p>
 
−
                            \[J=K_{m}A\dfrac{\left(c_{fo}-c_{fi}\right)}{\ln\left(\dfrac{c_{B}-c_{fi}}{c_{B}-c_{fo}}\right)}\]
 
−
                            <p>
 
−
                                Comparing this to the form of the mass equivalent of Newton's law of cooling, \(J = K_mA\Delta c\), we see that we have derived an effective average concentration difference across the entire mass exchanger. This is known as the log mean concentration difference. We used this result to determine an estimate for the total exchange area required for our system, <a href="https://2015.igem.org/Team:Oxford/Modeling#delivery-beads-mass-exchange">here</a>.
 
−
                            </p>
 
 
                         </div>
 
                         </div>
 +
                    </div>
 +
                    <div class="section-spacer">
 +
                    </div>
 +
                    <div class="section" id="conclusion">
 +
                        <h2>Conclusion</h2>
 +
                        <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.
 +
                        </p>
 
                     </div>
 
                     </div>
 +
                </div>
 +
                <div class="section" id="references">
 +
                    <h2>References</h2>
 +
                    <ol class="references">
 +
                        <li>Jeffrey B. Kaplan; Dispersin B polypeptides and uses thereof. Patent PI 8580551, Nov 12, 2013</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>"Physical Biology of the Cell", Rob Phillips, Jane Kondev and Julie Theriot (2009). Page 110</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</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>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>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>
 +
                    </ol>
 
                 </div>
 
                 </div>
 
             </div>
 
             </div>
 
             <div class="col-md-3 contents-sidebar">
 
             <div class="col-md-3 contents-sidebar">
−
                 <ul id="sidebar" class="nav nav-stacked" data-spy="affix">
+
                 <ul id="sidebar" class="nav nav-stacked affix-top sm-hidden xs-hidden" data-spy="affix">
−
                    <li><a href="#introduction">Introduction</a></li>
+
 
                     <li>
 
                     <li>
−
                         <a href="#gene-expression-networks">Gene Expression Networks</a>
+
                         <a href="#introduction">Introduction</a>
 +
                    </li>
 +
                    <li>
 +
                        <a href="#characterising-our-cells">Gene expression rates</a>
 
                         <ul class="nav nav-stacked">
 
                         <ul class="nav nav-stacked">
−
                             <li><a href="#gene-expression-networks-mass-action">Mass Action</a></li>
+
                             <li><a href="#characterising-our-cells-arab">Arabinose-induced promoter</a></li>
−
                            <li><a href="#gene-expression-networks-michaelis">Michaelis-Menten Kinetics</a></li>
+
−
                            <li><a href="#gene-expression-networks-stochastic">Stochastic</a></li>
+
 
                         </ul>
 
                         </ul>
 
                     </li>
 
                     </li>
 
                     <li>
 
                     <li>
−
                         <a href="#diffusion">Diffusion</a>
+
                         <a href="#delivery">Delivery</a>
 
                         <ul class="nav nav-stacked">
 
                         <ul class="nav nav-stacked">
−
                             <li><a href="#diffusion-beads">Diffusion out of our Beads</a></li>
+
                             <li><a href="#delivery-dispersin">Dispersin B</a></li>
−
                             <li><a href="#diffusion-mass-exchanger">Mass Exchanger</a></li>
+
                             <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>
 +
                        <a href="#references">References</a>
 
                     </li>
 
                     </li>
 
                 </ul>
 
                 </ul>

Revision as of 16:06, 18 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, is known as a double repressor. AraC binds to pBAD which represses transcription of mRNA. By introducing arabinose into the system, AraC will bind to arabinose and form the Arab:AraC compound. Transcription can then 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 we define the symbols as:

Symbol Definition Initial Value/Literature Value Fitted
\([Arab:AraC]\) The concentration of associated arabinose and AraC \(0\) -
\([mRNA]\) The concentration of mRNA \(0\) -
\([P]\) The concentration of our product \(0\) -
\(\alpha\) Translation rate \(15ntd\: s^{-1}\)/length of sequence [6] \(6.6ntd\: s^{-1}\)/length of sequence
\(\gamma_{1}\) Degradation rate of mRNA \(5.13\times10^{-4}s^{-1}\) [5] \(1.1\times10^{-2}s^{-1}\)
\(\gamma_{2}\) Degradation rate of product \(5.13\times10^{-4}s^{-1}\) [5] \(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 expression in E. coli to extract experimental values. Here are the details of our fitting function.

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 at. There are mutliple products being expressed, each of different sequence lengths. 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 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 the 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.

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 form 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.0238\) \(m^{2}\)
\(V_{b}\) Total volume of beads \(1.3463\times10^{-5}\) \(m^{3}\)
\(c_{bo}\) Initial concentration in beads \(0.02451107\) \(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 diffusion 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.

Our fitting script, detailed here, 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 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 correct this we need to make use of similarity. More specifically we take the Sherwood Numbers of the systems to be equal therefore:

\[\left(\dfrac{K_{m}R}{D}\right)_{crystal violet} = \left(\dfrac{K_{m}R}{D}\right)_{Dispersin B}\]
Symbol Definition Value Units
\(D_{crystal violet}\) Mass diffusivity of crystal violet in water \(2.8652\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}\)
\(R\) Radius of bead \(1.695\) \(mm\)

By rearranging this 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