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>Practices</h3>
+
         <h3>Modeling</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="practices-introduction">
+
                 <div class="section" id="introduction">
 
                     <div class="slim">
 
                     <div class="slim">
 
                         <h2>Introduction</h2>
 
                         <h2>Introduction</h2>
 
                         <p>
 
                         <p>
−
                             Our project relies on a three way conversation between the team, the public and experts. It touches every aspect of the project, from our choice of application to the details of our delivery system. We promoted Synthetic Biology and iGEM through outreach programs to inspire the next generation.
+
                             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.
 
                         </p>
 
                         </p>
 
                         <p>
 
                         <p>
−
                             Urinary tract infections are a huge problem globally with millions of cases reported each year. We’re producing a guide for everything you need to know about urinary tract infections, as well as a treatment to beat antibiotics, which are rapidly becoming ineffective.
+
                             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.
 
                         </p>
 
                         </p>
 
                     </div>
 
                     </div>
 
                 </div>
 
                 </div>
−
                 <div class=”section-spacer”></div>
+
                 <div class="section-spacer"></div>
−
                 <div class="section" id="practices-public-interaction">
+
                 <div class="section" id="gene-expression">
 
                     <div class="slim">
 
                     <div class="slim">
−
                         <h2>Public Interaction</h2>
+
                         <h2>Gene Expression</h2>
−
                    </div>
+
                        <div id="gene-expression-deterministic">
−
                    <div id="practices-public-interaction-project-choice">
+
                             <h3>Deterministic</h3>
−
                        <div class="slim">
+
−
                             <h3>Project Choice</h3>
+
 
                             <p>
 
                             <p>
−
                                 To decide on our project idea, we sent out an initial questionnaire to the public to hear about what they thought about synthetic biology. We asked what big problems they wanted solving. We took the questionnaire to schools, to the streets and to our friends.
+
                                 A deterministic model evolves according to ordinary or partial differential equations. We model the following system:
 
                             </p>
 
                             </p>
 +
                                \[\overset{\alpha_{1}}{\rightarrow}Arab\quad\overset{\alpha_{5}}{\rightarrow}AraC\]
 +
 +
                                \[Arab+AraC\mathrel{\mathop{\rightleftharpoons}^{\mathrm{\alpha_{2}}}_{\mathrm{\alpha_{3}}}}Arab:AraC\]
 +
 +
                                \[\overset{K}{\rightarrow}mRNA\overset{\alpha_{4}}{\rightarrow}P\]
 +
 +
                                \[Arab\overset{\gamma_{1}}{\rightarrow}\phi\quad AraC\overset{\gamma_{2}}{\rightarrow}\phi\quad mRNA\overset{\gamma_{3}}{\rightarrow}\phi\quad P\overset{\gamma_{4}}{\rightarrow}\phi\]
 
                             <p>
 
                             <p>
−
                                 Examples of their suggestions for the applications of synthetic biology include bacteria which:
+
                                 Using the following deterministic equations [<a href="#Ref4">4</a>]:
 
                             </p>
 
                             </p>
−
                            <ul>
+
                                 \[\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]\]
−
                                 <li>Remove carbon dioxide from the atmosphere</li>
+
 
−
                                <li>Target and kill cancerous cells</li>
+
                                 \[\dfrac{d\left[Arab:AraC\right]}{dt}=\alpha_{3}\left[Arab:AraC\right]-\alpha_{2}\left[Arab\right]\left[AraC\right]\]
−
                                 <li>Help treat Alzheimer's disease</li>
+
 
−
                                <li>Produce energy</li>
+
                                 \[\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]\]
−
                                 <li>Sew up holes in clothes</li>
+
 
−
                                <li>Produce teeth glue</li>
+
                                 \[\dfrac{d[mRNA]}{dt}=K_{max}\dfrac{[Arab:AraC]^{n}}{K_{half}^{n}+[Arab:AraC]^{n}}-\gamma_{3}[mRNA]\]
−
                                 <li>Indicate how long someone has been dead for</li>
+
 
−
                                 <li>Combat antibiotic resistance</li>
+
                                 \[\dfrac{d\left[P\right]}{dt}=\alpha_{4}\left[mRNA\right]-\gamma_{4}\left[P\right]\]
−
                            </ul>
+
 
                             <p>
 
                             <p>
−
                                 Of our responses, around 40 were related to Medicine and Health [<a href="#PRef1">1</a>]. This led us to choose that track for our project. However, it was our team member George Driscoll’s work at the UTI clinic in London which helped us to select UTIs as a specific cause. Due to the un-aesthetic nature of the infection, it often receives less attention with regard to research.
+
                                 Where we define the symbols
 +
                                <table border="1">
 +
                                    <tr>
 +
                                        <th>Symbol</th>
 +
                                        <th>Definition</th>
 +
                                        <th>Initial Value/Literature Value</th>
 +
                                        <th>Fitted</th>
 +
                                    </tr>
 +
                                    <tr>
 +
                                        <td>\([Arab]\)</td>
 +
                                        <td>The concentration of Arabinose</td>
 +
                                        <td>\(1\times10^{-4}M\)</td>
 +
                                        <td>-</td>
 +
                                    </tr>
 +
                                    <tr>
 +
                                        <td>\([AraC]\)</td>
 +
                                        <td>The concentration of AraC</td>
 +
                                        <td>\(1\times10^{-5}M\)</td>
 +
                                        <td>-</td>
 +
                                    </tr>
 +
                                    <tr>
 +
                                        <td>\([Arab:AraC]\)</td>
 +
                                        <td>The concentration of associated Arabinose and AraC</td>
 +
                                        <td>\(0\)</td>
 +
                                        <td>-</td>
 +
                                    </tr>
 +
                                    <tr>
 +
                                        <td>\([mRNA]\)</td>
 +
                                        <td>The concentration of mRNA</td>
 +
                                        <td>\(0\)</td>
 +
                                        <td>-</td>
 +
                                    </tr>
 +
                                    <tr>
 +
                                        <td>\([P]\)</td>
 +
                                        <td>The concentration of our product</td>
 +
                                        <td>\(0\)</td>
 +
                                        <td>-</td>
 +
                                    </tr>
 +
                                    <tr>
 +
                                        <td>\(\alpha_{1}\)</td>
 +
                                        <td>Basal production of Arabinose</td>
 +
                                        <td>???</td>
 +
                                        <td>?</td>
 +
                                    </tr>
 +
                                    <tr>
 +
                                        <td>\(\alpha_{2}\)</td>
 +
                                        <td>Association constant</td>
 +
                                        <td>\(2.8\times10^{7}s^{-1}\) [<a href="#Ref7">7</a>]</td>
 +
                                        <td>?</td>
 +
                                    </tr>
 +
                                    <tr>
 +
                                        <td>\(\alpha_{3}\)</td>
 +
                                        <td>Dissociation constant</td>
 +
                                        <td>\(0.022s^{-1}\) [<a href="#Ref7">7</a>]</td>
 +
                                        <td>?</td>
 +
                                    </tr>
 +
                                    <tr>
 +
                                        <td>\(\alpha_{4}\)</td>
 +
                                        <td>Translation rate</td>
 +
                                        <td>\(15ntd\: s^{-1}\)/length of sequence [<a href="#Ref6">6</a>]</td>
 +
                                        <td>?</td>
 +
                                    </tr>
 +
                                    <tr>
 +
                                        <td>\(\alpha_{5}\)</td>
 +
                                        <td>Basal production of AraC</td>
 +
                                        <td>???</td>
 +
                                        <td>?</td>
 +
                                    </tr>
 +
                                    <tr>
 +
                                        <td>\(\gamma_{1}\)</td>
 +
                                        <td> Degradation rate of Arabinose </td>
 +
                                        <td>\(5.13\times10^{-4}s^{-1}\) [<a href="#Ref5">5</a>]</td>
 +
                                        <td>?</td>
 +
                                    </tr>
 +
                                    <tr>
 +
                                        <td>\(\gamma_{2}\)</td>
 +
                                        <td>Degradation rate of AraC</td>
 +
                                        <td>\(5.13\times10^{-4}s^{-1}\) [<a href="#Ref5">5</a>]</td>
 +
                                        <td>?</td>
 +
                                    </tr>
 +
                                    <tr>
 +
                                        <td>\(\gamma_{3}\)</td>
 +
                                        <td>Degradation rate of mRNA</td>
 +
                                        <td>\(5.13\times10^{-4}s^{-1}\) [<a href="#Ref5">5</a>]</td>
 +
                                        <td>?</td>
 +
                                    </tr>
 +
                                    <tr>
 +
                                        <td>\(\gamma_{4}\)</td>
 +
                                        <td>Degradation rate of product</td>
 +
                                        <td>\(5.13\times10^{-4}s^{-1}\) [<a href="#Ref5">5</a>]</td>
 +
                                        <td>?</td>
 +
                                    </tr>
 +
                                    <tr>
 +
                                        <td>\(K_{max}\)</td>
 +
                                        <td>Maximal transcription rate</td>
 +
                                        <td>\(50ntd\: s^{-1}\)/length of sequence [<a href="#Ref6">6</a>]</td>
 +
                                        <td>?</td>
 +
                                    </tr>
 +
                                    <tr>
 +
                                        <td>\(K_{half}\)</td>
 +
                                        <td>Half-maximal transcription rate</td>
 +
                                        <td>\(160\mu M\) [<a href="#Ref8">8</a>]</td>
 +
                                        <td>?</td>
 +
                                    </tr>
 +
                                    <tr>
 +
                                        <td>\(n\)</td>
 +
                                        <td>Hill coefficient</td>
 +
                                        <td>\(2.65\) [<a href="#Ref3">3</a>]</td>
 +
                                        <td>?</td>
 +
                                    </tr>
 +
                                </table>
 
                             </p>
 
                             </p>
 
                             <p>
 
                             <p>
−
                                 A large proportion of our responses expressed concern for how Synthetic Biology would be used in society, with several references to the issues of contamination and exploitation for profit. With this in mind, we constructed a second questionnaire about our project, to test whether the public would get behind it.
+
                                 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.
 
                             </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>
−
                    </div>
+
                        <div id="gene-expression-stochastic">
−
                    <div id="practices-public-interaction-initial-feedback">
+
                             <h3>Stochastic</h3>
−
                        <div class="slim">
+
−
                             <h3>Initial Feedback</h3>
+
 
                             <p>
 
                             <p>
−
                                 We sent a second questionnaire to find out more about whether the public would use a Solution from synthetic biology to treat Urinary tract infections. We asked more about whether they had heard of genetic engineering or synthetic biology, and how much they trust a recommended treatment by a doctor. In collaboration with UCL, we also filmed some of these responses on the street. The results were overwhelmingly positive.
+
                                 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>].
 +
                            </p>
 +
                            <div class="algorithm">
 +
                                <h3>The Gillespie Algorithm</h3>
 +
                                <ol>
 +
                                    <li>Find random numbers \(r_{1}\) and \(r_{2}\) which are uniformly distributed between 0 and 1.</li>
 +
 
 +
                                    <li>Find propensities \(\alpha_{i}\) of \(N\) possible reactions. Find \(\alpha_{0}=\sum\limits_{i=1}^N \alpha_{i}\).</li>
 +
 
 +
                                    <li>Compute time \(\tau=\dfrac{1}{\alpha_{0}}\ln\left(\dfrac{1}{r_{1}}\right)\).</li>
 +
 
 +
                                    <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>
 +
 
 +
                                    <li>Implement reaction \(j\) and increment time \(t\) to \(t+\tau\). </li>
 +
 
 +
                                    <li>Repeat steps 1-5 until a sufficient time or computation budget has been exhausted.</li>
 +
                                </ol>
 +
 
 +
                            </div>
 +
                            <p>
 +
                                Turning the reaction rates in our model into propensities is simple - for example, the propensity for the association reaction
 +
                                \(Arab+AraC\rightarrow Arab:AraC\) is \(\alpha_{3}\left[Arab\right]\left[AraC\right]\).
 
                             </p>
 
                             </p>
 
                         </div>
 
                         </div>
Line 61: Line 197:
 
                 </div>
 
                 </div>
 
                 <div class="section-spacer"></div>
 
                 <div class="section-spacer"></div>
−
                 <div class="section" id="practices-medical-professionals">
+
                 <div class="section" id="diffusion">
 
                     <div class="slim">
 
                     <div class="slim">
−
                         <h2>Talking to Medical Professionals</h2>
+
                         <h2>Diffusion</h2>
−
                        <p>
+
                    </div>
−
                            Now that we were confident that the public would benefit from our idea, we set out designing our constructs. We visited nurses at the JR hospital in Oxford and the UTI clinic in London to get the Doctors on board; this discussion is essential for the success of the project.
+
                    <div id="diffusion-deterministic">
−
                        </p>
+
                         <div class="slim">
−
                         <div id="practices-medical-professionals-churchill">
+
                             <h3>Deterministic</h3>
−
                             <h3>Churchill Hospital, Oxford</h3>
+
 
                             <p>
 
                             <p>
−
                                 Our first visit to the hospital was to the outpatient clinic during which we spoke with Jan, one of the nurses on the ward. Jan told us about a case of a person getting septicaemia as a result of a urinary infection. The patient had received antibiotics for seven days and had come back for a check up. Their urine sample was clear and all seemed fine but then the patient had started to shake. The bacteria were now in their blood as it had travelled back up the ureter to the kidney. Even though this was a rare case, it was shocking to hear about such a serious case and made our project feel very relevant.
+
                                 The standard from for the diffusion equation,
 
                             </p>
 
                             </p>
 +
                                \[\dfrac{d[P]}{dt}=D\dfrac{d^{2}[P]}{dx^{2}}\]
 
                             <p>
 
                             <p>
−
                                 Jan also made the following points:
+
                                 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>
−
                            <ul>
+
                                 \[[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)\]
−
                                 <li>People with infections have a catheter because they need a way to empty the bladder; else the urine travels up the ureter and back into the bladder</li>
+
−
                                <li>If a patient becomes septic the catheter has to be removed or can be fatal</li>
+
−
                                <li>UTIs are not just contracted by the catheter and it is important to also consider community based UTIs</li>
+
−
                                <li>“UTIs are very common and can be quite painful”</li>
+
−
                                <li>No separate ward for UTIs – they are treated in every ward</li>
+
−
                                <li>The protocol for treatment is to take a urine sample, see if there is an infection, and prescribe antibiotics that the bacteria are most sensitive to</li>
+
−
                                <li>Elderly hospital wards are likely to have many cases of UTIs</li>
+
−
                            </ul>
+
 
                             <p>
 
                             <p>
−
                                 We took a lot from this initial conversation. We went onto investigating the pros and cons of the current methods of treating urinary infections and compared these to what Solution could offer. We realized that we needed to consider the catheter more from a hospital/medical perspective as up to this point we had confused its function, thinking it was more to do with administering medication rather than emptying the bladder. Following this meeting, the design of the catheter became an integral part of the project.
+
                                 where \(x_{i}\) is a generalised position co-ordinate and \(N\) is the dimensionality of the system. This solution only applies to a source of [P] at the position \(x=0\). We can also solve for a grid of 1 or 2D points using finite difference models given in the <a href="#modeling-tutorial">tutorial section</a>.
 
                             </p>
 
                             </p>
 
                         </div>
 
                         </div>
−
                        <div id="practices-medical-professionals-jr">
+
                    </div>
−
                             <h3>The John Radcliffe Hospital, Oxford</h3>
+
                    <div id="diffusion-stochastic">
 +
                        <div class="slim">
 +
                             <h3>Stochastic</h3>
 
                             <p>
 
                             <p>
−
                                 We still wanted to learn more about urinary infections as well as to get some feedback from nurses abourtour idea. We organized a trip to the Adams Ward (Geratology) to learn more about how urinary infections affect elderly people.
+
                                 We have the option of using three stochastic models for diffusion - random walk, random velocity and compartmental Gillespie.
 
                             </p>
 
                             </p>
−
                             <div class="algorithm" id="practices-medical-professionals-jr-laura-1">
+
                             <h4>Random Walk</h4>
−
                                <h3>First interview with Laura Evans, Adams Ward</h3>
+
                            <p>
−
                                 <p>
+
                                We initially take a collection of particles and, for each one, compute the number
−
                                    What is the procedure for treating UTIs?
+
                            </p>
−
                                </p>
+
                                 \[dr=\sqrt{2Ddt}\xi_{r}\]
−
                                 <ol class="interview-response">
+
                            <p>
−
                                    <li><em>Dip urine</em></li>
+
                                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\)].
−
                                     <li><em>If the test comes back as positive, treat with a wide spectrum antibiotic</em></li>
+
                            </p>
−
                                     <li><em>Whether or not the catheter is inserted with prophylactic antibiotic treatment is the doctor’s decision</em></li>
+
                            <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>
 
                                 </ol>
−
                                <p>
 
−
                                    What happens when a patient tests positive for a urinary infection?
 
−
                                </p>
 
−
                                <p class="interview-response">
 
−
                                    <em>Whether or not the catheter is removed if a patient tests positive for a urinary infection depends largely on the reason that the catheter has been fitted. On the whole, the catheter remains fitted and the patient is treated with a large dose of antibiotics.</em>
 
−
                                </p>
 
−
                                <p>
 
−
                                    Is antibiotic resistance a problem?
 
−
                                </p>
 
−
                                <p class="interview-response">
 
−
                                    <em>Yes, particularly on this ward. As we treat elderly patients with recurring infections, their UTIs are frequently resistant to antibiotic treatment. We try different combinations of antibiotics but recurring infections are a significant problem.</em>
 
−
                                </p>
 
−
                                <p>
 
−
                                    Our project involves designing a catheter that prevents the formation of a biofilm on its surface. What do you think of this idea?
 
−
                                </p>
 
−
                                <p class="interview-response">
 
−
                                    <em>A catheter like that would be useful, but it depends on how long your catheter would work for. Patients can have catheters fitted for 3 months or longer. Catheters are also removed for other reasons, for example if they become blocked. Blockage is particularly an issue when the patient is suffering from a urinary infection.</em>
 
−
                                </p>
 
−
                                <p>
 
−
                                    Laura’s response regarding how long a catheter remains in place spurred us into researching how we could keep our Solution bacteria alive. This is what we found.
 
−
                                </p>
 
 
                             </div>
 
                             </div>
−
                        </div>
+
                            </p>
−
                        <div id="practices-medical-professional-biochem-dept">
+
                             <h4>Compartmental diffusion; Gillespie</h4>
−
                             <h3>Biochemistry Department</h3>
+
 
                             <p>
 
                             <p>
−
                                 To gain a further insight into the feasibility of Solution, we gave two talks during the summer, one at the termly Corpus Christi College Biochemistry talks and another to a group of alumni from the Oxford Biochemistry department. Two important questions arose from these talks:
+
                                 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>
 
                             </p>
−
                            <ul>
 
−
                                <li>Have you considered whether the proteins you planning on secreting are immunogenic?</li>
 
−
                                <li>If you are to kill all of the pathogenic bacteria in the urinary tract, will that make fungal infections more likely?</li>
 
−
                            </ul>
 
 
                         </div>
 
                         </div>
 
                     </div>
 
                     </div>
 
                 </div>
 
                 </div>
 
                 <div class="section-spacer"></div>
 
                 <div class="section-spacer"></div>
−
                 <div class="section" id="practices-patients">
+
                 <div class="section" id="modeling-tutorial">
 
                     <div class="slim">
 
                     <div class="slim">
−
                         <h2>Talking to Patients</h2>
+
                         <h2>Modeling Tutorial</h2>
−
                         <div class="quote quote-right">
+
                    </div>
 +
                    <div id="modeling-tutorial-deterministic">
 +
                         <div class="slim">
 +
                            <h3>Deterministic</h3>
 
                             <p>
 
                             <p>
−
                                 the process is not very different to using antibiotics
+
                                 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>
 
                             </p>
 
                         </div>
 
                         </div>
−
                        <p>
 
−
                            During the process of talking to medical professionals, the importance of talking to patients as well was clear to us. On one occasion we spoke with Mavis, a patient staying on the Bedford Ward at the John Radcliffe Hospital. Her experience of urinary infections was that she had only suffered from them before having a catheter; since having a catheter fitted she hadn’t had any problems of infection. Her catheter was in place for up to 10 weeks. This enforced the importance of being able to keep our bacteria alive for a sustained period of time. When we asked her about treating infection with bacteria she said she would be happy to if it had been recommended to her by a doctor and told us that it is not dissimilar to using antibiotics.
 
−
                        </p>
 
 
                     </div>
 
                     </div>
−
                </div>
+
                    <div id="modeling-tutorial-stochastic">
−
                <div class="section-spacer"></div>
+
                        <div class="slim">
−
                <div class="section" id="practices-further-feedback">
+
                            <h3>Stochastic</h3>
−
                    <div class="slim">
+
                            <p>
−
                        <h2>Further Feedback</h2>
+
                                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>
+
                            </p>
−
                            We are still designing a delivery system for our UTI treatment. We plan on asking the public to select from a variety of options, and then to design our treatment around the most popular one.
+
                                \[\dfrac{dx}{dt}=c+r\]
−
                        </p>
+
                            <div class="image image-right">
−
                    </div>
+
                                <img src="https://static.igem.org/mediawiki/2015/1/1f/OxiGEM_Stochastic_Tutorial_01.jpg"/>
−
                </div>
+
                            </div>
−
                <div class="section-spacer"></div>
+
                            <p>
−
                <div class="section" id="practices-outreach">
+
                                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.
−
                    <div class="slim">
+
                            </p>
−
                        <h2>Outreach</h2>
+
                            <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:
−
                             To promote Synthetic Biology and iGEM, we’ve used a variety of approaches.
+
                            </p>
−
                        </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>
−
                     <div id="practices-outreach-uniq-workshop">
+
                     <div id="modeling-tutorial-finite-difference">
 
                         <div class="slim">
 
                         <div class="slim">
−
                             <h3>UNIQ Workshop</h3>
+
                             <h3>Finite Difference Models</h3>
 
                             <p>
 
                             <p>
−
                                 We met with 40 prospective Oxford students to teach them about Synthetic Biology. The students had in interest in Biochemistry but knew nothing about iGEM. We hammered home the key message of Synthetic Biology - that we achieve more progress by expanding a registry of standardised biological parts - through a 15 minute introductory presentation on BioBricks. We then split them into groups and gave each one a mentor from our iGEM team. We worked through questions to test their understanding in a tutorial style and asked them to explain the constructs of previous iGEM teams. They finished by presenting their findings to each other.
+
                                 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>
 
                             </p>
 
                         </div>
 
                         </div>
 
                     </div>
 
                     </div>
 
                 </div>
 
                 </div>
−
                 <div class="section-spacer">
+
                 <div class="section" id="references">
−
                </div>
+
                    <h2>References</h2>
−
                <div class="section" id="practices-videos">
+
                     <ol class="references">
−
                     <div class="slim">
+
                         <li>Synthetic Biology - A Primer, Paul S. Freemont, Richard I Kenne et al., 2012 Imperial College Press<a name="Ref1"></a></li>
−
                         <h2>Videos</h2>
+
 
−
                        <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>
−
                            We’ve produced the following videos to promote our project and help future teams:
+
 
−
                        </p>
+
                        <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>
−
                    </div>
+
 
−
                         <ul>
+
                        <li>Brian Ingalls,Mathematical Modeling in Systems Biology: An Introduction, 2012 University of Waterloo<a name="Ref4"></a></li>
−
                            <li>
+
 
−
                                Introduction to Oxford iGEM (<a href="https://www.youtube.com/watch?v=JQ9b_7fWrfk"> link</a>)<br>
+
                        <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>
−
                                <iframe width="560" height="315" src="https://www.youtube.com/embed/JQ9b_7fWrfk" frameborder="0" allowfullscreen></iframe>
+
 
−
                            </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>
+
 
−
                                Stochastic Modelling Tutorial (<a href="https://www.youtube.com/watch?v=xMy8oqP5Mis"> link</a>)<br>
+
                         <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>
−
                                <iframe src="https://www.youtube.com/watch?v=xMy8oqP5Mis"
+
 
−
                                width="560" height="315" frameborder="0" allowfullscreen>
+
                        <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>
−
                                </iframe>
+
                     </ol>
−
                            </li>
+
−
                         </ul>
+
−
                </div>
+
−
                <div class="section references" id="practices-references">
+
−
                    <div class="slim">
+
−
                        <h2>References</h2>
+
−
                         <ol>
+
−
                            <li>
+
−
                                [1] You can see all of our responses <a href="https://docs.google.com/document/d/1MAzVkPVVGAbAbFMIN9OIvZtC4vBpGkW2et9W3omnkzM/edit?usp=sharing">here</a>.
+
−
                            </li>
+
−
                        </ol>
+
−
                     </div>
+
 
                 </div>
 
                 </div>
 
             </div>
 
             </div>
Line 214: Line 392:
 
                 <ul id="sidebar" class="nav nav-stacked affix-top sm-hidden xs-hidden" data-spy="affix">
 
                 <ul id="sidebar" class="nav nav-stacked affix-top sm-hidden xs-hidden" data-spy="affix">
 
                     <li>
 
                     <li>
−
                         <a href="#practices-introduction">Introduction</a>
+
                         <a href="#introduction">Introduction</a>
 
                     </li>
 
                     </li>
 
                     <li>
 
                     <li>
−
                         <a href="#practices-public-interaction">Public Interaction</a>
+
                         <a href="#gene-expression">Gene Expression</a>
 
                         <ul class="nav nav-stacked">
 
                         <ul class="nav nav-stacked">
−
                             <li><a href="#practices-public-interaction-project-choice">Project Choice</a></li>
+
                             <li><a href="#gene-expression-deterministic">Deterministic</a></li>
−
                             <li><a href="#practices-public-interaction-initial-feedback">Initial Feedback</a></li>
+
                             <li><a href="#gene-expression-stochastic">Stochastic</a></li>
 
                         </ul>
 
                         </ul>
 
                     </li>
 
                     </li>
 
                     <li>
 
                     <li>
−
                         <a href="#practices-medical-professionals">Talking to Medical Professionals</a>
+
                         <a href="#diffusion">Diffusion</a>
 
                         <ul class="nav nav-stacked">
 
                         <ul class="nav nav-stacked">
−
                             <li><a href="#practices-medical-professionals-churchill">Churchill Hospital, Oxford</a></li>
+
                             <li><a href="#diffusion-deterministic">Deterministic</a></li>
−
                             <li>
+
                             <li><a href="#diffusion-stochastic">Stochastic</a></li>
−
                                <a href="#practices-medical-professionals-jr">John Radcliffe Hospital, Oxford</a>
+
−
                                <ul class="nav nav-stacked">
+
−
                                    <li><a href="#practices-medical-professionals-jr-laura-1">First Interview with Laura Evans, Adams Ward</a></li>
+
−
                                </ul>
+
−
                            </li>
+
−
                            <li><a href="#practices-medical-professional-biochem-dept">Biochemistry Department</a></li>
+
 
                         </ul>
 
                         </ul>
 
                     </li>
 
                     </li>
 
                     <li>
 
                     <li>
−
                         <a href="#practices-further-feedback">Further Feedback</a>
+
                         <a href="#modeling-tutorial">Modeling Tutorial</a>
−
                    </li>
+
−
                    <li>
+
−
                        <a href="#practices-outreach">Outreach</a>
+
 
                         <ul class="nav nav-stacked">
 
                         <ul class="nav nav-stacked">
−
                             <li><a href="#practices-outreach-uniq-workshop">UNIQ Workshop</a></li>
+
                             <li><a href="#modeling-tutorial-deterministic">Deterministic</a></li>
 +
                            <li><a href="#modeling-tutorial-stochastic">Stochastic</a></li>
 +
                            <li><a href="#modeling-tutorial-finite-difference">Finite Difference Models</a></li>
 
                         </ul>
 
                         </ul>
 
                     </li>
 
                     </li>
 
                     <li>
 
                     <li>
−
                         <a href="#practices-videos">Videos</a>
+
                         <a href="#references">References</a>
−
                    </li>
+
−
                    <li>
+
−
                        <a href="#practices-references">References</a>
+
 
                     </li>
 
                     </li>
 
                 </ul>
 
                 </ul>

Revision as of 13:36, 15 August 2015

Modeling

Introduction

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.

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.

Gene Expression

Deterministic

A deterministic model evolves according to ordinary or partial differential equations. We model the following 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\]

Using the following deterministic equations [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

Symbol Definition Initial Value/Literature Value Fitted
\([Arab]\) The concentration of Arabinose \(1\times10^{-4}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] ?

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.

Stochastic

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.[2, 4].

The Gillespie Algorithm

  1. Find random numbers \(r_{1}\) and \(r_{2}\) which are uniformly distributed between 0 and 1.
  2. Find propensities \(\alpha_{i}\) of \(N\) possible reactions. Find \(\alpha_{0}=\sum\limits_{i=1}^N \alpha_{i}\).
  3. Compute time \(\tau=\dfrac{1}{\alpha_{0}}\ln\left(\dfrac{1}{r_{1}}\right)\).
  4. 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}\).
  5. Implement reaction \(j\) and increment time \(t\) to \(t+\tau\).
  6. Repeat steps 1-5 until a sufficient time or computation budget has been exhausted.

Turning the reaction rates in our model into propensities is simple - for example, the propensity for the association reaction \(Arab+AraC\rightarrow Arab:AraC\) is \(\alpha_{3}\left[Arab\right]\left[AraC\right]\).

Diffusion

Deterministic

The standard from for the diffusion equation,

\[\dfrac{d[P]}{dt}=D\dfrac{d^{2}[P]}{dx^{2}}\]

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

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 tutorial section.

Stochastic

We have the option of using three stochastic models for diffusion - random walk, random velocity and compartmental Gillespie.

Random Walk

We initially take a collection of particles and, for each one, compute the number

\[dr=\sqrt{2Ddt}\xi_{r}\]

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\)].

Random Velocity

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.

Random Velocity Algorithm

  1. For each particle, generate a random speed \(s\) in a random direction.
  2. Compute the position of the particle after a timestep \(\Delta t\).
  3. Generate a random number \(r\) uniformly distributed between 0 and 1.
  4. If \(r<\dfrac{s^{2}}{2D}\Delta t\), calculate a new random direction of velocity for the particle.
  5. Repeat steps 2-5 for each particle until a certain time or computational budget is exhausted.

Compartmental diffusion; Gillespie

We can divide our region into a system of compartments

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

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

\[k=\dfrac{D}{(dx)^{2}}\]

where \(dx\) is the physical size of each compartment [2]. We now run the Gillespie algorithm but with \(2(N-1)\) reactions to consider.

Comparison of methods

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.

The challenge

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.

Modeling Tutorial

Deterministic

Our goal is to model a system of reactions such as the reversible reaction

\[A \mathrel{\mathop{\rightleftharpoons}^{\mathrm{k_{+}}}_{\mathrm{k_{-}}}} B\]

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

\[\dfrac{dB}{dt}=k_{+}A\]

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

\[\dfrac{dA}{dt}=k_{-}B-k_{+}A\] \[\dfrac{dB}{dt}=k_{+}A-k_{-}B\]

At this point, you can type these into your favourite computing program and solve them.

Stochastic

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:

\[\dfrac{dx}{dt}=c+r\]

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 do better than this when dealing with systems such as our example in the previous section. We use the Gillespie algorithm to find:

  1. The time taken between reactions
  2. The reaction that took place in that time.

as neither of these could be determined by our initial method [2, 4].

1) 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(t)=e^{-\alpha_{0}t}\]

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\).

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.

2) 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.

Finite Difference Models

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.

The 1D formula to iterate over for the finite difference model of diffusion is

\[P_{j}^{n+1}=(1-2r)P_{j}^{n}+r(P_{j-1}^{n}+P_{j+1}^{n})\]

The 2D formula to iterate over is

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

Similarly, for 3D:

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

we (arbitrarily) choose our boundary conditions to be continuous at the boundary.

References

  1. Synthetic Biology - A Primer, Paul S. Freemont, Richard I Kenne et al., 2012 Imperial College Press
  2. R. Erban, S. J. Chapman, and P. Maini. A practical guide to stochastic simulations of reaction- diffusion processes, 2007. Available here
  3. 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
  4. Brian Ingalls,Mathematical Modeling in Systems Biology: An Introduction, 2012 University of Waterloo
  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. 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
  8. 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