Team:UT-Tokyo/Modeling
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.