Difference between revisions of "Team:Oxford/Modeling"

Line 3: Line 3:
 
<html>
 
<html>
 
     <div class="container-fluid page-heading">
 
     <div class="container-fluid page-heading">
−
         <h3>Modeling</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>
−
                             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 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 focussed on building a better picture of the project to evaluate the effectiveness of initial designs, as well as to provide insight into how the system can (or must) be improved. To do this we have split our modelling efforts into three main sections: <em>Characterising our Cells</em>, <em>Interaction with the Environment</em>, and <em>Interaction with the Biofilm</em>. By combining the information gathered in each of these section we hope to ultimately answer the question: Is our system feasible? If not, where should the design be altered?
 
                         </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.
+
                             To help readers of all kinds and specialisations understand this page we have produced guides for all the modelling techniques used in this section which are available on the Modelling Tutorial page and will be linked to when relevant on this page.
 
                         </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>Characterising Our Cells</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?”, and “What are the main limiting rates/concentrations?”. The first will help us assess the feasibility of our project, ie are our cells too slow? The second will aid us in further optimising our design.
−
                         <div id="gene-expression-deterministic">
+
                        </p>
−
                             <h3>Deterministic</h3>
+
                         <div id="characterising-our-cells-arab">
 +
                             <h3>Arabinose Induced</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 arabinose induced promoter for the expression of a number of our proteins. This promoter can be modelled as the following chemical system:
 
                             </p>
 
                             </p>
−
                                \[\overset{\alpha_{1}}{\rightarrow}Arab\quad\overset{\alpha_{5}}{\rightarrow}AraC\]
+
                            \[\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\mathrel{\mathop{\rightleftharpoons}^{\mathrm{\alpha_{2}}}_{\mathrm{\alpha_{3}}}}Arab:AraC\]
  
−
                                \[\overset{K}{\rightarrow}mRNA\overset{\alpha_{4}}{\rightarrow}P\]
+
                            \[\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\]
+
                            \[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>]:
+
                                 You can find out more information about how this promoter works here.
 +
                            </p>
 +
                            <p>
 +
                                From this system of chemical reactions we can derive the following set of ODEs[<a href="#Ref4">4</a>]:
 
                             </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\right]}{dt}=\alpha_{1}+\alpha_{2}\left[Arab\right]\left[AraC\right]-\alpha_{3}\left[Arab:AraC\right]-\gamma_{1}\left[Arab\right]\]
Line 42: Line 46:
 
                                 \[\dfrac{d\left[Arab:AraC\right]}{dt}=\alpha_{3}\left[Arab:AraC\right]-\alpha_{2}\left[Arab\right]\left[AraC\right]\]
 
                                 \[\dfrac{d\left[Arab:AraC\right]}{dt}=\alpha_{3}\left[Arab:AraC\right]-\alpha_{2}\left[Arab\right]\left[AraC\right]\]
  
−
                                 \[\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[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[mRNA]}{dt}=K_{max}\dfrac{[Arab:AraC]^{n}}{K_{half}^{n}+[Arab:AraC]^{n}}-\gamma_{3}[mRNA]\]
 
                                 \[\dfrac{d[mRNA]}{dt}=K_{max}\dfrac{[Arab:AraC]^{n}}{K_{half}^{n}+[Arab:AraC]^{n}}-\gamma_{3}[mRNA]\]
Line 48: Line 52:
 
                                 \[\dfrac{d\left[P\right]}{dt}=\alpha_{4}\left[mRNA\right]-\gamma_{4}\left[P\right]\]
 
                                 \[\dfrac{d\left[P\right]}{dt}=\alpha_{4}\left[mRNA\right]-\gamma_{4}\left[P\right]\]
 
                             <p>
 
                             <p>
−
                                 Where we define the symbols
+
                                 Where we define the symbols as:
−
                                 <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>\([Arab]\)</td>
 
                                         <td>The concentration of Arabinose</td>
 
                                         <td>The concentration of Arabinose</td>
−
                                         <td>\(1\times10^{-4}M\)</td>
+
                                         <td>\(1\times10^{-5}M\)</td>
 
                                         <td>-</td>
 
                                         <td>-</td>
 
                                     </tr>
 
                                     </tr>
Line 161: Line 165:
 
                             </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. Later we will fit the parameters to the experimental data found by the wet lab team. For now though we can plot the expression graphs using the literature values. This will provide an estimate to the time scales involved.
 
                             </p>
 
                             </p>
−
                            <div class="image image-full">
 
−
                                <img src="https://static.igem.org/mediawiki/2015/7/74/OxiGEM_Modeling_Arabinose_01.png"/>
 
−
                            </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>].
+
                                 There are mutliple products being expressed using this inducer-promoter pair, each of different sequence lenghts. 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
 +
                                        </th>
 +
                                        <th>
 +
                                            Sequence Length (/bp)
 +
                                        </th>
 +
                                    </tr>
 +
                                 </thead>
 +
                                <tr>
 +
                                     <td>
 +
                                        pBAD HisB DNase DsbA
 +
                                    </td>
 +
                                    <td>
 +
                                        621
 +
                                    </td>
 +
                                </tr>
 +
                                <tr>
 +
                                    <td>
 +
                                        pBAD HisB DspB YebF
 +
                                    </td>
 +
                                    <td>
  
−
                                     <li>Find propensities \(\alpha_{i}\) of \(N\) possible reactions. Find \(\alpha_{0}=\sum\limits_{i=1}^N \alpha_{i}\).</li>
+
                                     </td>
 +
                                </tr>
 +
                                <tr>
 +
                                    <td>
 +
                                        pBAD HisB DspB
 +
                                    </td>
 +
                                    <td>
  
−
                                     <li>Compute time \(\tau=\dfrac{1}{\alpha_{0}}\ln\left(\dfrac{1}{r_{1}}\right)\).</li>
+
                                     </td>
 +
                                </tr>
 +
                                <tr>
 +
                                    <td>
 +
                                        pBAD HisB MccS
 +
                                    </td>
 +
                                    <td>
 +
                                        414
 +
                                    </td>
 +
                                </tr>
 +
                                <tr>
 +
                                    <td>
 +
                                        pBAD HisB DspB Fla
 +
                                    </td>
 +
                                    <td>
  
−
                                     <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>
+
                                     </td>
−
 
+
                                </tr>
−
                                     <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 Art-175 DsbA
−
                                 </ol>
+
                                    </td>
−
 
+
                                    <td>
−
                             </div>
+
                                        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>
−
                                 Turning the reaction rates in our model into propensities is simple - for example, the propensity for the association reaction
+
                                 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:
−
                                \(Arab+AraC\rightarrow Arab:AraC\) is \(\alpha_{3}\left[Arab\right]\left[AraC\right]\).
+
 
                             </p>
 
                             </p>
 +
                            <div class="image image-full">
 +
                                <img src="https://static.igem.org/mediawiki/2015/f/f6/Ox_arab_induced_proteins.png"/>
 +
                            </div>
 
                         </div>
 
                         </div>
 
                     </div>
 
                     </div>
−
                </div>
+
                    <div class="section-spacer"></div>
−
                <div class="section-spacer"></div>
+
                    <div class="section" id="interaction-with-environment">
−
                <div class="section" id="diffusion">
+
                        <h2>Interaction with the Environment</h2>
−
                    <div class="slim">
+
                         <p>
−
                         <h2>Diffusion</h2>
+
                            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. This involves modelling the diffusion of the products in different topologies, each associated with a potential physical design of the catheter. With this information we can provide a better estimate of the time scale that our project is working on and assess any need for optimisation.
 +
                        </p>
 +
                        <p>
 +
                            In addition, one of our systems relies on the detection of the biofilm to cause lysis. We will look at the how the quorum sensing signal moves to the cells.
 +
                        </p>
 
                     </div>
 
                     </div>
−
                     <div id="diffusion-deterministic">
+
                     <div class="section-spacer"></div>
−
                        <div class="slim">
+
                    <div class="section" id="interaction-with-biofilm">
−
                            <h3>Deterministic</h3>
+
                        <h2>Interaction with the Biofilm</h2>
−
                            <p>
+
                        <p>
−
                                The standard from for the diffusion equation,
+
                             Once the antibiofilm agents have arrived at the biofilm, what impact do they have? We will look into the required concentrations and rates of the antibiofilm agents at the biofilm to overcome its rate of growth.
−
                             </p>
+
                        </p>
−
                                \[\dfrac{d[P]}{dt}=D\dfrac{d^{2}[P]}{dx^{2}}\]
+
−
                            <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:
+
−
                            </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>
+
−
                        </div>
+
 
                     </div>
 
                     </div>
−
                     <div id="diffusion-stochastic">
+
                     <div class="section-spacer"></div>
−
                        <div class="slim">
+
                    <div class="section" id="overall">
−
                            <h3>Stochastic</h3>
+
                        <h2>Combined: What do we know about the system as a whole?</h2>
−
                            <p>
+
                        <p>
−
                                We have the option of using three stochastic models for diffusion - random walk, random velocity and compartmental Gillespie.
+
                             What we write here depends on the results of the previous section.
−
                            </p>
+
                        </p>
−
                            <h4>Random Walk</h4>
+
                        <div id="overall-time">
−
                            <p>
+
                             <h3>How quickly will the system run?</h3>
−
                                We initially take a collection of particles and, for each one, compute the number
+
−
                             </p>
+
−
                                \[dr=\sqrt{2Ddt}\xi_{r}\]
+
−
                            <p>
+
−
                                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>
+
−
                            <h4>Random Velocity</h4>
+
−
                            <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>
+
−
 
+
−
                                    <li>Generate a random number \(r\) uniformly distributed between 0 and 1.</li>
+
−
 
+
−
                                    <li>If \(r<\dfrac{s^{2}}{2D}\Delta t\), calculate a new random direction of velocity for the particle.</li>
+
−
 
+
−
                                    <li>Repeat steps 2-5 for each particle until a certain time or computational budget is exhausted.</li>
+
−
                                </ol>
+
−
                             </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 id="overall-quantity">
−
                </div>
+
                             <h3>How many engineered bacteria do we need?</h3>
−
                <div class="section-spacer"></div>
+
−
                <div class="section" id="modeling-tutorial">
+
−
                    <div class="slim">
+
−
                         <h2>Modeling Tutorial</h2>
+
−
                    </div>
+
−
                    <div id="modeling-tutorial-deterministic">
+
−
                        <div class="slim">
+
−
                             <h3>Deterministic</h3>
+
−
                            <p>
+
−
                                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>
−
                    </div>
+
                        <div id="overall-concentration">
−
                    <div id="modeling-tutorial-stochastic">
+
                             <h3>What is the minimum steady state concentration we need to reach to overcome the growth rate of the biofilm?</h3>
−
                        <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>
Line 370: Line 314:
 
                 <div class="section" id="references">
 
                 <div class="section" id="references">
 
                     <h2>References</h2>
 
                     <h2>References</h2>
−
                     <ol class="references">
+
                     <p>
−
                        <li>Synthetic Biology - A Primer, Paul S. Freemont, Richard I Kenne et al., 2012 Imperial College Press<a name="Ref1"></a></li>
+
                    </p>
−
 
+
−
                        <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>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>Brian Ingalls,Mathematical Modeling in Systems Biology: An Introduction, 2012 University of Waterloo<a name="Ref4"></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<a name="Ref5"></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<a name="Ref6"></a></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<a name="Ref8"></a></li>
+
−
                    </ol>
+
 
                 </div>
 
                 </div>
 
             </div>
 
             </div>
Line 395: Line 324:
 
                     </li>
 
                     </li>
 
                     <li>
 
                     <li>
−
                         <a href="#gene-expression">Gene Expression</a>
+
                         <a href="#characterising-our-cells">Characterising Our Cells</a>
 
                         <ul class="nav nav-stacked">
 
                         <ul class="nav nav-stacked">
−
                             <li><a href="#gene-expression-deterministic">Deterministic</a></li>
+
                             <li><a href="#characterising-our-cells-arab">Arabinose Induced</a></li>
−
                            <li><a href="#gene-expression-stochastic">Stochastic</a></li>
+
 
                         </ul>
 
                         </ul>
 
                     </li>
 
                     </li>
 
                     <li>
 
                     <li>
−
                         <a href="#diffusion">Diffusion</a>
+
                         <a href="#interaction-with-environment">Interaction with the Environment</a>
−
                        <ul class="nav nav-stacked">
+
                    </li>
−
                            <li><a href="#diffusion-deterministic">Deterministic</a></li>
+
                    <li>
−
                            <li><a href="#diffusion-stochastic">Stochastic</a></li>
+
                        <a href="#interaction-with-biofilm">Interaction with the Biofilm</a>
−
                        </ul>
+
 
                     </li>
 
                     </li>
 
                     <li>
 
                     <li>
−
                         <a href="#modeling-tutorial">Modeling Tutorial</a>
+
                         <a href="#overall">Combined: What do we know about the system as a whole?</a>
 
                         <ul class="nav nav-stacked">
 
                         <ul class="nav nav-stacked">
−
                             <li><a href="#modeling-tutorial-deterministic">Deterministic</a></li>
+
                             <li><a href="#overall-time">How quickly will the system run?</a></li>
−
                             <li><a href="#modeling-tutorial-stochastic">Stochastic</a></li>
+
                             <li><a href="#overall-quantity">How many engineered bacteria do we need?</a></li>
−
                             <li><a href="#modeling-tutorial-finite-difference">Finite Difference Models</a></li>
+
                             <li><a href="#overall-concentration">What is the minimum steady state concentration we need to reach to overcome the growth rate of the biofilm?</a></li>
 
                         </ul>
 
                         </ul>
 
                     </li>
 
                     </li>

Revision as of 12:57, 4 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 focussed on building a better picture of the project to evaluate the effectiveness of initial designs, as well as to provide insight into how the system can (or must) be improved. To do this we have split our modelling efforts into three main sections: Characterising our Cells, Interaction with the Environment, and Interaction with the Biofilm. By combining the information gathered in each of these section we hope to ultimately answer the question: Is our system feasible? If not, where should the design be altered?

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

Characterising Our Cells

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?”, and “What are the main limiting rates/concentrations?”. The first will help us assess the feasibility of our project, ie are our cells too slow? The second will aid us in further optimising our design.

Arabinose Induced

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:

\[\overset{\alpha_{1}}{\rightarrow}Arab\quad\overset{\alpha_{5}}{\rightarrow}AraC\] \[Arab+AraC\mathrel{\mathop{\rightleftharpoons}^{\mathrm{\alpha_{2}}}_{\mathrm{\alpha_{3}}}}Arab:AraC\] \[\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\]

You can find out more information about how this promoter works here.

From this system of chemical reactions we can derive the following set of ODEs[4]:

\[\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\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[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]\]

Where we define the symbols as:

Symbol Definition Initial Value/Literature Value Fitted
\([Arab]\) The concentration of Arabinose \(1\times10^{-5}M\) -
\([AraC]\) The concentration of AraC \(1\times10^{-5}M\) -
\([Arab:AraC]\) The concentration of associated Arabinose and AraC \(0\) -
\([mRNA]\) The concentration of mRNA \(0\) -
\([P]\) The concentration of our product \(0\) -
\(\alpha_{1}\) Basal production of Arabinose ??? ?
\(\alpha_{2}\) Association constant \(2.8\times10^{7}s^{-1}\) [7] ?
\(\alpha_{3}\) Dissociation constant \(0.022s^{-1}\) [7] ?
\(\alpha_{4}\) Translation rate \(15ntd\: s^{-1}\)/length of sequence [6] ?
\(\alpha_{5}\) Basal production of AraC ??? ?
\(\gamma_{1}\) Degradation rate of Arabinose \(5.13\times10^{-4}s^{-1}\) [5] ?
\(\gamma_{2}\) Degradation rate of AraC \(5.13\times10^{-4}s^{-1}\) [5] ?
\(\gamma_{3}\) Degradation rate of mRNA \(5.13\times10^{-4}s^{-1}\) [5] ?
\(\gamma_{4}\) Degradation rate of product \(5.13\times10^{-4}s^{-1}\) [5] ?
\(K_{max}\) Maximal transcription rate \(50ntd\: s^{-1}\)/length of sequence [6] ?
\(K_{half}\) Half-maximal transcription rate \(160\mu M\) [8] ?
\(n\) Hill coefficient \(2.65\) [3] ?

This table contains literature values for the parameters, found from a number of sources. Later we will fit the parameters to the experimental data found by the wet lab team. For now though we can plot the expression graphs using the literature values. This will provide an estimate to the time scales involved.

There are mutliple products being expressed using this inducer-promoter pair, each of different sequence lenghts. Here is a table showing the relevant proteins and sequence lengths:

Product Sequence Length (/bp)
pBAD HisB DNase DsbA 621
pBAD HisB DspB YebF
pBAD HisB DspB
pBAD HisB MccS 414
pBAD HisB DspB Fla
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:

Interaction with the Environment

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. This involves modelling the diffusion of the products in different topologies, each associated with a potential physical design of the catheter. With this information we can provide a better estimate of the time scale that our project is working on and assess any need for optimisation.

In addition, one of our systems relies on the detection of the biofilm to cause lysis. We will look at the how the quorum sensing signal moves to the cells.

Interaction with the Biofilm

Once the antibiofilm agents have arrived at the biofilm, what impact do they have? We will look into the required concentrations and rates of the antibiofilm agents at the biofilm to overcome its rate of growth.

Combined: What do we know about the system as a whole?

What we write here depends on the results of the previous section.

How quickly will the system run?

How many engineered bacteria do we need?

What is the minimum steady state concentration we need to reach to overcome the growth rate of the biofilm?

References