Difference between revisions of "Team:UT-Tokyo/Modeling"
(Prototype team page) |
|||
(7 intermediate revisions by 3 users not shown) | |||
Line 1: | Line 1: | ||
− | {{UT-Tokyo}} | + | {{UT-Tokyo/Header}} |
+ | |||
<html> | <html> | ||
+ | <head> | ||
+ | <meta name='viewport' content='width=device-width, initial-scale=1, user-scalable=no'> | ||
+ | <meta http-equiv='Content-Type' content='text/html; charset=UTF-8'> | ||
+ | <title>The Pattern-Formation Game</title> | ||
+ | <link id='css' rel='stylesheet' href='./css/reset.css?v=2' /> | ||
+ | <link id='css' rel='stylesheet' href='./css/style.css?v=2' /> | ||
+ | <script type="text/javascript" src='./js/include.js'></script> | ||
+ | <script type="text/javascript" | ||
+ | src="http://cdn.mathjax.org/mathjax/latest/MathJax.js?config=TeX-AMS_HTML"> | ||
+ | </script> | ||
+ | <script type="text/javascript" src="./js/jaxconfig.js"></script> | ||
+ | </head> | ||
+ | <body> | ||
+ | <script type="text/javascript"> | ||
+ | include("./header.html"); | ||
+ | </script> | ||
− | <h2> | + | <div class='main'> |
+ | <div class='panel' id='top'> | ||
+ | <h2>TOP</h2> | ||
+ | <h3>Overview</h3> <br> | ||
+ | <p>To generate Turing pattern using synthetic biology, we thought quite a few systems. We analyzed some of them and checked if they can generate a pattern and how the pattern can be. We used some differential equations and multi-agent simulation for this. From modeling, it was suggested that it is difficult to find the parameters that certainly satisfies the condition on pattern formation. However, some trends on parameters were observed. </p> | ||
− | < | + | <h3>The aim of the modeling</h3> <br> |
− | < | + | |
− | <p>In order to | + | <p> Turing mechanism is based on a mathematical model (<a href = " " name="mathematics">mathematics</a> page). We modeled our system in order to understand it mathematically as well as to find the parameter that satisfies the condition for forming a pattern.</p> |
+ | |||
+ | <h3>Background</h3> | ||
+ | |||
+ | <p>In order to understand our system and analysis, some basic knowledge for Turing pattern is needed. We’d like to introduce some mathematical basics of Turing Pattern to make it easier to understand our project in the <a id = "mathematics">mathematics</a> page. </p> | ||
+ | |||
+ | <h3>The flow of whole modeling project</h3> <br> | ||
+ | |||
+ | <h4>Analysis and explanation of pattern generation</h4> <br> | ||
+ | |||
+ | <p>We first showed that the first system cannot form a pattern. To solve this problem, we decided to use colicin. Why we choosed to use colicin is qualitatively explained here.</p> | ||
+ | <h4>System Analysis</h4> | ||
+ | <p>As we saw above, it is difficult to satisfy the conditions for forming Turing pattern. Thus, at first we tried to weaken this conditions, to find the paramerer range which would form the pattern. <br> | ||
+ | We succeeded in determining the range, however, due to the fructuations, it turnd out that they no longer certainly generate Turnig pattern.</p> | ||
+ | <h4>How the pattern will be in the last</h4> | ||
+ | |||
+ | <p>In this part, we simulated how the pattern wil be. we constructed a mathematical model which describe our system, and applied the sets of parameters found in the previous analysis.<br> | ||
+ | we did simulation on on dimension and two dimension case.</p> | ||
+ | |||
+ | <p>one dimension: <br> | ||
+ | The analysis on one dimension. Assuming that <i><i>E. coli</i></i> are in a long and thin container. Stripe pattern would appear, if the optimization of part 1 is well conducted. We also checked how long it would take to form the pattern, and whether a periodical pattern would appear. <br> | ||
+ | <br> | ||
+ | two dimension: <br> | ||
+ | Generally, the model showed various patterns: spot, stripe, competition,... etc. in two dimension case. <br> | ||
+ | We used the set of parameters obtained in the first analysis. <br> | ||
+ | in this model, we assumed the container large enough to ignore effescs from the boundary: i.e. the diffusion equation holds for this system. <br> | ||
+ | The simulation under this assumption ended only in activator spot pattern.</p> | ||
+ | |||
+ | |||
+ | <h3>New model </h3> | ||
+ | |||
+ | <p>We also conducted the analysis above to a new model. It is inspired by zebrafish and includes the effect of colicin.</p> | ||
+ | |||
+ | CODE | ||
+ | see this GitHub page. We modeled first and second sections by MATLAB, and the last one by Java. | ||
+ | |||
+ | <p>[1] Instabilities in spatially extended predator–prey systems Spatio-temporal patterns in the neighborhood of Turing–Hopf bifurcations, Baurmann et. al 2006 Journal of Theoretial Biology 245 (2007) 220–229</p> | ||
</div> | </div> | ||
− | <p> | + | <div class='panel' id='system'> |
+ | <h2>System</h2> | ||
+ | <p>First, we will explain why the first model is not suitable.</p> | ||
+ | <h3>The reason why the first model cannnot seem to form pattern</h3> | ||
+ | <p>The interactions between <i><i>E. coli</i></i> and AHL are expressed as two differential equations below.<br> | ||
+ | $$ | ||
+ | \frac{du}{dt}=d_u\nabla^2u+u(1-u)-{\frac{a}{a+K}}{\delta}u=d_u\nabla^2u+f(u,a)\\ | ||
+ | \frac{da}{dt}=d_a\nabla^2a+u-{\gamma}a=d_a\nabla^2a+g(u,a) | ||
+ | $$ | ||
+ | u: concentration of <i><i>E. coli</i></i><br> | ||
+ | a: concentration of AHL<br> | ||
+ | d_u: the diffusion coefficient of u.<br> | ||
+ | d_a: the diffusion coefficient of a.<br> | ||
+ | {\gamma}: degradation coefficient of AHL<br> | ||
+ | (Note that the parameters are non-dimensionized.)<br><br> | ||
+ | </p> | ||
+ | <h3>equation 1</h3><br> | ||
<p> | <p> | ||
− | Here are a | + | The concentration of <i><i>E. coli</i></i> increases by multiplication in proportion to its concentration, but if the concentration is too high, the growth rate decreases because of the lack of nutrition(+u(1-u)). If <i><i>E. coli</i></i> receives enough AHL, it is killed by Barnase {\frac{a}{a+K}}{\delta}u.The reaction term {\frac{a}{a+K}}{\delta} can be approximated by Hill equation because the system of <i><i>E. coli</i></i> death is similar to ligand-receptor system of protein.</p> |
+ | |||
+ | <h3>equation2</h3><br> | ||
+ | <p>The concentration of AHL increases in proportion to that of <i><i>E. coli</i></i>, and AHL degrades in proportion to -{\gamma}.<br> | ||
+ | We concluded that this model may not be able to form a pattern because this model doesn’t satisfy required mathematical conditions of Turing pattern.</p> | ||
+ | |||
+ | <h3>Proof</h3> | ||
+ | <p>Following two equations should be satisfied.<br> | ||
+ | |||
+ | \begin{equation} | ||
+ | f_u+g_a\lt 0 | ||
+ | \label{eq:vern1} | ||
+ | \end{equation} | ||
+ | \begin{equation} | ||
+ | f_ug_a-f_vg_u\gt 0 | ||
+ | \label{eq:vern2} | ||
+ | \end{equation} | ||
+ | |||
+ | First, there are two possible equilibrium points \(\frac{du}{dt}= \frac{da}{dt}= 0\).<br> | ||
+ | \((u,a)=(0,0)\) and \((u,a)=(u_0,a_0)\) <br> | ||
+ | \(u_0\) and \(a_0\) satisfies the relations,\(u_0 = \gamma a_0\) and \(f(a_0)=1− u_0\)<br> | ||
+ | \((u,a)=(0,0)\) doesn’t satisfy equation [\ref{eq:vern1}] and [\ref{eq:vern2} ].<br> | ||
+ | When \((u,a)=(u_0,a_0)\), <br> | ||
+ | $$ | ||
+ | f_u=1 – 2u_0 – {\frac{a}{a+K}}{\delta}=1 – 2u_0 – (1 – u_0)= –u_0 \lt 0 | ||
+ | $$ | ||
+ | $$ | ||
+ | f_a= –u_0・d_u{\frac{f(a)}{da}}|a_0 \lt 0 \qquad (\mathrm{for} \ d_u\frac{f(a)}{d_a} \gt 0)$$ | ||
+ | because \(f(a)\) is approximated by Hill’s equation.<br> | ||
+ | |||
+ | \(g_u= 1\)<br> | ||
+ | \(g_a= –{\gamma}\lt0\)<br> | ||
+ | Putting these values into equation \ref{eq:vern1} and \ref{eq:vern2} ,<br> | ||
+ | \(f_u+g_a= – u_0 –\gamma\lt 0\)<br> | ||
+ | \(f_ug_a−f_ag_u = u_0(\gamma+ d_u\frac{f(a)}{d_a}|a_0)\gt0\)<br> | ||
+ | So, \((u,a)=(u_0,a_0)\) satisfies those requirements.<br> | ||
+ | To form Turing pattern, however, there is still one more condition to be satisfied two equations below.<br> | ||
+ | \((d_ug_a + f_u)^2 – 4d_u(f_ug_a−f_ag_u) \gt 0\)<br> | ||
+ | \(d_ug_a + d_af_u\gt 0\)<br> | ||
+ | For \(g_a\lt 0, f_u\lt 0\), equation 6 is never satisfied whatever value d_u takes.<br> | ||
+ | Requirement 2 unsatisfied, this model cannot form Turing pattern.<br> | ||
+ | To meet the equation 6, \(f_u\) must be bigger than zero. And to satisfy this, the positive feedback <br> | ||
+ | Therefore, there must be a system that has allee effect, or density dependant effect. The former means that "When <i>E. coli</i> increases, the growth rate of <i>E. coli</i> also increases."Allee effect is a positive correlation between population size and indivisual fitness. Mimura and Murray(1972) said that Allee effect is important for some pattern formation on predator-prey systems[1]. However, there is no Allee effect in <i>E. coli</i>. <i>E. coli</i> increase however little amount they are(one colony would be formed from one cell). The multiplication rate of <i><i>E. coli</i></i> is small when it’s introduced into culture medium, for it’s not adapted to the environment. But that doesn’t affect the formation of Turing pattern, because the equilibrium point of Turing pattern occurs on log phase or stationary phase.The latter is saying that "When there are very much <i>E. coli</i>, the effect of AHL is not proportional to <i>E. coli</i>'s density"; the effect will be lower. This is seen in the model of "Prey-Predator model". Since predator eats no more than they can eat, when preys are rich, the probablity of being killed of an indivisual prey would be decreased. However, since AHL is much smaller and richer than <i>E. coli</i>, there seems little chance to expect it to happen.<br> | ||
+ | It seemed that both cannot be satisfied only by one type of <i>E. coli</i>. Therefore, we decided to solve this problem by using two types of <i>E. coli</i> and made positive-feedback-like system. Most positive feedback system have bistability (which does not exist in the system of population of <i>E. coli</i>). In other words, the space would be separated by two areas, according to the convergence value. The convergence value is polarized. This is a kind of the system of the competition exclusion. To learn more about the principle of competition exclusion, see Murray, mathematical biology volume 1.<br> | ||
+ | We wanted to realize competition exclusion system. We chose the system that utilizes colicin. See the system page for detail. | ||
+ | Why colicin? Because colicin does not diffuse less quickly than AHL. Short-range positive feedback and long-range negative feedback are necessary for Turing pattern. Hence, bistability should be achieved locally. <br> | ||
+ | Chao and Levin says that this colicin-using system (no AHL effect on our system) exhibit the bistability[3]. We show this in the <a href = "mathematics" >mathematics</a> page.<br> | ||
+ | We considered two cases.<br> | ||
+ | Case 1. <i>E. coli</i> that emit AHL(shown u in figure) and <i>E. coli</i> that emit colicin(shown v in figure) are different<br> | ||
+ | formula<br> | ||
+ | $$ | ||
+ | \frac{du}{dt} = u({\alpha} - (u + v) - {\kappa}c - \frac{a}{a + K}{\delta}) + du{\nabla}^2u\\ | ||
+ | \frac{dv}{dt} = v({\beta} - (u + v)\) + du{\nabla}^2v \\ | ||
+ | \frac{da}{dt} = u - {\gamma}a + da{\nabla}^2a \\ | ||
+ | \frac{dc}{dt} = v - {\mu}c + dc{\nabla}^2c \\ | ||
+ | $$ | ||
+ | Here, AHL, E. coli, colicin are all scaled, and unnecessary parameters are deleted.<br> | ||
+ | Case 2. E. coli that emit AHL(shown u in figure) and E. coli that is killed by colicin(shown v in figure) are same.<br> | ||
+ | formula<br> | ||
+ | $$ | ||
+ | \frac{du}{dt} = u({\alpha} - (u + v) - \frac{a}{a + K}{\delta}) + du{\nabla}^2u \\ | ||
+ | \frac{dv}{dt} = v(\beta - (u + v) - {\kappa} c) + du{\nabla}^2v \\ | ||
+ | \frac{da}{dt} = u - {\gamma}a + da{\nabla}^2a \\ | ||
+ | \frac{dc}{dt} = u - {\mu}c + dc{\nabla}^2c \\ | ||
+ | $$ | ||
+ | We first investigated "which parameter set can lead to form pattern" by stability analysis.<br> | ||
+ | First, we defined parameter range. See the "parameter"<link> page to see the detail.<br> | ||
+ | Second, we varied parameters according to parameter range, and constructed some parameter sets.<br> | ||
+ | Third, for each parameter set, we saw if it satisfies the condition of Turing pattern, mentioned here (<a id="mathematics">mathematics</a>). If it is satisfied, it is registered as a "passed parameter set". Each parameter set is also characterized by its largest eigenvalue of stability matrix, with and without diffusion term, sensitivity, and the wavelength that takes the maximum eigenvalue.<br> | ||
+ | Then, after seeking parameters, we saw some trend of parameters. We also ranked parameter sets, acording to their stability matrixes' eigenvalue. Since <i>E. coli</i> increases in the order of hour, it is considered to take a lot of time for pattern formation. However, if it takes much time for pattern formation, there may arise some difficulty on the experiment that the nutrients for <i><i>E. coli</i></i> would be run out until the pattern is formed.<br> | ||
+ | There are 10 parameters in the last formulas. However, there were some parameters that do not seem to have been found. There were also some parameters that can be changed by experiment; strength of the constitutive promoter, a copy number of the promoter, etc... For those parameters, we defined a parameter range so that each parameter can affect forming a pattern. See parameter page for detail.</p> | ||
+ | <h4>Result</h4> | ||
+ | <p> two cases on the same parameter range, it turned out that more patterns were formed in Case 1 than Case 2. And the biggest eigenvalue of Case 1 was smaller than that of Case 2.</p> | ||
+ | |||
+ | |||
+ | <h3>The trend of parameters</h3> | ||
+ | <p>Case 1<br> | ||
+ | ・Set same parameter range, Case 1 is less likely to make pattern than Case 2.<br> | ||
+ | ・alpha/beta > 1 should be satisfied. In case 2, however, alhpa should not be very large compared to beta. <br> | ||
+ | ・K should not be so large, nor so small. <br> | ||
+ | ・delta should be larger than 1, so that it can kill u or v. This means that <br> | ||
+ | ・kappa should be large, especially if you want to make pattern fast (from the relationship between the eigenvalue).<br> | ||
+ | ・there seems small sensitivity in gamma.<br> | ||
+ | ・mu should be large, so that the order is near alpha and beta. mu should not be so large(larger than 0.7).<br> | ||
+ | ・du, dc da have no trend as long as its range is narrow enough<br> | ||
+ | ・i think there would be a relationship between K and delta. I will check it.<br> | ||
+ | ・ delta should be large and K should be small.<br> | ||
+ | ・ K and delta is related. If K is large, delta must not be so small.<br> | ||
+ | <img src="https://static.igem.org/mediawiki/2015/4/4b/UT_TOKYO_Modeling_alphabeta1.jpg" ></br> | ||
+ | x-axis alpha, y-axis beta for case 2.</p> | ||
+ | |||
+ | |||
+ | |||
+ | <h3>Sensitivity and Robustness</h3> | ||
+ | <p>We've tested whether pattern forms or not if we change the parameters, and what parameters would exert effects on the formation. Here are the details:<br> | ||
+ | There are two points of view in sensitivity; local or global. This classification was suggested by iGEM Waterloo 2014 team. Local sensitivity is a test about how the result turns by changing only one parameter. Global one sees the same thing by changing various parameters. We've done the two tests and examined what parameter sets would satisfy the conditions of Turing Pattern.<br></p> | ||
+ | |||
+ | <h3>case 2</h3> | ||
+ | <p>We observed strong polarization along the eigenvalue gets bigger or smaller.</br> | ||
+ | Polarization of 50 of those, which have small eigenvalue was not so stronger than previous examination. The more eigenvalue get bigger, the more strong polarization occurs. In other words, it can be said that <i>something which makes eigenvalue bigger</i> exists in the place where alpha's max and beta's minimum. </br> | ||
+ | The max place where each parameters take would be the boundary line of whether pattern is formed or not. If the parameter go beyond the boundary line, parameter sets couldn't have stable point or could not be stable, because of no occurrence of diffusion. Thus, we believe there would be trade-off between size of eigenvalue and robustness.</br> | ||
+ | We analyzed the stability of the random numbers around the parameter (parameter*0.5~parameter*1.5) and examined if there was deviation of distribution when pattern formed. The distribution should be linearly condition if there was no deviation.</br> | ||
+ | |||
+ | |||
+ | Second, we simulated Turing pattern numerically.</br> | ||
+ | |||
+ | in two dimension, and saw if it forms a pattern.</br> | ||
+ | See this figure.</br> | ||
+ | <img src="https://static.igem.org/mediawiki/2015/8/8b/UT_TOKYO_Modeling_2D.png"></br> | ||
+ | Here, you can see pattern is formed. | ||
</p> | </p> | ||
− | < | + | </div> |
− | + | ||
− | + | ||
− | + | ||
+ | <div class='panel' id='system2'> | ||
+ | <h2>System2</h2> | ||
+ | <h3>Another system</h3> | ||
+ | <p>Here, we explain the "third" model. This is similar to the colicin-using model. There is some difference.</br> | ||
+ | |||
+ | The third strategy is inspired by the pattern formation mechanism of zebrafish[1]. </br> | ||
+ | Two types of <i><i>E. coli</i></i> (activator and inhibitor) react to each other in two manners according to the distance between them.</br> | ||
+ | </br> | ||
+ | When <i><i>E. coli</i></i> (activator) and <i><i>E. coli</i></i> (inhibitor) are close, they inhibit each other in the same manner as in the conventional system. This control loop functions as local activation.</br> | ||
+ | </br> | ||
+ | In addition, <i><i>E. coli</i></i> (activator) activates <i><i>E. coli</i></i> (inhibitor) through AHL in long distance (detailed scheme is explained in the chapter of System). The control loop composed of this activation and the short-range inhibition of <i><i>E. coli</i></i> (activator) by activated <i><i>E. coli</i></i> (inhibitor) functions as lateral inhibition. </br> | ||
+ | </br> | ||
+ | "Activator" <i><i>E. coli</i></i> promotes the multiplication of "Inhibitor" <i><i>E. coli</i></i>, and "Inhibitor" inhibits the multiplication of "Activator" near of "Inhibitor"</br> | ||
+ | </br> | ||
+ | |||
+ | First, we show that there are two cases for this modeling.</br> | ||
+ | </br> | ||
+ | |||
+ | Case 3: <i>E. coli</i> that emit AHL(shown u in figure) and <i>E. coli</i> that emit colicin(shown v in figure) are different.</br> | ||
+ | formula</br> | ||
+ | $$ | ||
+ | \frac{du}{dt} = u({\alpha} - (u + v) - {\kappa}c) + du{\nabla}^2u\\ | ||
+ | \frac{dv}{dt} = v({\beta} - (u + v)\) + \frac{a}{a + K}{\delta}) + du{\nabla}^2v \\ | ||
+ | \frac{da}{dt} = u - {\gamma}a + da{\nabla}^2a \\ | ||
+ | \frac{dc}{dt} = v - {\mu}c + dc{\nabla}^2c \\ | ||
+ | $$ | ||
+ | |||
+ | |||
+ | Case 4; <i>E. coli</i> that emit AHL(shown u in figure) and <i>E. coli</i> that emit colicin are same.</br> | ||
+ | formula</br> | ||
+ | $$ | ||
+ | \frac{du}{dt} = u({\alpha} - (u + v)) + du{\nabla}^2u \\ | ||
+ | \frac{dv}{dt} = v({\beta} - (u + v) - {\kappa}c + \frac{a}{a + K}{\delta}) + du{\nabla}^2v \\ | ||
+ | \frac{da}{dt} = u - {\gamma}a + da{\nabla}^2a \\ | ||
+ | \frac{dc}{dt} = u - {\mu}c + dc{\nabla}^2c \\ | ||
+ | $$ | ||
+ | |||
+ | The parameter range which is different from strategy 2 is delta. delta + beta must be lower than max value of alpha or beta, because barstar or colI can't restore the growth rate of <i>E. coli</i> more than its max value. Hence, the sum of ${\delta}$ and ${\beta}$ must be the same or below the maximum growth rate.</br> | ||
+ | The method is not different from the first one.</br> | ||
+ | </br> | ||
+ | [2]Nakamasu, A., Takahashi, G., Kanbe, A., & Kondo, S. (2009). Interactions between zebrafish pigment cells responsible for the generation of Turing patterns. Proceedings of the National Academy of Sciences, 106(21), 8429-8434. | ||
+ | </p> | ||
</div> | </div> | ||
+ | |||
+ | <div class='panel' id='mathematics'> | ||
+ | <p> | ||
+ | In order to understand our system and analysis, some basic knowledge for Turing pattern is needed. We’d like to introduce some mathematical basics of Turing Pattern to make it easier to understand our project. | ||
+ | </p> | ||
+ | <h3>1.Stability</h3> | ||
+ | <p> The idea of “stability” is so important in our modeling.Consider conditions of | ||
+ | $ \frac{du}{dt} = f(u, v), f(u_0, v_0)=0\\ | ||
+ | \frac{dv}{dt} = g(u, v), g(u_0, v_0)=0 $ | ||
+ | where u is the concentration of activator, v is of inhibitor. So both u and v are functions of time(t) and location. | ||
+ | </p> | ||
+ | <p> | ||
+ | Suppose $(u, v) = (u_0, v_0)$ and a fluctuation then occurs to the point. | ||
+ | Let $(u_0 + u_s, v_0 + v_s)$ be the point after the fluctuation occurred. | ||
+ | Here we can write the differential equations as | ||
+ | |||
+ | $ \frac{du_s}{dt} = f(u_0 + u_s, v_0 + v_s) \\ | ||
+ | \frac{dv_s}{dt} = g(u_0 + u_s, v_0 + v_s) $ | ||
+ | <p> | ||
+ | Given $u_s$, $v_s$ are small, we can expand $f$,$g$ around $(u_0, v_0)$ excluding the terms whose degree is over one. | ||
+ | </p> | ||
+ | <p> | ||
+ | That is | ||
+ | |||
+ | $ \frac{du_s}{dt} = f(u_0 + u_s, v_0 + v_s) = f_uu_s + f_vv_s\\ | ||
+ | \frac{dv_s}{dt} = g(u_0 + u_s, v_0 + v_s) = g_uu_s + g_vv_s $ | ||
+ | where $ f_u = f_u(u_0, v_0),\\ | ||
+ | f_v = f_v(u_0, v_0)\\ | ||
+ | g_u = g_u(u_0, v_0),\\ | ||
+ | g_v = g_v(u_0, v_0) $ | ||
+ | |||
+ | (From now on, we mention $f_u$, $f_v$, $g_u$, $g_v$ as $f_u(u_0, v_0)$, $f_v(u_0, v_0)$, $g_u(u_0, v_0)$, $g_v(u_0, v_0)$ until the end for convenience) | ||
+ | </p> | ||
+ | <p> | ||
+ | In another way, | ||
+ | \begin{eqnarray} | ||
+ | \left( | ||
+ | \begin{array}{cccc} | ||
+ | \frac{du_s}{dt}\\ | ||
+ | \frac{dv_s}{dt} \end{array} | ||
+ | \right) | ||
+ | &=& | ||
+ | \left( | ||
+ | \begin{array}{cc} | ||
+ | {f_u} & {f_v}\\ | ||
+ | {g_u} & {g_v} | ||
+ | \end{array} | ||
+ | \right) | ||
+ | \left( | ||
+ | \begin{array}{cc} | ||
+ | {u_s}\\ | ||
+ | {v_s} | ||
+ | \end{array} | ||
+ | \right) | ||
+ | ...(A) | ||
+ | \end{eqnarray} | ||
+ | </p> | ||
+ | <p> | ||
+ | We want to find the solutions of $(u_s, v_s)$. | ||
+ | </p> | ||
+ | <p> | ||
+ | Given both $u_s$ and $v_s$ are in proportion to $e^{{\lambda}t}$, it can be deduced that ${\lambda}$ is equal to the eigenvalues of the matrix in (A). And as the matrix in (A) has two eigenvalues, it can be said that $u_s$ and $v_s$ are described as a linear combination of $e^{{\lambda}t}s$. If all of ${\lambda}$s are negative, $e^{{\lambda}t}$s get smaller and smaller as time passes. It means every time when (u,v) feels any small fluctuations, (u,v) returns to the point $(u_0, v_0)$. In other words, if all ${\lambda}$s are negative, $(u_0, v_0)$ is a “stable” point. On the contrary, if any of ${\lambda}$s is positive, as $e^{{\lambda}t}$s whose ${\lambda}$ is positive gets larger and larger, it can be said that (u, v) won’t return to $(u_0, v_0)$. Then $(u_0, v_0)$ isn’t a “stable“ point. Conditions for a point $(u_0, v_0)$ to be stable are known as | ||
+ | |||
+ | $ f_u + g_v < 0\\ | ||
+ | f_ug_v - f_vg_u > 0 $ | ||
+ | These conditions can be expanded to cases of three and more factors. | ||
+ | </p> | ||
+ | |||
+ | <h3>2. Turing Pattern</h3> | ||
+ | <p> | ||
+ | For Turing Pattern, such points are required that satisfy the conditions of… | ||
+ | <ul> | ||
+ | <li>・stable when dilution isn’t taken into account </li> | ||
+ | <li>・not stable when dilution is taken into account </li> | ||
+ | </ul> | ||
+ | When you consider dilution, the equations change into | ||
+ | |||
+ | $ \frac{du}{dt} = f(u, v) + D_u{\nabla}^2u, f(u_0, v_0) = 0\\ | ||
+ | \frac{dv}{dt} = g(u, v) + D_v{\nabla}^2v, g(u_0, v_0) = 0 $ | ||
+ | |||
+ | these conditions can be changed into | ||
+ | $ \frac{du_s}{dt} = f_uu_s + f_vv_s + D_u{\nabla}^2u_s\\ | ||
+ | \frac{dv_s}{dt} = g_uu_s + g_vv_s + D_v{\nabla}^2v_s $ | ||
+ | |||
+ | In the same way as 1.Stability, given u,v are proportional to $e^{ikx}$, it can be deduced that | ||
+ | $ \frac{du_s}{dt} = (f_u - D_uk^2)u_s + f_vv_s\\ | ||
+ | \frac{dv_s}{dt} = g_uu_s + (g_v - D_vk^2)v_s $ | ||
+ | |||
+ | that is | ||
+ | |||
+ | \begin{eqnarray} | ||
+ | \left( | ||
+ | \begin{array}{cc} | ||
+ | \frac{du_s}{dt}\\ | ||
+ | \frac{dv_s}{dt} | ||
+ | \end{array} | ||
+ | \right) | ||
+ | = | ||
+ | \left( | ||
+ | \begin{array}{cc} | ||
+ | {f_u - D_uk^2} & {f_v}\\ | ||
+ | {g_u} & {g_v - D_vk^2} | ||
+ | \end{array} | ||
+ | \right) | ||
+ | \left( | ||
+ | \begin{array}{cc} | ||
+ | {u_s}\\ | ||
+ | {v_s} | ||
+ | \end{array} | ||
+ | \right) | ||
+ | \end{eqnarray} | ||
+ | </p> | ||
+ | <p> | ||
+ | For the point of $(u_0, v_0)$ to be unstable, it’s necessary that | ||
+ | $ (D_ug_v + D_vf_u)^2 - 4D_uD_v(f_ug_v - f_vg_u)> 0 \\ | ||
+ | D_ug_v + D_vf_u > 0 $ | ||
+ | </p> | ||
+ | <p> | ||
+ | This means that the larger the ratio of $D_u$ and $D_v$ is, the more the point considered is likely to be unstable. To explain more detail, it's deduced that | ||
+ | |||
+ | \begin{eqnarray} | ||
+ | \left( | ||
+ | \begin{array}{cc} | ||
+ | {f_u} & {f_v}\\ | ||
+ | {g_u} & {g_v} | ||
+ | \end{array} | ||
+ | \right) | ||
+ | = | ||
+ | \left( | ||
+ | \begin{array}{cc} | ||
+ | {+} & {+}\\ | ||
+ | {-} & {-} | ||
+ | \end{array} | ||
+ | \right) | ||
+ | \ or | ||
+ | \left( | ||
+ | \begin{array}{cc} | ||
+ | {+} & {-}\\ | ||
+ | {+} & {-} | ||
+ | \end{array} | ||
+ | \right) | ||
+ | \end{eqnarray} | ||
+ | |||
+ | There're two possible situations as above,but we consider only the | ||
+ | \begin{eqnarray} | ||
+ | \left( | ||
+ | \begin{array}{cc} | ||
+ | {f_u} & {f_v}\\ | ||
+ | {g_u} & {g_v} | ||
+ | \end{array} | ||
+ | \right) | ||
+ | = | ||
+ | \left( | ||
+ | \begin{array}{cc} | ||
+ | {+} & {-}\\ | ||
+ | {+} & {-} | ||
+ | \end{array} | ||
+ | \right) | ||
+ | \end{eqnarray} | ||
+ | </p> | ||
+ | <p> | ||
+ | situation because...<br> | ||
+ | <br> | ||
+ | $f_u$ must be positive (this condition corresponds to the positive feedback in Turing Pattern).<br> | ||
+ | <br> | ||
+ | $g_u$ must be positive and f_v must be negative(this condition corresponds to the negative feedback in Turing Pattern<br> | ||
+ | |||
+ | </p> | ||
+ | |||
+ | <h3>Explanation about equation</h3> | ||
+ | |||
+ | <p> | ||
+ | We’ll explain the colicin model. For example, case1 is dimentionalized as following. | ||
+ | </p> | ||
+ | |||
+ | <p> | ||
+ | Case 1 | ||
+ | |||
+ | $$ | ||
+ | |||
+ | \frac{du}{dt} = u\{{\alpha} - b(u + v)\} - \frac{a}{a + K_A}{\Delta}u - {\kappa}cu + du{\nabla}^2u \\ | ||
+ | |||
+ | \frac{dv}{dt} = v\{{\beta} - b(u + v)\} + du{\nabla}^2dv \\ | ||
+ | |||
+ | \frac{da}{dt} = {\lambda}_au - {\gamma}a + da{\nabla}^2da \\ | ||
+ | |||
+ | \frac{dc}{dt} = {\lambda}_cv - {\mu}c + dc{\nabla}^2dc \\ | ||
+ | |||
+ | $$ | ||
+ | |||
+ | $u$, $v$, $a$ and $c$ mean activator, co-activator, AHL and colicin, respectively. | ||
+ | </p> | ||
+ | |||
+ | <p> | ||
+ | The first terms of $u$ and $v$ indicate logistically increasing of $u$ and $v$. Alpha and beta are max growth rates of $u$ and $v$, respectively. $b$ indicates decreasing effect caused by $u$ and $v$ scramble for nutrition and we assumed the effects against $u$ and $v$ are same. | ||
+ | </p> | ||
+ | <p> | ||
+ | This is the model based on [1]. The second term in the equation for $u$ indicates death of $u$. | ||
+ | </p> | ||
+ | <p> | ||
+ | We defined the Hill coefficient as 1 because we simplified the Hill coefficient of pLux expression level against AHL, modeled by ETH Zurich 2014, in order to simply model. 4th equation is term of death caused by colicin. We assumed this term is simply in proportion to the concentration of colicin. | ||
+ | </p> | ||
+ | <p> | ||
+ | Death rates of <i><i>E. coli</i></i> caused by AHL and colicin are in proportion to the concentration of <i><i>E. coli</i></i> because <i><i>E. coli</i></i> dies with a certain probability regardless to its concentration. | ||
+ | </p> | ||
+ | <p> | ||
+ | $a$ and $c$ increase proportionally to both concentrations of $u$ and $v$ ach and are degraded proportionally to its concentration. | ||
+ | </p> | ||
+ | <p> | ||
+ | We simplified the equation by deleting ${\lambda}_a, {\lambda}_c, b$. We could nondimensionalize $x$ and $t$ but we didn’t because we wanted to leave the scale of time and location. Variables with hat are new ones. We regarded we nondimensionalized $x$ and $t$ by multiply them by 1. Parameters ${\alpha}, {\beta}, {\gamma}, {\mu}, {\du}, {\dc}, {\da} and {\delta}$ don’t change after nondimentionalized because they’re only depending on $x$ and $t$ | ||
+ | </p> | ||
+ | |||
+ | </div> | ||
+ | |||
+ | |||
+ | <div class='panel' id='parameters'> | ||
+ | <h2>PARAMETERS</h2> | ||
+ | <h3>Colicin model</h3> | ||
+ | <p> At first, we should say that we did not define a strict "parameter", but a "parameter range" for each parameter. As mentioned at the top page, the condition for Turing pattern cannot be satisfied easily. Thus, if a parameter is defined strictly, there is little chance to form pattern. We considered the changeable or fluctuating parameters <br> | ||
+ | Therefore, we defined a paramter range widely. There are no parameters by which patterns are well formed when the parameter is too large or too small , excluding the ratio of diffusion coefficient (the more different, the more likely pattern emerge). All parameters should not be too large, nor too small. Hence we define some range, and within the range, vary parameters, conduct stability analysis, and "seek" the parameter set that generates pattern. <br> | ||
+ | But there needs some calculation time for stability analysis. The system we want to analyze have four elements and ten parameters, and all parameters are related complicatedly. Before stability analysis, we want to assume suitable parameters. We focused on the condition to have a stable point.<br> | ||
+ | For example, on case 2, the condition of having a fixed point ($u_0$, $v_0$, $a_0$, $c_0$) satisfies the condition below.<br> | ||
+ | $$ | ||
+ | {\alpha} - u_0 - v_0 - \frac{a_0}{a_0+K}{\delta} = 0 \\ | ||
+ | {\beta} - u_0 - v0 - {\kappa}c_0 = 0 \\ | ||
+ | u_0 - {\gamma}a_0 = 0 \\ | ||
+ | u_0 - {\mu}c_0 = 0 \\ | ||
+ | u_0, v_0, a_0, c_0 > 0 \\ | ||
+ | $$ | ||
+ | Then, u0 satisfies the condition below. <br> | ||
+ | $$ | ||
+ | -\frac{\kappa}{\mu}u_0^2 +({\alpha} - {\beta} - \frac{\kappa}{\mu}{\gamma}K - {\delta})u_0 + ({\alpha} - {\beta}){\gamma}K = 0 \\ | ||
+ | u_0 > 0 \\ | ||
+ | {\beta} - (1 + \frac{\kappa}{\mu}) * u_0 > 0 (from v_0 > 0) | ||
+ | $$ | ||
+ | From this condition, below can be assumed.<br> | ||
+ | I. ${\alpha}$, ${\beta}$, ${\delta}$ are the same order, from the second term. The term would be useless if the parameters are too different.<br> | ||
+ | II. ${\kappa}$ and ${\mu}$ are the same order, since these two values is always used as ratio.</p><br> | ||
+ | |||
+ | |||
+ | <h3>parameter values:</h3> | ||
+ | The parameters we estimated are on this table.<br> | ||
+ | <table> | ||
+ | <tbody> | ||
+ | <tr class='top’><td>Parameter</td><td>Value</td></tr> | ||
+ | <tr><td> \alpha </td><td> 0.01 ~ 0.7 </td></tr> | ||
+ | <tr><td> \beta </td><td> 0.01 ~ 0.7 </td></tr> | ||
+ | <tr><td> \gamma </td><td> 0.274 ~ 1.19 </td></tr> | ||
+ | <tr><td> \mu </td><td> 0.01 ~ 1.0 </td></tr> | ||
+ | <tr><td> d_a </td><td> 3.0 ~ 3.6 </td></tr> | ||
+ | <tr><td> d_u </td><td> 0.033 ~ 0.034 </td></tr> | ||
+ | <tr><td> d_c </td><td> 0.15 </td></tr> | ||
+ | <tr><td> K </td><td> {(\alpha*\beta)^(1/2)*(0.01 ~ 100.0)}/\gamma </td></tr> | ||
+ | <tr><td> \delta </td><td> (\alpha*\beta)^(1/2)*(0.01 ~ 100.0) </td></tr> | ||
+ | <tr><td> \kappa </td><td> \mu*(0.1 ~ 10.0) </td></tr> | ||
+ | </tbody> | ||
+ | </table> | ||
+ | |||
+ | <p> | ||
+ | The explanations are below.<br> | ||
+ | {\alpha} and {\beta} are the largest growth rate of activator and co-activator respectively. We used least squares method and determined them from the result of experiments. <br> | ||
+ | Max growth rate is (average of plac_barnase without inhibition of barnase).We tooke the logarithm to linearly approximate and gained proportionality constant.The proportionality constant was about 0.18. However, this experiment was under room temperature. Otherwise, the constant would be larger. It was estimated that the constant is 0.7 on this paper[1].<br> | ||
+ | However, those can be changed by introducing reasonable construction; for example, inducing cell death constitutively by killing gene. Scott et. al(2010) said that we can inhibit cell growth by only expressing protein[6]. Controlling the resource of ribosome by expressing gene, growth inhibition would be caused. We defined the minimum of parameter to be 0.01.The concentration is uniformly distributed in the range.<br> | ||
+ | <br> | ||
+ | |||
+ | K is dissociation constant of Hill equation of death term of activator caused by AHL.<br> | ||
+ | <br> | ||
+ | This parameter also can be changed by changing the expression rate of LuxR(if we use LuxAHL), copy number of the promoter, etc...<br> | ||
+ | This can be explained by the origin of Hill's equation, explained by the law of mass action.<br> | ||
+ | <br> | ||
+ | Consider the reaction formula below.<br> | ||
+ | LuxR + AHL ⇆ LuxR-AHL ... 1<br> | ||
+ | LuxR-AHL + promoter ⇆ LuxR-AHL-promoter ... 2<br> | ||
+ | (Naturally, LuxR-AHL dimer binds to promoter and activate this[4] but we can consider this model because we assumed Hill constant is 1.) | ||
+ | K increases in proportion to LuxR.<br> | ||
+ | If we assume there is 1's equilibrium, [LuxR][AHL]/[LuxR-AHL] = const<br> | ||
+ | At the same time, there is also 2's equlibrium, [LuxR][AHL][promoter]/[<br> | ||
+ | <br> | ||
+ | (We indicated that K becomes smaller inversely to the numbers of LuxR and promoter on the condition)<br> | ||
+ | [LuxR-AHL-promoter]/[LuxR]<br> | ||
+ | <br> | ||
+ | There is also the selection of AHL(Las, Lux, H6). ←K is different from this selection. AHL:Refer to experience[2]<br> | ||
+ | <br> | ||
+ | This parameter can also be changed by expression rate of LuxI.<br> | ||
+ | LuxI generate AHL, so increasing LuxI means increasing AHL itself. <br> | ||
+ | On the dimensionalized form(<a name="mathematics">mathematics</a>),<br> | ||
+ | $K = \frac{b}{{\lambda}_a}K_A$ | ||
+ | K_A is dimensionalized and steady value regardless of the amount of AHL if the concentration of LuxR is steady.<br> | ||
+ | Increasing the expression rate of AHL means increasing ${\lambda}_A$ and K decreases.<br> | ||
+ | <br> | ||
+ | Therefore, many variations can be considered.<br> | ||
+ | We considered the parameter from {\gamma}, {\alpha} and {\beta} that are already found to some extent. <br> | ||
+ | specifically,divide the geometric means of values of {\alpha} and {\beta} of each parameter set by {\gamma}. And then we mulitiplied it by random numbers which range from 0.1 to 100. But the random numbers here aren't unifomal because the ratio is important this time.<br> | ||
+ | <br> | ||
+ | <h4>delta</h4> | ||
+ | delta is the maximum value of <i>E. coli</i> death rate caused by AHL<br> | ||
+ | Varying promoters such as pLux or pLas, you can control the value of this parameter. According to ETH Zurich 2014, crosstalk by promoters happens[2]. Because the expression caused by crosstalk is relatively weak, you can adjust the expression of delta by taking advantage of this.<br> | ||
+ | <br> | ||
+ | According to the result of the experiment, we succeded in reducing the reproduction rate to its half value by the growth inhibition of barnase on highly expressed pLac. We can reduce the rate more by inducing other lethal genes.<br> | ||
+ | <br> | ||
+ | We defined the range of this value as geometric mean of alpha and beta multiplied by a random number from 0.1 to 1.0.<br> | ||
+ | |||
+ | |||
+ | |||
+ | <h4>kappa</h4> | ||
+ | |||
+ | The parameter corresponds to the constant part of the term representing the death of <i>E. coli</i> by Colicin.<br> | ||
+ | This can be changed by expression rate of colicin (just as like K).<br> | ||
+ | We defined the range of the parameter as mu multiplied by 0.1 ~ mu multiplied 10. We chose the value from this range randomly.<br> | ||
+ | |||
+ | <h4>gamma</h4> | ||
+ | |||
+ | This can be changed by pH or temperature[10]. Weiss et. al[9] estimated this value to be between 0.274 to 1.19,<br> | ||
+ | within the range of pH where the reproduction of <i>E. coli</i> is not interfered.We chose the value of the parameter from this range randomly.<br> | ||
+ | |||
+ | <h4>mu</h4> | ||
+ | |||
+ | The degradation rate of colicin<br> | ||
+ | We could not find the value of this. We also have not certificate if this can be changed. iGEM Heidelberg 2008[5] said that the degradation of colicin can be ignored; however, its scale was second, not hour. They also said that the death rate of the AHL can also be ignored. However, as shown above, it cannot be ignored in a scale of So we set this value as the same or less than AHL. This parameter also may be changed by adding protease to the medium.<br> | ||
+ | <br> | ||
+ | We want to take the range of mu as broad as possible, because only the value of the ratio to kappa is provided. Too large or too small value is not appropriate for the pattern formation. Considering this, we chose the value of the parameter from 0.01 to 1.0 randomly.<br> | ||
+ | |||
+ | <h4>$d_a$</h4> | ||
+ | |||
+ | the diffusion coefficient of AHL<br> | ||
+ | According to Berg, Brown(1972), 3.2 [mm^2/h]. Around its value (some fluctuation on the condition of experiment) →[3.0, 3.6]<br> | ||
+ | |||
+ | <h4>$d_u$</h4> | ||
+ | |||
+ | the diffusion coefficient of <i>E. coli</i><br> | ||
+ | we estimated it from experiment ; that was0.0337<br> | ||
+ | We defined 0.033~0.034 as parameter range.<br> | ||
+ | |||
+ | <h4>$d_c$</h4> | ||
+ | the diffusion coefficient of colicin<br> | ||
+ | <br> | ||
+ | According to Schwartz, Helinski, 1971, 0.15 [mm^2/h] However, it is not the condition of semi-agarose gel. Hence, we decided the parameter range 0.05~0.5, around its value.</p> | ||
+ | |||
+ | |||
+ | |||
+ | <h3>Future work: further model estimation</h3> | ||
+ | <p> | ||
+ | There are a lot of parameters that have not been determined yet. However, some of them can be determined by experiment.<br> | ||
+ | Parameters on AHL(K, delta) can be calculated just like below.<br> | ||
+ | |||
+ | <h3>Method</h3> | ||
+ | 1. Construct the first strategy. <br> | ||
+ | 2. Varying pH (or construction itself), measure the concentration at the equilibrium point.<br> | ||
+ | In a article, they estimated the parameter gamma in this way [9] . However, its model is different from ours.<br> | ||
+ | |||
+ | If we can find this point from the experiment, some parameters can be estimated. In the experiment, we vary the expression rate of colicin by the selection of promoter. This relates varying ${\kappa}$; from dimensionalized form, when the expression rate increase, ${\kappa}$ would be decreased. Prepare two or more systems that have different promoter, and there appears two equation, which kappa is different from each other. ${\kappa}$ is inverse proportion to the strength of the promoter, so the relationship between those two different ${\kappa}$ can be estimated. Since , there arise three equations, and we can determine those three values.<br> | ||
+ | |||
+ | |||
+ | </div> | ||
+ | |||
+ | <div class='posnav' id='breadcrumbs-one'></div> | ||
+ | </div> | ||
+ | </body> | ||
+ | <script type="text/javascript" src="http://ajax.googleapis.com/ajax/libs/jquery/1.11.1/jquery.min.js"></script> | ||
+ | <script type='text/javascript' src='./js/isotope.pkgd.js'></script> | ||
</html> | </html> |
Latest revision as of 03:55, 19 September 2015
TOP
Overview
To generate Turing pattern using synthetic biology, we thought quite a few systems. We analyzed some of them and checked if they can generate a pattern and how the pattern can be. We used some differential equations and multi-agent simulation for this. From modeling, it was suggested that it is difficult to find the parameters that certainly satisfies the condition on pattern formation. However, some trends on parameters were observed.
The aim of the modeling
Turing mechanism is based on a mathematical model (mathematics page). We modeled our system in order to understand it mathematically as well as to find the parameter that satisfies the condition for forming a pattern.
Background
In order to understand our system and analysis, some basic knowledge for Turing pattern is needed. We’d like to introduce some mathematical basics of Turing Pattern to make it easier to understand our project in the mathematics page.
The flow of whole modeling project
Analysis and explanation of pattern generation
We first showed that the first system cannot form a pattern. To solve this problem, we decided to use colicin. Why we choosed to use colicin is qualitatively explained here.
System Analysis
As we saw above, it is difficult to satisfy the conditions for forming Turing pattern. Thus, at first we tried to weaken this conditions, to find the paramerer range which would form the pattern.
We succeeded in determining the range, however, due to the fructuations, it turnd out that they no longer certainly generate Turnig pattern.
How the pattern will be in the last
In this part, we simulated how the pattern wil be. we constructed a mathematical model which describe our system, and applied the sets of parameters found in the previous analysis.
we did simulation on on dimension and two dimension case.
one dimension:
The analysis on one dimension. Assuming that E. coli are in a long and thin container. Stripe pattern would appear, if the optimization of part 1 is well conducted. We also checked how long it would take to form the pattern, and whether a periodical pattern would appear.
two dimension:
Generally, the model showed various patterns: spot, stripe, competition,... etc. in two dimension case.
We used the set of parameters obtained in the first analysis.
in this model, we assumed the container large enough to ignore effescs from the boundary: i.e. the diffusion equation holds for this system.
The simulation under this assumption ended only in activator spot pattern.
New model
We also conducted the analysis above to a new model. It is inspired by zebrafish and includes the effect of colicin.
CODE see this GitHub page. We modeled first and second sections by MATLAB, and the last one by Java.[1] Instabilities in spatially extended predator–prey systems Spatio-temporal patterns in the neighborhood of Turing–Hopf bifurcations, Baurmann et. al 2006 Journal of Theoretial Biology 245 (2007) 220–229
System
First, we will explain why the first model is not suitable.
The reason why the first model cannnot seem to form pattern
The interactions between E. coli and AHL are expressed as two differential equations below.
$$
\frac{du}{dt}=d_u\nabla^2u+u(1-u)-{\frac{a}{a+K}}{\delta}u=d_u\nabla^2u+f(u,a)\\
\frac{da}{dt}=d_a\nabla^2a+u-{\gamma}a=d_a\nabla^2a+g(u,a)
$$
u: concentration of E. coli
a: concentration of AHL
d_u: the diffusion coefficient of u.
d_a: the diffusion coefficient of a.
{\gamma}: degradation coefficient of AHL
(Note that the parameters are non-dimensionized.)
equation 1
The concentration of E. coli increases by multiplication in proportion to its concentration, but if the concentration is too high, the growth rate decreases because of the lack of nutrition(+u(1-u)). If E. coli receives enough AHL, it is killed by Barnase {\frac{a}{a+K}}{\delta}u.The reaction term {\frac{a}{a+K}}{\delta} can be approximated by Hill equation because the system of E. coli death is similar to ligand-receptor system of protein.
equation2
The concentration of AHL increases in proportion to that of E. coli, and AHL degrades in proportion to -{\gamma}.
We concluded that this model may not be able to form a pattern because this model doesn’t satisfy required mathematical conditions of Turing pattern.
Proof
Following two equations should be satisfied.
\begin{equation}
f_u+g_a\lt 0
\label{eq:vern1}
\end{equation}
\begin{equation}
f_ug_a-f_vg_u\gt 0
\label{eq:vern2}
\end{equation}
First, there are two possible equilibrium points \(\frac{du}{dt}= \frac{da}{dt}= 0\).
\((u,a)=(0,0)\) and \((u,a)=(u_0,a_0)\)
\(u_0\) and \(a_0\) satisfies the relations,\(u_0 = \gamma a_0\) and \(f(a_0)=1− u_0\)
\((u,a)=(0,0)\) doesn’t satisfy equation [\ref{eq:vern1}] and [\ref{eq:vern2} ].
When \((u,a)=(u_0,a_0)\),
$$
f_u=1 – 2u_0 – {\frac{a}{a+K}}{\delta}=1 – 2u_0 – (1 – u_0)= –u_0 \lt 0
$$
$$
f_a= –u_0・d_u{\frac{f(a)}{da}}|a_0 \lt 0 \qquad (\mathrm{for} \ d_u\frac{f(a)}{d_a} \gt 0)$$
because \(f(a)\) is approximated by Hill’s equation.
\(g_u= 1\)
\(g_a= –{\gamma}\lt0\)
Putting these values into equation \ref{eq:vern1} and \ref{eq:vern2} ,
\(f_u+g_a= – u_0 –\gamma\lt 0\)
\(f_ug_a−f_ag_u = u_0(\gamma+ d_u\frac{f(a)}{d_a}|a_0)\gt0\)
So, \((u,a)=(u_0,a_0)\) satisfies those requirements.
To form Turing pattern, however, there is still one more condition to be satisfied two equations below.
\((d_ug_a + f_u)^2 – 4d_u(f_ug_a−f_ag_u) \gt 0\)
\(d_ug_a + d_af_u\gt 0\)
For \(g_a\lt 0, f_u\lt 0\), equation 6 is never satisfied whatever value d_u takes.
Requirement 2 unsatisfied, this model cannot form Turing pattern.
To meet the equation 6, \(f_u\) must be bigger than zero. And to satisfy this, the positive feedback
Therefore, there must be a system that has allee effect, or density dependant effect. The former means that "When E. coli increases, the growth rate of E. coli also increases."Allee effect is a positive correlation between population size and indivisual fitness. Mimura and Murray(1972) said that Allee effect is important for some pattern formation on predator-prey systems[1]. However, there is no Allee effect in E. coli. E. coli increase however little amount they are(one colony would be formed from one cell). The multiplication rate of E. coli is small when it’s introduced into culture medium, for it’s not adapted to the environment. But that doesn’t affect the formation of Turing pattern, because the equilibrium point of Turing pattern occurs on log phase or stationary phase.The latter is saying that "When there are very much E. coli, the effect of AHL is not proportional to E. coli's density"; the effect will be lower. This is seen in the model of "Prey-Predator model". Since predator eats no more than they can eat, when preys are rich, the probablity of being killed of an indivisual prey would be decreased. However, since AHL is much smaller and richer than E. coli, there seems little chance to expect it to happen.
It seemed that both cannot be satisfied only by one type of E. coli. Therefore, we decided to solve this problem by using two types of E. coli and made positive-feedback-like system. Most positive feedback system have bistability (which does not exist in the system of population of E. coli). In other words, the space would be separated by two areas, according to the convergence value. The convergence value is polarized. This is a kind of the system of the competition exclusion. To learn more about the principle of competition exclusion, see Murray, mathematical biology volume 1.
We wanted to realize competition exclusion system. We chose the system that utilizes colicin. See the system page for detail.
Why colicin? Because colicin does not diffuse less quickly than AHL. Short-range positive feedback and long-range negative feedback are necessary for Turing pattern. Hence, bistability should be achieved locally.
Chao and Levin says that this colicin-using system (no AHL effect on our system) exhibit the bistability[3]. We show this in the mathematics page.
We considered two cases.
Case 1. E. coli that emit AHL(shown u in figure) and E. coli that emit colicin(shown v in figure) are different
formula
$$
\frac{du}{dt} = u({\alpha} - (u + v) - {\kappa}c - \frac{a}{a + K}{\delta}) + du{\nabla}^2u\\
\frac{dv}{dt} = v({\beta} - (u + v)\) + du{\nabla}^2v \\
\frac{da}{dt} = u - {\gamma}a + da{\nabla}^2a \\
\frac{dc}{dt} = v - {\mu}c + dc{\nabla}^2c \\
$$
Here, AHL, E. coli, colicin are all scaled, and unnecessary parameters are deleted.
Case 2. E. coli that emit AHL(shown u in figure) and E. coli that is killed by colicin(shown v in figure) are same.
formula
$$
\frac{du}{dt} = u({\alpha} - (u + v) - \frac{a}{a + K}{\delta}) + du{\nabla}^2u \\
\frac{dv}{dt} = v(\beta - (u + v) - {\kappa} c) + du{\nabla}^2v \\
\frac{da}{dt} = u - {\gamma}a + da{\nabla}^2a \\
\frac{dc}{dt} = u - {\mu}c + dc{\nabla}^2c \\
$$
We first investigated "which parameter set can lead to form pattern" by stability analysis.
First, we defined parameter range. See the "parameter" page to see the detail.
Second, we varied parameters according to parameter range, and constructed some parameter sets.
Third, for each parameter set, we saw if it satisfies the condition of Turing pattern, mentioned here (mathematics). If it is satisfied, it is registered as a "passed parameter set". Each parameter set is also characterized by its largest eigenvalue of stability matrix, with and without diffusion term, sensitivity, and the wavelength that takes the maximum eigenvalue.
Then, after seeking parameters, we saw some trend of parameters. We also ranked parameter sets, acording to their stability matrixes' eigenvalue. Since E. coli increases in the order of hour, it is considered to take a lot of time for pattern formation. However, if it takes much time for pattern formation, there may arise some difficulty on the experiment that the nutrients for E. coli would be run out until the pattern is formed.
There are 10 parameters in the last formulas. However, there were some parameters that do not seem to have been found. There were also some parameters that can be changed by experiment; strength of the constitutive promoter, a copy number of the promoter, etc... For those parameters, we defined a parameter range so that each parameter can affect forming a pattern. See parameter page for detail.
Result
two cases on the same parameter range, it turned out that more patterns were formed in Case 1 than Case 2. And the biggest eigenvalue of Case 1 was smaller than that of Case 2.
The trend of parameters
Case 1
・Set same parameter range, Case 1 is less likely to make pattern than Case 2.
・alpha/beta > 1 should be satisfied. In case 2, however, alhpa should not be very large compared to beta.
・K should not be so large, nor so small.
・delta should be larger than 1, so that it can kill u or v. This means that
・kappa should be large, especially if you want to make pattern fast (from the relationship between the eigenvalue).
・there seems small sensitivity in gamma.
・mu should be large, so that the order is near alpha and beta. mu should not be so large(larger than 0.7).
・du, dc da have no trend as long as its range is narrow enough
・i think there would be a relationship between K and delta. I will check it.
・ delta should be large and K should be small.
・ K and delta is related. If K is large, delta must not be so small.
x-axis alpha, y-axis beta for case 2.
Sensitivity and Robustness
We've tested whether pattern forms or not if we change the parameters, and what parameters would exert effects on the formation. Here are the details:
There are two points of view in sensitivity; local or global. This classification was suggested by iGEM Waterloo 2014 team. Local sensitivity is a test about how the result turns by changing only one parameter. Global one sees the same thing by changing various parameters. We've done the two tests and examined what parameter sets would satisfy the conditions of Turing Pattern.
case 2
We observed strong polarization along the eigenvalue gets bigger or smaller. Polarization of 50 of those, which have small eigenvalue was not so stronger than previous examination. The more eigenvalue get bigger, the more strong polarization occurs. In other words, it can be said that something which makes eigenvalue bigger exists in the place where alpha's max and beta's minimum. The max place where each parameters take would be the boundary line of whether pattern is formed or not. If the parameter go beyond the boundary line, parameter sets couldn't have stable point or could not be stable, because of no occurrence of diffusion. Thus, we believe there would be trade-off between size of eigenvalue and robustness. We analyzed the stability of the random numbers around the parameter (parameter*0.5~parameter*1.5) and examined if there was deviation of distribution when pattern formed. The distribution should be linearly condition if there was no deviation. Second, we simulated Turing pattern numerically. in two dimension, and saw if it forms a pattern. See this figure. Here, you can see pattern is formed.
System2
Another system
Here, we explain the "third" model. This is similar to the colicin-using model. There is some difference. The third strategy is inspired by the pattern formation mechanism of zebrafish[1]. Two types of E. coli (activator and inhibitor) react to each other in two manners according to the distance between them. When E. coli (activator) and E. coli (inhibitor) are close, they inhibit each other in the same manner as in the conventional system. This control loop functions as local activation. In addition, E. coli (activator) activates E. coli (inhibitor) through AHL in long distance (detailed scheme is explained in the chapter of System). The control loop composed of this activation and the short-range inhibition of E. coli (activator) by activated E. coli (inhibitor) functions as lateral inhibition. "Activator" E. coli promotes the multiplication of "Inhibitor" E. coli, and "Inhibitor" inhibits the multiplication of "Activator" near of "Inhibitor" First, we show that there are two cases for this modeling. Case 3: E. coli that emit AHL(shown u in figure) and E. coli that emit colicin(shown v in figure) are different. formula $$ \frac{du}{dt} = u({\alpha} - (u + v) - {\kappa}c) + du{\nabla}^2u\\ \frac{dv}{dt} = v({\beta} - (u + v)\) + \frac{a}{a + K}{\delta}) + du{\nabla}^2v \\ \frac{da}{dt} = u - {\gamma}a + da{\nabla}^2a \\ \frac{dc}{dt} = v - {\mu}c + dc{\nabla}^2c \\ $$ Case 4; E. coli that emit AHL(shown u in figure) and E. coli that emit colicin are same. formula $$ \frac{du}{dt} = u({\alpha} - (u + v)) + du{\nabla}^2u \\ \frac{dv}{dt} = v({\beta} - (u + v) - {\kappa}c + \frac{a}{a + K}{\delta}) + du{\nabla}^2v \\ \frac{da}{dt} = u - {\gamma}a + da{\nabla}^2a \\ \frac{dc}{dt} = u - {\mu}c + dc{\nabla}^2c \\ $$ The parameter range which is different from strategy 2 is delta. delta + beta must be lower than max value of alpha or beta, because barstar or colI can't restore the growth rate of E. coli more than its max value. Hence, the sum of ${\delta}$ and ${\beta}$ must be the same or below the maximum growth rate. The method is not different from the first one. [2]Nakamasu, A., Takahashi, G., Kanbe, A., & Kondo, S. (2009). Interactions between zebrafish pigment cells responsible for the generation of Turing patterns. Proceedings of the National Academy of Sciences, 106(21), 8429-8434.
In order to understand our system and analysis, some basic knowledge for Turing pattern is needed. We’d like to introduce some mathematical basics of Turing Pattern to make it easier to understand our project.
1.Stability
The idea of “stability” is so important in our modeling.Consider conditions of $ \frac{du}{dt} = f(u, v), f(u_0, v_0)=0\\ \frac{dv}{dt} = g(u, v), g(u_0, v_0)=0 $ where u is the concentration of activator, v is of inhibitor. So both u and v are functions of time(t) and location.
Suppose $(u, v) = (u_0, v_0)$ and a fluctuation then occurs to the point. Let $(u_0 + u_s, v_0 + v_s)$ be the point after the fluctuation occurred. Here we can write the differential equations as $ \frac{du_s}{dt} = f(u_0 + u_s, v_0 + v_s) \\ \frac{dv_s}{dt} = g(u_0 + u_s, v_0 + v_s) $
Given $u_s$, $v_s$ are small, we can expand $f$,$g$ around $(u_0, v_0)$ excluding the terms whose degree is over one.
That is $ \frac{du_s}{dt} = f(u_0 + u_s, v_0 + v_s) = f_uu_s + f_vv_s\\ \frac{dv_s}{dt} = g(u_0 + u_s, v_0 + v_s) = g_uu_s + g_vv_s $ where $ f_u = f_u(u_0, v_0),\\ f_v = f_v(u_0, v_0)\\ g_u = g_u(u_0, v_0),\\ g_v = g_v(u_0, v_0) $ (From now on, we mention $f_u$, $f_v$, $g_u$, $g_v$ as $f_u(u_0, v_0)$, $f_v(u_0, v_0)$, $g_u(u_0, v_0)$, $g_v(u_0, v_0)$ until the end for convenience)
In another way, \begin{eqnarray} \left( \begin{array}{cccc} \frac{du_s}{dt}\\ \frac{dv_s}{dt} \end{array} \right) &=& \left( \begin{array}{cc} {f_u} & {f_v}\\ {g_u} & {g_v} \end{array} \right) \left( \begin{array}{cc} {u_s}\\ {v_s} \end{array} \right) ...(A) \end{eqnarray}
We want to find the solutions of $(u_s, v_s)$.
Given both $u_s$ and $v_s$ are in proportion to $e^{{\lambda}t}$, it can be deduced that ${\lambda}$ is equal to the eigenvalues of the matrix in (A). And as the matrix in (A) has two eigenvalues, it can be said that $u_s$ and $v_s$ are described as a linear combination of $e^{{\lambda}t}s$. If all of ${\lambda}$s are negative, $e^{{\lambda}t}$s get smaller and smaller as time passes. It means every time when (u,v) feels any small fluctuations, (u,v) returns to the point $(u_0, v_0)$. In other words, if all ${\lambda}$s are negative, $(u_0, v_0)$ is a “stable” point. On the contrary, if any of ${\lambda}$s is positive, as $e^{{\lambda}t}$s whose ${\lambda}$ is positive gets larger and larger, it can be said that (u, v) won’t return to $(u_0, v_0)$. Then $(u_0, v_0)$ isn’t a “stable“ point. Conditions for a point $(u_0, v_0)$ to be stable are known as $ f_u + g_v < 0\\ f_ug_v - f_vg_u > 0 $ These conditions can be expanded to cases of three and more factors.
2. Turing Pattern
For Turing Pattern, such points are required that satisfy the conditions of…
- ・stable when dilution isn’t taken into account
- ・not stable when dilution is taken into account
For the point of $(u_0, v_0)$ to be unstable, it’s necessary that $ (D_ug_v + D_vf_u)^2 - 4D_uD_v(f_ug_v - f_vg_u)> 0 \\ D_ug_v + D_vf_u > 0 $
This means that the larger the ratio of $D_u$ and $D_v$ is, the more the point considered is likely to be unstable. To explain more detail, it's deduced that \begin{eqnarray} \left( \begin{array}{cc} {f_u} & {f_v}\\ {g_u} & {g_v} \end{array} \right) = \left( \begin{array}{cc} {+} & {+}\\ {-} & {-} \end{array} \right) \ or \left( \begin{array}{cc} {+} & {-}\\ {+} & {-} \end{array} \right) \end{eqnarray} There're two possible situations as above,but we consider only the \begin{eqnarray} \left( \begin{array}{cc} {f_u} & {f_v}\\ {g_u} & {g_v} \end{array} \right) = \left( \begin{array}{cc} {+} & {-}\\ {+} & {-} \end{array} \right) \end{eqnarray}
situation because...
$f_u$ must be positive (this condition corresponds to the positive feedback in Turing Pattern).
$g_u$ must be positive and f_v must be negative(this condition corresponds to the negative feedback in Turing Pattern
Explanation about equation
We’ll explain the colicin model. For example, case1 is dimentionalized as following.
Case 1 $$ \frac{du}{dt} = u\{{\alpha} - b(u + v)\} - \frac{a}{a + K_A}{\Delta}u - {\kappa}cu + du{\nabla}^2u \\ \frac{dv}{dt} = v\{{\beta} - b(u + v)\} + du{\nabla}^2dv \\ \frac{da}{dt} = {\lambda}_au - {\gamma}a + da{\nabla}^2da \\ \frac{dc}{dt} = {\lambda}_cv - {\mu}c + dc{\nabla}^2dc \\ $$ $u$, $v$, $a$ and $c$ mean activator, co-activator, AHL and colicin, respectively.
The first terms of $u$ and $v$ indicate logistically increasing of $u$ and $v$. Alpha and beta are max growth rates of $u$ and $v$, respectively. $b$ indicates decreasing effect caused by $u$ and $v$ scramble for nutrition and we assumed the effects against $u$ and $v$ are same.
This is the model based on [1]. The second term in the equation for $u$ indicates death of $u$.
We defined the Hill coefficient as 1 because we simplified the Hill coefficient of pLux expression level against AHL, modeled by ETH Zurich 2014, in order to simply model. 4th equation is term of death caused by colicin. We assumed this term is simply in proportion to the concentration of colicin.
Death rates of E. coli caused by AHL and colicin are in proportion to the concentration of E. coli because E. coli dies with a certain probability regardless to its concentration.
$a$ and $c$ increase proportionally to both concentrations of $u$ and $v$ ach and are degraded proportionally to its concentration.
We simplified the equation by deleting ${\lambda}_a, {\lambda}_c, b$. We could nondimensionalize $x$ and $t$ but we didn’t because we wanted to leave the scale of time and location. Variables with hat are new ones. We regarded we nondimensionalized $x$ and $t$ by multiply them by 1. Parameters ${\alpha}, {\beta}, {\gamma}, {\mu}, {\du}, {\dc}, {\da} and {\delta}$ don’t change after nondimentionalized because they’re only depending on $x$ and $t$
PARAMETERS
Colicin model
At first, we should say that we did not define a strict "parameter", but a "parameter range" for each parameter. As mentioned at the top page, the condition for Turing pattern cannot be satisfied easily. Thus, if a parameter is defined strictly, there is little chance to form pattern. We considered the changeable or fluctuating parameters
Therefore, we defined a paramter range widely. There are no parameters by which patterns are well formed when the parameter is too large or too small , excluding the ratio of diffusion coefficient (the more different, the more likely pattern emerge). All parameters should not be too large, nor too small. Hence we define some range, and within the range, vary parameters, conduct stability analysis, and "seek" the parameter set that generates pattern.
But there needs some calculation time for stability analysis. The system we want to analyze have four elements and ten parameters, and all parameters are related complicatedly. Before stability analysis, we want to assume suitable parameters. We focused on the condition to have a stable point.
For example, on case 2, the condition of having a fixed point ($u_0$, $v_0$, $a_0$, $c_0$) satisfies the condition below.
$$
{\alpha} - u_0 - v_0 - \frac{a_0}{a_0+K}{\delta} = 0 \\
{\beta} - u_0 - v0 - {\kappa}c_0 = 0 \\
u_0 - {\gamma}a_0 = 0 \\
u_0 - {\mu}c_0 = 0 \\
u_0, v_0, a_0, c_0 > 0 \\
$$
Then, u0 satisfies the condition below.
$$
-\frac{\kappa}{\mu}u_0^2 +({\alpha} - {\beta} - \frac{\kappa}{\mu}{\gamma}K - {\delta})u_0 + ({\alpha} - {\beta}){\gamma}K = 0 \\
u_0 > 0 \\
{\beta} - (1 + \frac{\kappa}{\mu}) * u_0 > 0 (from v_0 > 0)
$$
From this condition, below can be assumed.
I. ${\alpha}$, ${\beta}$, ${\delta}$ are the same order, from the second term. The term would be useless if the parameters are too different.
II. ${\kappa}$ and ${\mu}$ are the same order, since these two values is always used as ratio.