Revision as of 07:56, 8 October 2018

Model

As the aim of our project is to construct a modular toolkit using co-cultures with three organism, we would have to conduct many experimental iterations to reach conclusive and reliable results. To avoid this tedious process, we, as many other scientists working in this field, turned our attention to computer based biological modeling to predict the behavior of our system.
In the realm of mathematical and computer based biology, scientists use algorithms, data structure and data visualization to create approximations of existing biological systems. This serves two important goals; One, to check wether the theoretical understanding of a process is correct. Of course, this can be done by constructing a model and checking if the cultures behave within reasonable correctness as the observations in the lab indicate.
The other goal of such modeling approach is to predict wether the behavior of the system are changed or other external influences are applied. This latter role is why we are using modeling in our project.
However, many of our team members lacked experience in this field. As such, we approached the Institute of Quantitative and Theoretical Biology at our university (Headed by Professor Dr. Oliver Ebenhöh) to discuss our approach.
The behavior of a coculture is complex. Modeling in it’s entirety would be far too complex. As such we had to find a way to reduce this complexity. To achieve this, we began our work with laying out a few assumptions that would reduce the complexity of our system and give us a clear indicator of what went wrong if a prediction turned out to be incorrect. In our model we assumed that:

The consumption of metabolites, especially for growth, occurs in discreet units.
All organisms have equal access to the metabolites.
There are no organism specific limits to metabolite uptake or release.
Every metabolite are constant and optimal except the ones we control with.
Growth as well as metabolic reactions are limited by low avalability as well as the slowest enzyme reaction, respectively.
There is no change in the speed of reactions.
The capacity of the system for each organism is set and not calculated dynamically.
While phosphate and nitrogen are important for growth, only a complete lack of carbon (as a stand-in for sugar in our model) will result in starvation of the organisms.

Assumption 1 is a consequence of the way our model is calculated, the others are there to reduce complexity. As stated above, this reduced complexity also reflects a reduction in realism, however, our model has produced results that in our mind, make a more complex approach unnecessary for now.

Our model is iterative and goes though the following equations for each organism and metabolite each iteration:

$a d d = {\begin{matrix} α * ((m_{c} * (Y_{e} [i] / Y_{e} [0]) - i_{c}) / m_{c}) & C > α * ((m_{c} * (Y_{e} [i] / Y_{e} [0]) - i_{c}) / m_{c}) \\ C & C \leq α * ((m_{c} * (Y_{e} [i] / Y_{e} [0]) - i_{c}) / m_{c}) \end{matrix}$
$i_{c} = i_{c} + a d d$
$c = c - a d d$
$g r o w t h f_{c} = {\begin{matrix} i_{c} / (u s e_{c} * (Y_{e} [i] / Y_{e} [0])) & i_{c} / (u s e_{c} * (Y_{e} [i] / Y_{e} [0])) \leq 1 \\ 1 & i_{c} / (u s e_{c} * (Y_{e} [i] / Y_{e} [0])) > 1 \end{matrix}$
$i_{c} = i_{c} - g r o w t h f_{c} * u s e_{c} * (Y_{e} [i] / Y_{e} [0])$
$g r o w t h = {\begin{matrix} g r o w t h f_{c} & g r o w t h f_{p} > g r o w t h f_{c} < g r o w t h f_{n} \\ g r o w t h f_{n} & g r o w t h f_{p} > g r o w t h f_{n} < g r o w t h f_{c} \\ g r o w t h f_{p} & g r o w t h f_{c} > g r o w t h f_{p} < g r o w t h f_{n} \end{matrix}$
$Y_{e} [i + 1] = Y_{e} [i] + g r o w t h * ((c a p a c i t y - Y_{e} [i]) / c a p a c i t y)$

Here, add represents the additional metabolite influx per iteration. &alpha is the organism and metabolite specific uptake speed, m_c is the maximum possible concentration of the specific metabolite in the organism, Y_e[i] is the current dry weight (in mol) of the organism and Y_e[0] is the initial dry weight. i_c is the current metabolite mass inside the organism (in mol) and C is the mass of metabolite (in mol) inside the medium. use_c is the specific metabolite mass the organisms the organism uses each iteration in it’s growth, whereas growthf_c, growthf_n, and growthf_p are the metabolite specific growth factors. The lowest of these is chosen as the overall growth factor growth and then used in the calculation of the new population size Y_e[i+1].

As stated above, these calculation are performed for each organism and each metabolite. the resulting populations lists are then plotted via the Python package Matplotlib (as you can see in the code we provided at the end of the page).
This produces a graph that shows what the growth and population sizes of our system look like within a predetermined time frame. This tells us which ratio of metabolite concentration to population sizes is necessary for the system to remain stable over a given period of time. This gives us a range of concentrations that we can then verify in the lab. This reduces the amount of experiments we have to conduct to reach our goal of establishing a stable, purpose-engineered three organism coculture.
Of course observation in the lab take precedent over our model. And data that is in conflict with the prediction of our model will always be given priority. Until now however, our model has been correct. Thus far, it has revealed that the same problem is the supplying the system with sugar, as S. elongatus is has the lowest doubling time and appears to be the main bottleneck. Below a few example models are supplied, with varying amounts of metabolite supplied at the beginning.

Here are a few examples where we used our three model organisms E. coli, S. cerevisae and S. elongatus and plotted their growth over the course of 48 hours. Changed are the concentrations of the metabolites we control with, in our example this is carbon, nitrogen and phosphor, in the medium, but everything else was kept at the same value;

Every metabolite is set at 10 Mol initially:
Every metabolite is set at 50 Mol initially:
Every metabolite is set at 100 Mol initially:

As stated above, each of our assumption attempts to reduce complexity also limits realism. If given more time we could further refine our model and start working without our assumptions. This would increase the accuracy of our predictions and would further our understanding of cocultures and microbe interaction.
It would however also necessitate more experiments to verify our predictions, which was something we didn’t have the time for.
Another possibility we could not pursue was interactivity and the option to include other organisms in our model. We hoped that we could construct a model that would receive a selection of organisms, metabolites that bind the organisms together and a time frame and from that calculate the optimal medium composition as well as initial population sizes to keep the coculture stable.
We believe that such a thing is possible with our model and if we would have had more time to verify our modeling of our lab work, we could have realized this idea.

Code

Provided below is the Python code that makes up our model:
(Note that we used the jupyter Notebook Editor for our work and we recommend using it as well if you want to try out our model!)

import matplotlib.pyplot as plt
import numpy as np
import scipy.integrate

def get_itc():
global itc
return itc
def set_itc(new):
global itc
itc = new
return

def get_itn():
global itn
return itn
def set_itn(new):
global itn
itn = new
return

def get_itp():
global itp
return itp
def set_itp(new):
global itp
itp = new
return

def get_itm():
global itm
return itm
def set_itm(new):
global itm
itm = new
return

def get_itpa():
global itpa
return itpa
def set_itpa(new):
global itpa
itpa = new
return

def get_i_c_e():
global i_c_e
return i_c_e
def set_i_c_e(new):
global i_c_e
i_c_e = new
return

def get_i_c_s():
global i_c_s
return i_c_s
def set_i_c_s(new):
global i_c_s
i_c_s = new
return

def get_i_c_se():
global i_c_se
return i_c_se
def set_i_c_se(new):
global i_c_se
i_c_se = new
return

def get_i_n_e():
global i_n_e
return i_n_e
def set_i_n_e(new):
global i_n_e
i_n_e = new
return

def get_i_n_s():
global i_n_s
return i_n_s
def set_i_n_s(new):
global i_n_s
i_n_s = new
return

def get_i_n_se():
global i_n_se
return i_n_se
def set_i_n_se(new):
global i_n_se
i_n_se = new
return

def get_i_p_e():
global i_p_e
return i_p_e
def set_i_p_e(new):
global i_p_e
i_p_e = new
return

def get_i_p_s():
global i_p_s
return i_p_s
def set_i_p_s(new):
global i_p_s
i_p_s = new
return

def get_i_p_se():
global i_p_se
return i_p_se
def set_i_p_se(new):
global i_p_se
i_p_se = new
return

def get_i_m():
global i_m
return i_m
def set_i_m(new):
global i_m
i_m = new
return

def get_i_pa():
global i_pa
return i_pa
def set_i_pa(new):
global i_pa
i_pa = new
return

def get_m_c_e():
global m_c_e
return m_c_e
def set_m_c_e(new):
global m_c_e
m_c_e = new
return

def get_m_c_s():
global m_c_s
return m_c_s
def set_m_c_s(new):
global m_c_s
m_c_s = new
return

def get_m_c_se():
global m_c_se
return m_c_se
def set_m_c_se(new):
global m_c_se
m_c_se = new
return

def get_m_n_e():
global m_n_e
return m_n_e
def set_m_n_e(new):
global m_n_e
m_n_e = new
return

def get_m_n_s():
global m_n_s
return m_n_s
def set_m_n_s(new):
global m_n_s
m_n_s = new
return

def get_m_n_se():
global m_n_se
return m_n_se
def set_m_n_se(new):
global m_n_se
m_n_se = new
return

def get_m_p_e():
global m_p_e
return m_p_e
def set_m_p_e(new):
global m_p_e
m_p_e = new
return

def get_m_p_s():
global m_p_s
return m_p_s
def set_m_p_s(new):
global m_p_s
m_p_s = new
return

def get_m_p_se():
global m_p_se
return m_p_se
def set_m_p_se(new):
global m_p_se
m_p_se = new
return

def get_m_m():
global m_m
return m_m
def set_m_m(new):
global m_m
m_m = new
return

def get_m_pa():
global m_pa
return m_pa
def set_m_pa(new):
global m_pa
m_pa = new
return

def get_sec():
global sec
return sec
def set_sec(new):
global sec
sec = new
return

def get_ssc():
global ssc
return ssc
def set_ssc(new): global ssc
ssc = new
return

def get_sel():
global sel
return sel
def set_sel(new):
global sel
sel = new
return

def get_expsec():
global expsec
return expsec
def set_expsec(new):
global expsec
expsec = new
return

def get_expssc():
global expssc
return expssc
def set_expssc(new):
global expssc
expssc = new
return

def get_expsel():
global expsel
return expsel
def set_expsel(new):
global expsel
expsel = new
return

def get_alpha_c_e():
global alpha_c_e
return alpha_c_e
def set_alpha_c_e(new):
global alpha_c_e
alpha_c_e = new
return

def get_alpha_c_s():
global alpha_c_s
return alpha_c_s
def set_alpha_c_s(new):
global alpha_c_s
alpha_c_s = new
return

def get_alpha_c_se():
global alpha_c_se
return alpha_c_se
def set_alpha_c_se(new):
global alpha_c_se
alpha_c_se = new
return

def get_alpha_n_e():
global alpha_n_e
return alpha_n_e
def set_alpha_n_e(new):
global alpha_n_e
alpha_n_e = new
return

def get_alpha_n_s():
global alpha_n_s
return alpha_n_s
def set_alpha_n_s(new):
global alpha_n_s
alpha_n_s = new
return

def get_alpha_n_se():
global alpha_n_se
return alpha_n_se
def set_alpha_n_se(new):
global alpha_n_se
alpha_n_se = new
return

def get_alpha_p_e():
global alpha_p_e
return alpha_p_e
def set_alpha_p_e(new):
global alpha_p_e
alpha_p_e = new
return

def get_alpha_p_s():
global alpha_p_s
return alpha_p_s
def set_alpha_p_s(new):
global alpha_p_s
alpha_p_s = new
return

def get_alpha_p_se():
global alpha_p_se
return alpha_p_se
def set_alpha_p_se(new):
global alpha_p_se
alpha_p_se = new
return

def get_alpha_m():
global alpha_m
return alpha_m
def set_alpha_m(new):
global alpha_m
alpha_m = new
return

def get_alpha_pa():
global alpha_pa
return alpha_pa
def set_alpha_pa(new):
global alpha_pa
alpha_pa = new
return

def get_use_c_e():
global use_c_e
return use_c_e
def set_use_c_e(new):
global use_c_e
use_c_e = new
return

def get_use_c_s():
global use_c_s
return use_c_s
def set_use_c_s(new):
global use_c_s
use_c_s = new
return

def get_use_c_se():
global use_c_se
return use_c_se
def set_use_c_se(new):
global use_c_se
use_c_se = new
return

def get_use_n_e():
global use_n_e
return use_n_e
def set_use_n_e(new):
global use_n_e
use_n_e = new
return

def get_use_n_s():
global use_n_s
return use_n_s
def set_use_n_s(new):
global use_n_s
use_n_s = new
return

def get_use_n_se():
global use_n_se
return use_n_se
def set_use_n_se(new):
global use_n_se
use_n_se = new
return

def get_use_p_e():
global use_p_e
return use_p_e
def set_use_p_e(new):
global use_p_e
use_p_e = new
return

def get_use_p_s():
global use_p_s
return use_p_s
def set_use_p_s(new):
global use_p_s
use_p_s = new
return

def get_use_p_se():
global use_p_se
return use_p_se
def set_use_p_se(new):
global use_p_se
use_p_se = new
return

def get_use_m():
global use_m
return use_m
def set_use_m(new):
global use_m
use_m = new
return

def get_use_pa():
global use_pa
return use_pa
def set_use_pa(new):
global use_pa
use_pa = new
return

def get_ye0():
global ye0
return ye0
def get_ys0():
global ys0
return ys0
def get_yse0():
global yse0
return yse0

def get_capacity_e():
global capacity_e
return capacity_e
def get_capacity_s():
global capacity_s
return capacity_s
def get_capacity_se():
global capacity_se
return capacity_se

def get_t():
global t
return t

def trinity():
#übergebe Zeit
t = get_t()

#Übergibt Initiale Organismen-Konzentrationen
e_coli = get_ye0()
s_cerevisiae = get_ys0()
s_elongatus = get_yse0()

#Übergibt Initiale Stoff-Konzentrationen (extern)
c = get_itc() #kohlenstoff/zucker
n = get_itn() #Stickstoff
p = get_itp() #Phosphor
m = get_itm() #Melanin
pa = get_itpa() #Phosphit

#Übergibt Initiale Stoff-Konzentrationen (intern)
i_c_e = get_i_c_e() #kohlenstoff/zucker, e.coli
i_c_s = get_i_c_s() #kohlenstoff/zucker, s.cervisae
i_c_se = get_i_c_se() #kohlenstoff/zucker, s.elongatus
i_n_e = get_i_n_e() #Stickstoff, e.coli
i_n_s = get_i_n_s() #Stickstoff, s.cervisae
i_n_se = get_i_n_se() #Stickstoff, s.elongatus
i_p_e = get_i_p_e() #Phosphor, e.coli
i_p_s = get_i_p_s() #Phosphor, s.cervisae
i_p_se = get_i_p_se() #Phosphor, s.elongatus
i_m = get_i_m() #Melanin
i_pa = get_i_pa() #Phosphit

#Übergibt Maximale Stoffkonzentrationen (intern)
m_c_e = get_m_c_e()
m_c_s = get_m_c_s()
m_c_se = get_m_c_se()
m_n_e = get_m_n_e()
m_n_s = get_m_n_s()
m_n_se = get_m_n_se()
m_p_e = get_m_p_e()
m_p_e = get_m_p_e()
m_p_se = get_m_p_se()
m_m = get_m_m()
m_pa = get_m_pa()

#Übergibt Metaboliten Geschwindigkeiten
sec = get_sec()
ssc = get_ssc()
sel = get_sel()
expsec = get_expsec()
expssc = get_expssc()
expsel = get_expsel()

#Übergibt Absorbationsgeschwindingikeiten
alpha_c_e = get_alpha_c_e()
alpha_c_s = get_alpha_c_s()
alpha_c_se = get_alpha_c_se()
alpha_n_e = get_alpha_n_e()
alpha_n_s = get_alpha_n_s()
alpha_n_se = get_alpha_n_se()
alpha_p_e = get_alpha_p_e()
alpha_p_s = get_alpha_p_s()
alpha_p_se = get_alpha_p_se()
alpha_m = get_alpha_m()
alpha_pa = get_alpha_pa()

#Übergibt Metaboliten Nutzung
use_c_e = get_use_c_e()
use_c_s = get_use_c_s()
use_c_se = get_use_c_se()
use_n_e = get_use_n_e()
use_n_s = get_use_n_s()
use_n_se = get_use_n_se()
use_p_e = get_use_p_e()
use_p_s = get_use_p_s()
use_p_se = get_use_p_se()

i = 0
Y_e = [e_coli]
Y_s = [s_cerevisiae]
Y_se = [s_elongatus]

capacity_e = get_capacity_e()
capacity_s = get_capacity_s()
capacity_se = get_capacity_se()

for j in t:

add = min(alpha_m * ((m_m*(Y_e[i]/Y_e[0]) - i_m)/m_m), m)
i_m = i_m+add
m = m - add
ret = min(i_m, sec*(Y_e[i]/Y_e[0]))
n = n + ret
i_m = i_m - ret

add = min(alpha_pa * ((m_pa*(Y_s[i]/Y_s[0])-i_pa)/m_pa), pa)
i_pa = i_pa+add
pa = pa - add
ret = min(i_pa, ssc*(Y_s[i]/Y_s[0]))
p = p + ret
i_pa = i_pa - ret

add = min(alpha_c_e * ((m_c_e*(Y_e[i]/Y_e[0]) - i_c_e)/m_c_e), c)
i_c_e = i_c_e+add
growthf_c_e = min(i_c_e/(use_c_e*(Y_e[i]/Y_e[0])), 1)
i_c_e = i_c_e-growthf_c_e*use_c_e*(Y_e[i]/Y_e[0])
c = c-add

add = min(alpha_c_s * ((m_c_s*(Y_s[i]/Y_s[0]) - i_c_s)/m_c_s), c)
i_c_s = i_c_s+add
growthf_c_s = min(i_c_s/(use_c_s*(Y_s[i]/Y_s[0])), 1)
i_c_s = i_c_s - growthf_c_s*use_c_s*(Y_s[i]/Y_s[0])
c = c-add

add = min(alpha_c_se * ((m_c_se*(Y_se[i]/Y_se[0]) - i_c_se)/m_c_se), c)
i_c_se = i_c_se+add + sel - sel*expsel
growthf_c_se = min(i_c_se/(use_c_se*(Y_se[i]/Y_se[0])), 1)
i_c_se=i_c_se- growthf_c_se*use_c_se*(Y_se[i]/Y_se[0])
c = c-add + sel*expsel*(Y_se[i]/Y_se[0])

add = min(alpha_n_e * ((m_n_e*(Y_e[i]/Y_e[0]) - i_n_e)/m_n_e), n)
i_n_e = i_n_e+add
growthf_n_e = min(i_n_e/(use_n_e*(Y_e[i]/Y_e[0])), 1)
i_n_e = i_n_e-growthf_n_e*use_n_e*(Y_e[i]/Y_e[0])
n = n-add

add = min(alpha_n_s * ((m_n_s*(Y_s[i]/Y_s[0]) - i_n_s)/m_n_s), n)
i_n_s = i_n_s+add
growthf_n_s = min(i_n_s/(use_n_s*(Y_s[i]/Y_s[0])), 1)
i_n_s=i_n_s-growthf_n_s*use_n_s*(Y_s[i]/Y_s[0])
n = n-add

add = min(alpha_n_se * ((m_n_se*(Y_se[i]/Y_se[0]) - i_n_se)/m_n_se), n)
i_n_se = i_n_se+add
growthf_n_se = min(i_n_se/(use_n_se*(Y_se[i]/Y_se[0])), 1)
i_n_se = i_n_se-growthf_n_se*use_n_se*(Y_se[i]/Y_se[0])
n = n-add

add = min(alpha_p_e * ((m_p_e*(Y_e[i]/Y_e[0]) - i_p_e)/m_p_e), p)
i_p_e = i_p_e+add
growthf_p_e = min(i_p_e/(use_p_e*(Y_e[i]/Y_e[0])), 1)
i_p_e=i_p_e-growthf_p_e*use_p_e*(Y_e[i]/Y_e[0])
p = p-add

add = min(alpha_p_s * ((m_p_s*(Y_s[i]/Y_s[0]) - i_p_s)/m_p_s), p)
i_p_s = i_p_s+add
growthf_p_s = min(i_p_s/(use_p_s*(Y_s[i]/Y_s[0])), 1)
i_p_s=i_p_s-growthf_p_s*use_p_s*(Y_s[i]/Y_s[0])
p = p-add

add = min(alpha_p_se * ((m_p_se*(Y_se[i]/Y_se[0]) - i_p_se)/m_p_se), p)
i_p_se = i_p_se+add
growthf_p_se = min(i_p_se/(use_p_se*(Y_se[i]/Y_se[0])), 1)
i_p_se=i_p_se-growthf_p_se*use_p_se*(Y_se[i]/Y_se[0])
p = p-add

if (growthf_c_e == 0 or growthf_n_e == 0 or growthf_p_e == 0):
growth_e = -0.1
else:
growth_e = min(growthf_c_e, growthf_n_e, growthf_p_e)

if (growthf_c_s == 0 or growthf_n_s == 0 or growthf_p_s == 0):
growth_s = -0.1
else:
growth_s = min(growthf_c_s, growthf_n_s, growthf_p_s)

if (growthf_c_se == 0 or growthf_n_se == 0 or growthf_p_se == 0):
growth_se = -0.1
else:
growth_se = min(growthf_c_se, growthf_n_se, growthf_p_se)

Y_e = Y_e + [Y_e[i]+growth_e * ((capacity_e - e_coli) / capacity_e)]
Y_s = Y_s + [Y_s[i]+growth_s * ((capacity_s - s_cerevisiae) / capacity_s)*(1/8)]
Y_se = Y_se + [Y_se[i]+growth_se * ((capacity_se - s_elongatus) / capacity_se)*(1/24)]

e_coli = Y_e[i]
s_cerevisiae = Y_s[i]
s_elongatus = Y_se[i]
i = i+1

Y_e.pop();
Y_s.pop();
Y_se.pop();
return Y_e, Y_s, Y_se;

fig = plt.figure()
ax = plt.axes()

Y = trinity();
for i, species in enumerate(["$\t{E. coli}$", "$\t{S. cerevisiae}$", "$\t{S. elongatus}$"]):
plt.plot(t, Y[i], label=species)
plt.legend(frameon=False, loc='upper right', ncol=1, labelspacing=1)
plt.xlabel("Time in Hours")
plt.ylabel("Dryweight in Mol")
plt.show()

@@ Line 1,504: / Line 1,504: @@
 plt.show()<br>
 </p>
-  				</div>
+</div>
-<button class="tablink" onclick="openCity('Info', this, 'gray')" id="defaultOpen">Info</button>
-<button class="tablink" onclick="openCity('Code', this, 'gray')">Code</button>
 </article>
 </body>

Difference between revisions of "Team:Duesseldorf/Model"

Revision as of 07:56, 8 October 2018

Model

Code