Team:Buenos Aires/Project/Model

From 2012.igem.org

(Difference between revisions)
(Modeling)
(Overview)
 
(203 intermediate revisions not shown)
Line 2: Line 2:
= Modeling a synthetic ecology =
= Modeling a synthetic ecology =
-
To gain insight into the behavior of our crossfeeding design, to understand which  parameters of the system are important, and to study the feasibility of the project, we decided to implement a model of the system. We found this interesting per se, as the model we need is that of a microbial community, an ecology.  
+
To gain insight into the behavior of our crossfeeding design, to understand which  parameters of the system are important, and to study the feasibility of the project, we decided to implement a model of the system. We found this interesting per se, as the model we need is that of a microbial community and ecology.  
The mathematical modeling of ecological system is an active field of research with a very rich and interesting history; from the works of Lotka and Volterra, who modeled predator-prey dynamics with ordinary differential equations, to Robert May’s chaotic logistic maps. In fact a lot of synthetic systems of interacting bacteria created in recent years were used to study these sort of models.
The mathematical modeling of ecological system is an active field of research with a very rich and interesting history; from the works of Lotka and Volterra, who modeled predator-prey dynamics with ordinary differential equations, to Robert May’s chaotic logistic maps. In fact a lot of synthetic systems of interacting bacteria created in recent years were used to study these sort of models.
Line 15: Line 15:
The model has four variables:
The model has four variables:
-
[Nhis-] : the concentration of histidine dependent (Trp producing) yeast cells, in cells per ml
+
*[Nhis-] : the concentration of histidine dependent (Trp producing) yeast cells, in cells per ml
-
[Ntrp-] :  the total amount of tryptophan dependent (His producing) yeast cells, in cell per ml
+
*[Ntrp-] :  the total amount of tryptophan dependent (His producing) yeast cells, in cell per ml
-
[his]: the concentration of histidine in the medium, in molecules per ml
+
*[his]: the concentration of histidine in the medium, in molecules per ml
-
[trp]: the concentration of tryptophan in the medium, in molecules per ml.
+
*[trp]: the concentration of tryptophan in the medium, in molecules per ml.
To build the model we did the following assumptions:
To build the model we did the following assumptions:
-
# The growth rate of each strain depends on the concentration in the medium of the AA they can’t produce. As the concentration of this AA in the medium increases, so does the growth rate of the strain, until it settles at the growth rate observed in optimal conditions (doubling time of 90 min for yeast).
+
* The growth rate of each strain depends on the concentration in the medium of the AA they can’t produce. As the concentration of this AA in the medium increases, so does the growth rate of the strain, until it settles at the growth rate observed in optimal conditions (doubling time of 90 min for yeast).
-
# There is a maximal density of cells the medium can support
+
* There is a maximal density of cells the medium can support (carrying capacity)
-
# Each cell has a fixed probability of dying per time unit
+
* Each cell has a fixed probability of dying per time interval
-
# Each cell releases to the medium the AA it produces at a fixed rate
+
* Each cell releases to the medium the AA it produces at a fixed rate
-
# Each strain only consumes the AA it doesn't produce
+
* Each strain only consumes the AA it doesn't produce
-
# The system reaches steady state.
+
* The system reaches steady state.
 +
[[File:Bsas2012-model1.png]] (model)
-
{| class="wikitable" border="1"
+
In the equation for the temporal evolution of each population we take the growth rate as a Hill function of the amount of AA to which the strain is auxotrophic (it can't produce). The function captures the increase of the growth rate with the relevant AA concentration, until it reaches a plateau at the maximal growth rate, that corresponds to a doubling time of 90 minutes (doubling time = ln(2)/kmax). These functions are of the form
-
|-
+
-
| [[File:Bsas2012-model-latex1.png|620px]]
+
-
|-
+
-
| [[File:Bsas2012-model-latex2.png|620px]]
+
-
|}
+
-
 
+
-
{| class="wikitable" border="1"
+
-
|-
+
-
| [[File:Bsas2012-model-latex4.png|383px]]
+
-
|-
+
-
| [[File:Bsas2012-model-latex3.png|383px]]
+
-
|}
+
-
 
+
-
 
+
-
 
+
-
==== Subsubsection ====
+
-
 
+
-
In the equation for the temporal evolution of each population we take the growth rate as a Hill function of the amount of AA to which the population is auxotrophic. The function represents the stabilization of the rate to a 90 minutes cell division cycle (kmax).
+
-
 
+
-
<math>
+
-
k=  \frac{ k_{max} [AA]^{n}}{ [AA]^{n} + K_{AA}}
+
-
</math>
+
-
 
+
-
KAA is the concentration of AA in the medium where k = 0.5  k max.
+
-
 
+
-
The second term accounts for other factors, not explicit in our model that can limit the concentration of cells in a given volume (capacity) enabling the population to reach equilibrium. Given that our model relies on nutrient deficiency we can't rule out starvation so we added a term related to cell death.
+
-
 
+
-
The terms in the equations for the evolution of each [AA] in the medium are the fluxes of AA entering and leaving the medium. Each rate '''p''' is a measure of the AA produced and exported in the form of a protein by a cell. In our simulations it is defined as a fraction of the calculated maximum of elongation events.
+
-
 
+
-
The second term is the flux leaving the medium and entering the auxotrophic cells. We have assumed that the total amount of AA absorbed, '''d''', during the cell cycle time &tau;  is the required to ''build'' the daughter. One of the hypotheses built into our model is that &tau;  will vary greatly with the concentration of nutrients available. We've used
+
-
 
+
-
<math>
+
-
\tau =k^{-1} = \frac{ [AA]^{n} + K_{AA}}{ k_{max} [AA]^{n}}</math>
+
-
The amount of AA entering the cells that produce it was considered negligible.
+
[[File:bsas2012-modeling-eq1.png|200px]](1)
 +
where '''Kaa''' is the effective concentration of AA at which half maximal growth rate is obtained, and '''l''' is the Hill coefficient, that describes how "steep" the curve is.
-
===== Steady State Solution =====
+
The term (1-(Nhis+Ntrp)/Cc) of the model accounts for other factors, not explicit in our model that can limit the concentration of cells in a given volume (carrying capacity), enabling the population to reach equilibrium. We included a term related to cell death (- death Ntrp), which is makes biologically sense (cell do die) and is critical for a correct behavior of the model.
-
The 4 equations in [] were equaled to zero and the non-linear algebraic system solved using MATEMATICA to find their equilibrium values.
+
The terms in the equations for the evolution of each [AA] in the medium are the fluxes of AA entering and leaving the medium. The first term (Ptrp*Nhis-) is a measure of the AA produced and exported in the form of peptides by a cell, for example the rate of tryptophan secretion by the histidine dependent cell.  
-
Here we change the notation a little to two generic strains '''a''' and '''b''' where population '''Na''' produces amino acid '''a''', etc. We are interested in the value of the cell populations in equilibrium, more importantly in the fraction of each population. Thus a change was made to more convenient variables: the sum of both population and a strain's molar fraction.
+
The second term is the flux leaving the medium and entering the auxotrophic cells. Each cell has a more or less fixed amount of each amino acid, which we call '''d'''. If a cell is to replicate itself and cannot produce and amino acid, it will have to absorb it from the medium. Therefore the flux of AA entering the cell, is '''d''' divided by the doubling time &tau; (the time it takes to "construct" a new cell).
 +
One of the hypotheses built into our model is that &tau; will vary greatly with the concentration of nutrients available. We've used
-
{N_{a}  ;  N_{b}} → { N_{t}  ,  x_{a} }
+
[[File:bsas2012-modeling-eq2.png|200px]](2)
-
Three solutions were found:
+
The amount of AA entering the cells that produce it was considered negligible. This means that a cell that produces Trp doesn't import it from the medium in significant quantities. This is likely to hold when AA concentrations are limiting. 
-
                #        {{N_t→0},
+
=== Steady State Solution ===
-
#{Nt→Cc+(Cc√da √db death)/(√pa √pb),          xa→(√db √pb)/(-√da √pa+√db √pb),
+
The 4 equations in the model were equaled to zero and the non-linear algebraic system solved using [http://www.wolfram.com/mathematica/ Mathematica] to find their equilibrium values.
-
AAa→〖(-(Kb√pa √pb)/(√da √db kmax+√pa √pb))〗^(1/m),AAb→〖(-(Kb√pa √pb)/(√da √db kmax+√pa √pb))〗^(1/n)}
+
Here we change the notation a little bit to two generic strains '''a''' and '''b''' where population '''Na''' produces amino acid '''a''' and requires '''b''', and population '''Nb''' produces '''b''' and requires '''a'''. We are interested in the value of the cell populations in equilibrium, or more importantly the fraction of each population in steady state. Thus a change was made to more convenient variables: the sum of both population and a strain's fraction.
-
#{Nt→Cc-(Cc√da √db death)/(√pa √pb),          xa→(√db √pb)/(√da √pa+√db √pb),
+
Nt = Na + Nb
-
AAa→〖((Kb√pa √pb)/(√da √db kmax-√pa √pb))〗^(1/m),AAb→〖((Kb√pa √pb)/(√da √db kmax-√pa √pb))〗^(1/n)}
+
Xa = Na / Nt
-
We’ve only found 3 fixed points for the system {Nt, Xa, AAa, AAb}. That means that for some solutions it’s not true that the four variables reach a steady state or equilibrium. There are initial conditions for which some of the variables will not stop growing!
+
Only 3 fixed points were found for the system (See [[Team:Buenos_Aires/Project/ModelAdvance#Appendix | Appendix]]). There are some solutions for which it isn't true that the four variables reach a steady state (SS) or equilibrium. There are initial conditions for which the AA concentrations will not stop growing! So caution is required when assuming steady state.
-
*The first solution is clear. The culture will not thrive under certain conditions… to be determined!
+
*In the first solution we found the culture doesn't thrive, for example because one strain was missing or the initial density was to low.
-
*The second one has no physical relevance since it yields negative concentrations for the Amino Acids.
+
*The second one has no biological relevance since it yields negative concentrations for the Amino Acids.
-
*However the last one yields :
+
*The third one is the significant one (see equations below). However for some parameters the concentrations of Nt and AA can be negative, therefore restricting the parameter space in which the model works.
-
N_t→Cc-(Cc √(da ) √db  death)/(√pa  √pb),
 
-
x_a→1/(√((da pa)/(db pb))+1)
+
[[File:bsas2012-modeling-eqsol3.png]](3)
=== Regulation ===
=== Regulation ===
-
Note that the percentage of each population in the community is a function of √((da pa)/(db pb))  
+
Note that the fraction of each strain in the community is a function of the AA secretion rates ('''pa''' and '''pb''') and the amount of AA required to "construct" a cell ('''da''' and '''db''').
-
In particular
+
[[File:bsas2012-modeling-eq4.png]] (4)
-
x_a→1/(√((da pa)/(db pb))+1)
+
For positive values of '''d''' and '''p''' the range of the function is (0; 1), consistent with what we expect for a fraction.
-
For positive values of d and p the range of the function is (0; 1], so far it’s consistent with what we expect:
+
The fraction of each strain in the culture in equilibrium are regulated by the production and export of each AA – represented in the model through '''pa''' and '''pb'''. These fractions are independent of initial conditions! This means that the system auto-regulates itself, as intended. 
-
The regulation of the production and export of AA, allows us to vary parameters pa and pb .This exploration of the parameter space determines the percentages of the strains in the culture in equilibrium; independent of initial conditions!!
+
-
The next figure illustrates how the molar fraction <math> x_{a} = N_{a}/N_{t}</math> varies with the variable &gamma;,  where <math> \frac{p_{a}}{p_{b}}= 10^{\gamma} </math>
+
In fact we only need to control the ratio between these parameters to control the culture composition:
-
 
+
   
-
  Figure 1.
+
[[File:bsas2012-modeling-eq5revised.png]](5)
-
 
+
-
.. shows the results for two values of the ratio da/db; 1 on the left and 10 on the right. The ratio is fixed for two A.A. in our model and determines the range of x_a and therefore percentages of each strain accessible for any selection of biobricks that regulate pa/pb.
+
-
 
+
-
Some basic properties were observed by simplifying the system; we wanted to check that our results were sound.
+
-
 
+
-
Nt=Cc-(Cc √da √db death)/(√pa √pb)  when {pa→p,pb→p,da→d,db→d}
+
-
Nt=Cc  (1-(d death)/p)
+
The next figure illustrates how the fraction Xa varies with the variable &epsilon;.
-
The condition for culture growth is    =  p /  ( d  death  )  >  1.  This ratio between two rates determines whether the system is able to produce and absorbed the nutrients needed for the cell cycle before its death. If is not true then the total amount of cells is negative, ie the culture dies.
+
[[File:bsas2012-modeling-fig1revised.png]]
 +
Figure 1. Fraction Xa as a fuction of the ratio of protein export for values of D=0.1,1,10 (in green, purple and  red).
-
  Figure 2.
+
To programme the percentage we need a second set of biobricks that can sense an external stimulus and transduce this signal  to modify &epsilon;. The range of percentages we can control is determined by the range of &epsilon; accessible with this second device and the ratio D that is set for any two A.A.
=== Parameters ===
=== Parameters ===
-
The parameter estimation process is shown in xx section(linked)
+
The parameter estimation process is detailed in [[Team:Buenos_Aires/Project/Model#Parameter_selection | parameter selection]]
{| class="wikitable"
{| class="wikitable"
Line 179: Line 142:
|-
|-
|style="color:white; background-color: purple;"|d trp
|style="color:white; background-color: purple;"|d trp
-
|2.630e7
+
|2.630e7  
|AA/cell
|AA/cell
|}
|}
-
   Table 1.
+
   Table 1. Model parameters used for the simulations
-
== Model Analysis ==
+
=== Parameter selection ===
-
There was a major difficulty in the analysis of the model because of the dimension of the problem - 4 variables. The analytical results yielded by the usual tools (nullclines, linearization) could not be simplified so that all the relevant properties of the system could be understood in terms of the parameters.
+
-
We’ve relied on numerical simulations to explore possible behaviors.
+
We estimatated values for all the parameter in the model, doing dedicated experiments or using values from the literature. This allowed to check the feasibility of the system.
-
===Simulations===
+
-
To run the simulations all that is left is to choose values {p,d, k death} from an interval around the estimated value P discussed before- since this is merely a framewok we have taken values p  from a log scale from 10^(-3) to 10. Note that epsilon controls the fraction of each population, in fact according to the equation for x_a should only depend on epsilon.  
+
-
==== Die or thrive?====
+
*The '''Kaa''' found in the equations are related to the concentration of AAs in the medium required to reach half the maximum rate of growth. Curves of OD600 vs [AA](t=0) were measured experimentally for each amino acid and the data fitted to a Hill equation (6) using MATLAB’s TOOLBOX: Curve Fitting Tool.
-
The first test is under which initial conditions do the populations thrive or die?  That is given values p from the interval which is the minimum initial concentration of cells required to do so. The model accounts for cell death, we’ll make a stronger restraint to the allowed values; requiring that the systems has a lag time smaller than 70 hrs.
+
 +
[[File:Bsas2012-modeling-eqfit.png|200px]](6)
-
First we’ve checked the validity of:  
+
This data from [[Team:Buenos_Aires/Results/Strains#Growth dependence on the Trp and His concentrations | experimental results]] and the best fit obtained for each are shown below in Figure 2.
-
  Nt=Cc ( 1 - ( √da √db  death)/(√pa √pb))
+
[[File:Bsas2012-modeling-fig_aux1.png]]
 +
[[File:Bsas2012-modeling-fig_aux2.png]]
 +
  Figure 2. Single strain culture density after over night growth vs the initial concentration of AA in the medium (dilutions 1:X).
 +
  The best fit to the data is also shown.
-
We’ve also noticed that the culture’s survival is dependent on the Ks even though is not explicit in the formula above; this is reasonable since it introduces another time scale. Not only have the cells to produce and absorb the required amounts of AA before they die, but their growth is restricted by the concentration of AA in the medium.  The cells must produce amino acids not in the order of d but K in the time frame to insure their prosperity.  
+
The values taken from the fits are  
-
  Figure 3. Evolution in time of two cultures that only differ in the value of the Ks for n_ hill =1 (left)and n_hill >1 (right).
+
His- :
-
  When n_ hill =1 only the one with lower Ks reaches a non-trivial steady state (ss). The other slowly dies. However for n_ hill =1.5
+
      K_dil =     0.2255 
-
   there is no difference between the cases with high and low Ks. Note that it reaches the ss faster as well.
+
      n =         1.52
 +
   R-square: 0.997  
-
The survival of an auxotrophic community is dependent on the concentration of amino acids in the medium required to reach half the maximum rate of growth (Ks) and the Hill coefficient as shown clearly on Figure 3.  
+
Trp-  :
-
So when creating an auxotrophy, pick a nutrient with an appropriate Hill coefficient and K. Check the range and stability of your model estimating the values  of the rest of the parameters.
+
      K_dil =    0.0469
 +
      n =          1.895
 +
  R-square: 0.998
-
Our case is shown in figure 3 with both hill coefficients are > 1 and we haven’t found the same variability in the behavior. It’s insensitive to change in the same range! And it’s faster!
+
 +
The results were then converted to the units chosen for the simulations.
-
===Parameter Space and Solutions ===
+
K = K_dil * [AA 1x] * Navog / (Molar mass AA)
-
The analytical solutions found for the system do not exhaust all possible ones - there are  fewer than the dimension of the problem. While running the simulations it became clear that there were cases where each auxotrophic strain reached a steady state but the fraction is not consistent with eq. X.
+
where [AA 1x] = 0.02 mg/ml is the concentration of the AA (His or Trp)  in the 1x medium and Navog is Avogadro's number. We get
-
Example
+
*K trp= .0469* 5.88e16 AA / ml =  2.76e15 molecules/ml
 +
*K his= 0.2255* 7.8e16  AA / ml = 1.76e16 molecules /ml
-
  Figure 4. For values p=0.1 and 0.5 the AA reach ss and the x_a tends to the predict one by the fomula;
+
The competition within a strain and with the other were taken as equal. The system's general ''carrying capacity'' considered in the model is the one often used here, in a lab that works with these yeast strains:
-
  fa =    0.6872. For p=1 the relationship found for x_a does not hold.
+
-
We can show this effect by varying the value of p.
+
**Cc= 3e7  cell/ml.
-
Next we see for another set of parameters how as p decreases the fraction x_a in ss tends to the theoretical one x_a=0.9076 (left). It’s also important to consider that the relaxation time will depend on p for those values where the fraction is the expected one (right). There must be a compromise; high enough that our time scale is a reasonable one, low enough that the formula and regulation of the community we seek holds true.
+
-
  Figure 5.
 
-
How do we explain this?
+
*The death rate was taken between 3 and 7 days.
-
For certain parameters [AA] does not reach a limit but continues to grow; saturating the medium. The growth rate then tends to the theoretical 90 min and the AA dependence disappears from the equations for each population:
 
-
  Eqs.
+
*The rate of “production and export” was given an upper bound value (P_MAX) of 1% of the maximum estimated number of protein elongation events (peptidil transferase reactions) for a yeast cell per hour.  That is, if all the biosynthetic capacity of the cell was used to create the AA rich exportation peptide, the export rate would equal P_MAX. Of course this is not possible because the cell has to do many other things, therefore we considered 1% of P_MAX as a reasonable upper bound for '''p'''.
-
The new system is degenerate, which in this case means solutions are set by the initial conditions and not by parameters subject of regulation. Thus the formula for the fraction x_a does not ALWAYS hold. It’d be great if the subset of parameter for which it does could be analytically determined. It’s also important to determine the regions in the parameter space that separate these different types of solutions
+
From [http://www.biomedcentral.com/1752-0509/2/87| von der Haar 2008]  we get an estimate for the total number of elongation events (peptidil transferase reactions):  6e6  1 / cell sec [ ]
-
However, falling that, it was done numerically.
+
* P_MAX = 2.16e8 1/ cell hour
-
=== Relaxation time===
+
These parameters will be regulated in the simulation by &epsilon;  and '''p''', to alter the fraction between populations.
-
Not only are we interested in predicting the fraction of each strain in the culture, we want it within a time frame.
+
*p trp = '''p'''* '''&epsilon;'''*2.16e8  AA / cell hour
-
Creating a mesh we explore the parameter space {p, epsilon} for different initial conditions to visualize the subspace where the model prediction and the numerical result of the simulations --differ in an amount lower than a given error; that means that the difference  stays within a band-- given by the error. We want this to happen before 50 hours.
+
-
We’ve taken 5% to 10% of the value as error since is not necessary to ask for a precision that can’t be later measured experimentally.
+
-
To explore the parameter space we have set the values pa  and pb for different initial conditions. We plot whether the theoretical results agree with the numerical for the fraction x_a - within an error of X %. Thus we have the next binary graph where (red)== in agreement,   (blue) == not in agreement. The results to the left are for a run time of 250 hrs.
+
*p his = '''p'''*2.16e8   AA / cell hour
-
  Figure 6.
+
where '''p''' < 1 controls the fraction of the P_MAX value and &epsilon; the ratio between the export of each amino acid.
-
=== Oscillations===
 
-
One interesting effect appears if the initial concentration of AA in the medium is of the order of the Ks.  Damp oscillations were observed!!
+
'''d''' is the total number of amino acids in a yeast cell –same as the ones needed to create a daughter -per cell:
-
+
-
  Figure 7.
+
-
The next figure present the oscillations for different initial dilutions. We found that the more diluded the initial culture, the more damped the oscillations.
+
'''d''' = # A.A. per cell = (mass of protein per yeast cell ) * relative abundance of AA * Navog / AA's molar mass.
-
  Figure 8.
+
*d trp= 2.630e7 AA / cell
-
== Parameter selection ==
+
*d his= 6.348e8 AA / cell
-
We've estimatated values for all the parameter in the model to check is viability.
+
==Numerical Simulations==
 +
:Initially we relied on numerical simulations (NS) performed with [http://www.mathworks.com/products/matlab/| MATLAB] to explore possible behaviors of the model, ODEs rarely have solutions that can be expressed in a closed form.
-
*The Ks found in the equations are a measure of the concentration of A.A. in the medium required to reach half the maximum rate of growth. Curves OD vs [AA] were measured experimentally for each amino acid; K was taken as the concentration for which the OD falls to half the saturation value.
+
:To run the simulations all that is left is to choose values {p, &epsilon;, initial conditions}.
-
Figure  aux his and  aux trp.
+
* Since this is merely a framewok we have taken values '''p''' from a log scale from -3 to 1.
 +
* '''&epsilon;''' was ranged from 0.0001 to 4; which varies the fraction of each population from 10% to 90%.
 +
* Cells initial concentration (i.c.c.): dilutions from 1:10000 to 10x of the carrying capacity (Cc).
 +
* Amino acids (AA) initial concentration: null, unless stated otherwise.
-
The values were taken from a fit to f(x) = ODmax*x^n/(x^n + K_dil) with MATLAB’s TOOLBOX: Curve Fitting Tool
+
== Die or thrive?==
 +
:Some basic properties were observed by simplifying the system; we wanted to check that our model was sound. We considered a "symmetric interaction" in which both secretion rates and amino acid requirement were equivalent.
-
His-  :
+
[[File:bsas2012-modeling-eq6.png]](7)
-
      K_dil =    0.3279 
+
-
      ODmax =      62.53 
+
-
      n =          1.243
+
-
  R-square: 0.9376 
+
-
Trp-  :
+
where we can define
-
      K_dil =    0.1770
+
-
      ODmax =      57.46 
+
-
      n =          1.636 
+
-
  R-square: 0.9905
+
-
+
[[File:bsas2012-modeling-eq7.png]] (8)
-
The results were then converted to the units chosen for the simulations.
+
:This parameter &lambda; has an intuitive interpretation. The ratio '''d'''/'''p''' is the time it takes a cell to export enough AA for a cell of the other strain to build a daughter. On the other hand, 1/death is the average life span of a cell. Therefore &lambda; is the ratio between the "life span" of a cell and the time it takes to export sufficient AA to build a cell, thus &lambda; represents the amount of cells that can be constructed with the AA secreted by a cell in its lifetime. If this value is grated than one (&lambda;  &gt;  1) the culture grows, otherwise it dies (This is completely analogous to the "force of infection" as defined in epidemiology).    
-
K = K_dil * [AA 1x] # Avog / (Molar mass AA), where [AA 1x] = 0.02 mg/ml
+
:Keeping '''d''' and '''death''' constant, the total steady state number of cells in the culture depends on the secretion rate '''p''' as shown in the following figure 3. There is a threshold value p given by &lambda;=1; with lower values the culture dies.
-
**K trp= 0.1770 * 5.88e16 AA / ml.
+
[[File:bsas2012-modeling-fig3revised.png]]
-
**K his= 0.3279 * 7.8e16  AA / ml.
+
Figure 3. Total number of cells in the mix vs the strengh '''p''' of the production and export of AAs for a set of parameters Cc,d, death.  
-
*The ''capacity'' of the system taken in the model is the one often used in the lab that works with these yeast strains:
+
:Another condition is related to extra-cellular concentration of amino acids. Looking at the equations for the evolution of the AA in the model, we can identify the parameter &delta; as follows
 +
[[File:bsas2012-modeling-eq8.png]](9)
-
**Cc=  3e7  #cell/ml.
+
: As before, the ratio '''&radic;(da db)/&radic;(pa pb)''' is the average time required for a cell to export enough AA to construct another cell. '''1/kmax''' is the time it takes a cell to replicate in optimal growth conditions. &delta; can be interpret as the amount of cells that can be made with the material secreted a by a single cell before it divides (in optimal conditions). If &delta; &gt; 1 the amount of nutrients produced exceeds the consumption, the medium gets saturated with AA and the regulation fails. So the system is auto-regulated only if &delta; &lt; 1.
 +
:Combining these two conditions (&lambda;  &gt;  1 and &delta; &lt; 1), we conclude that the production rate has to be  in a defined region for the system to work. To low and the culture dies out, to big and it gets out of control.
-
*The death rate was taken between 3 and 7 days.
+
== Parameter Space and Solutions ==
 +
:Uniting the conditions for &lambda; and &delta; we see that in general '''p''' and &epsilon;  must be bound for our solution to make sense (i.e. all concentrations &gt; 0):
 +
[[File:bsas2012-modeling-eq14.png]](10)
-
*The rate of “production and export” was given a top value; 1%  of  the maximum estimated number of elongation events for a yeast cell per hour. Our model does not explicitly include time as a variable, that means that the production is switched ON all the time.
 
-
Elongation events (peptidil transferase reactions)6e6 1 / cell sec [ref]
+
[[File:bsas2012-modeling-fig3.png]]
 +
  Figure 4. Regions in the parameter space that present different types of solutions. Only in Region III the system works as intended.
-
These parameters will be regulated experimentally ('''&epsilon;''' , ''p'') to alter the fraction between populations.
 
-
**p trp = ''p''*6e4*3600                    aa / cell hour
+
:Now, let's classify these regions numerically to see whether the ideas we just presented are supported. Taking random values from each region we found that all four variables are always positive, but each region is associated with a different kind of solution. The conditions shape the behavior and not the existence of a solution.
-
**p his = ''p''* '''&epsilon;'''*6e4*3600  aa / cell hour
+
:Taking AA(t=0)=0:
-
where ''p'' < 1 controls the fraction of the top value and '''&epsilon;''' the ratio between  the export of each amino acid.
+
::I. The culture doesn't grow, it decays with varying velocities for any initial concentration of cells (i.c.c.). This is consistent with solution 1, where Nt= 0 as t &rarr; &infin;. This is the &lambda; &lt; 1 case.
-
*'''d''' is the total number of amino acids in a yeast cell –same as the ones needed to create a daughter -per cell:
+
::II. The extracellular AA concentration keeps growing and the medium gets saturated. The culture reaches equilibrium for every i.c.c. no matter how low, but there is no regulation of the composition of the culture. This is the &delta; &gt; 1 case.
-
'''d''' = # A.A. per cell = (mass of protein per yeast cell [ref]) * relative abundance of AA * # Avog / AA's molar mass.
+
::The growth rate tends to the theoretical 90 min and the AA dependence disappears from the equations for each population. Remarkably the fraction of each strain is still independent of the i.c.c.
-
**d trp= 6.348e8  AA / cell 
+
::III. The four variables reach SS. The mole fraction is solely dependent on &epsilon; according to the formula on equation (5). The total population of the community (Nt) is the same as in Region II.
-
**d his= 2.630e7  AA / cell
+
::'''This is the region where the auxotrophy leads to regulation of the community.''' However the culture doesn’t grow for all i.c.c. – will get back to this point shortly.
-
= References =
+
===Conclusions===
 +
::Some conclusions we've drawn from these numerical simulations:
 +
 
 +
# We had postulated that the solution doesn't hold in region II and expected a dramatic change in the system’s behavior. However there is a smooth transition between regions II and III. For a fixed &epsilon; we can go from region III to region II by increasing '''p'''. Once we cross the threshold we notice that the total population of cells (Nt) doeen't vary from one region to the next. Though the AA continues to grow over time (and accumulate) it does so slowly; more so than we'd expect for a vertical asymptote.
 +
# More so, the formula (5) is a good approximation for the fraction of each strain in the region II close to the limit with region III; the difference between the two increases gradually as '''p''' gets further away from this threshold value. The formula is not true in II; thus regulation fails.
 +
# There is no dramatic change caused by the asymptote, in fact the two regions can't be distinguished experimentally through a strain's fraction measurement in each side of the frontier; so small is the difference between the values across.
 +
# Taking into account that the threshold is based on estimations; perhaps is safer to choose points well inside Region III to ensure proper regulation. This has a cost, for a fixed percentage the production and export of AAs is now lower and the dynamic of the culture slower.
 +
 
 +
:The original purpose of this analysis was to undertand '''under which initial conditions do the populations thrive or die'''. This has to be analyzed in terms of Figure 4, not merely the conditions for &lambda;.
 +
 
 +
:We noticed that the culture's survival in region III also depends on the values of '''Kaa''', even though is not explicit in any formula so far. Remember that '''Kaa''' is the concentration of an AA in the medium required for half maximal growth rate. For a fixed value of '''Kaa''' there is a critical initial concentration of cells below which the culture doesn't prosper. This effect is important when the percentage desired is extreme (close 100%), where the '''p''' value required for the strain gets really low.
 +
 +
:This is reasonable considering that another time scale is introduced. Not only have the cells to produce and absorb the required amounts of AA before they die, they also need to modify the AA concentration of the medium. The cells must produce enough amino acids so that the concentration of AA approximates '''Kaa'', before the original cells die out. If the initial cell density is to low, this won't happen.
 +
 
 +
 
 +
== Lag Time ==
 +
 
 +
{|
 +
|Let's finally look at some numerical simulations to see how the system evolves in time. We present a sample simulation with parameters from Region III to the rigth.  There's a new characeristic that our previous analysis didn't show: the system reaches the desired steady state, however there will be significant lag phase.
 +
 
 +
We found numerically that it's related to time it takes for the Amino Acids in the medium to reach a certain concentration; one similar to the parameters '''Kaa''' and '''Kbb''' respectively. To this effect we note that it's possible to reduce this ''lag time'' by
 +
 
 +
*Using higher initial concentration for each population.
 +
 
 +
*Lowering parameters  '''Kaa''' and '''Kbb'''.
 +
 
 +
Given the relationship between these parameters and the AA absortion rates, we deviced a way to facilitate absortion. 
 +
We decided to work with
 +
 
 +
'''Troyan Peptides'''
 +
|
 +
[[File:timeevol.jpg|350px]]
 +
|-
 +
|This lag time was quantified using the slowest strain's third cell division as a lag time.
 +
 
 +
The unit of the color scale is hours, it's clear than for lower initial concentrations and export rates
 +
|
 +
[[File:BsAs2012_11LagK.jpg.jpg|350px]]
 +
|
 +
|}

Latest revision as of 03:13, 27 October 2012

Contents

Modeling a synthetic ecology

To gain insight into the behavior of our crossfeeding design, to understand which parameters of the system are important, and to study the feasibility of the project, we decided to implement a model of the system. We found this interesting per se, as the model we need is that of a microbial community and ecology.

The mathematical modeling of ecological system is an active field of research with a very rich and interesting history; from the works of Lotka and Volterra, who modeled predator-prey dynamics with ordinary differential equations, to Robert May’s chaotic logistic maps. In fact a lot of synthetic systems of interacting bacteria created in recent years were used to study these sort of models.

The crossfeeding model

In the crossfeeding design, each strain produces and releases to the medium an aminoacid the other strain(s) need to grow. The growth rate of each strain depends on the aminoacid (AA) concentration of those AA it can’t produce. In turn these concentrations depend on the abundance of the other strains. Therefore the growth of the strains is interdependent.

For simplicity and to be consistent with the experimental work, we decided to model two interacting strains, one that produces tryptophan (Trp) but requires histidine (His), and another that produces His but requires Trp.

We decided to use ordinary differential equations to model the system. This is a normal approach when studying population dynamics of this kind.

The model has four variables:

  • [Nhis-] : the concentration of histidine dependent (Trp producing) yeast cells, in cells per ml
  • [Ntrp-] : the total amount of tryptophan dependent (His producing) yeast cells, in cell per ml
  • [his]: the concentration of histidine in the medium, in molecules per ml
  • [trp]: the concentration of tryptophan in the medium, in molecules per ml.

To build the model we did the following assumptions:

  • The growth rate of each strain depends on the concentration in the medium of the AA they can’t produce. As the concentration of this AA in the medium increases, so does the growth rate of the strain, until it settles at the growth rate observed in optimal conditions (doubling time of 90 min for yeast).
  • There is a maximal density of cells the medium can support (carrying capacity)
  • Each cell has a fixed probability of dying per time interval
  • Each cell releases to the medium the AA it produces at a fixed rate
  • Each strain only consumes the AA it doesn't produce
  • The system reaches steady state.

Bsas2012-model1.png (model)

In the equation for the temporal evolution of each population we take the growth rate as a Hill function of the amount of AA to which the strain is auxotrophic (it can't produce). The function captures the increase of the growth rate with the relevant AA concentration, until it reaches a plateau at the maximal growth rate, that corresponds to a doubling time of 90 minutes (doubling time = ln(2)/kmax). These functions are of the form

Bsas2012-modeling-eq1.png(1)

where Kaa is the effective concentration of AA at which half maximal growth rate is obtained, and l is the Hill coefficient, that describes how "steep" the curve is.

The term (1-(Nhis+Ntrp)/Cc) of the model accounts for other factors, not explicit in our model that can limit the concentration of cells in a given volume (carrying capacity), enabling the population to reach equilibrium. We included a term related to cell death (- death Ntrp), which is makes biologically sense (cell do die) and is critical for a correct behavior of the model.

The terms in the equations for the evolution of each [AA] in the medium are the fluxes of AA entering and leaving the medium. The first term (Ptrp*Nhis-) is a measure of the AA produced and exported in the form of peptides by a cell, for example the rate of tryptophan secretion by the histidine dependent cell.

The second term is the flux leaving the medium and entering the auxotrophic cells. Each cell has a more or less fixed amount of each amino acid, which we call d. If a cell is to replicate itself and cannot produce and amino acid, it will have to absorb it from the medium. Therefore the flux of AA entering the cell, is d divided by the doubling time τ (the time it takes to "construct" a new cell). One of the hypotheses built into our model is that τ will vary greatly with the concentration of nutrients available. We've used

Bsas2012-modeling-eq2.png(2)

The amount of AA entering the cells that produce it was considered negligible. This means that a cell that produces Trp doesn't import it from the medium in significant quantities. This is likely to hold when AA concentrations are limiting.

Steady State Solution

The 4 equations in the model were equaled to zero and the non-linear algebraic system solved using [http://www.wolfram.com/mathematica/ Mathematica] to find their equilibrium values.

Here we change the notation a little bit to two generic strains a and b where population Na produces amino acid a and requires b, and population Nb produces b and requires a. We are interested in the value of the cell populations in equilibrium, or more importantly the fraction of each population in steady state. Thus a change was made to more convenient variables: the sum of both population and a strain's fraction.

Nt = Na + Nb

Xa = Na / Nt

Only 3 fixed points were found for the system (See Appendix). There are some solutions for which it isn't true that the four variables reach a steady state (SS) or equilibrium. There are initial conditions for which the AA concentrations will not stop growing! So caution is required when assuming steady state.

  • In the first solution we found the culture doesn't thrive, for example because one strain was missing or the initial density was to low.
  • The second one has no biological relevance since it yields negative concentrations for the Amino Acids.
  • The third one is the significant one (see equations below). However for some parameters the concentrations of Nt and AA can be negative, therefore restricting the parameter space in which the model works.


Bsas2012-modeling-eqsol3.png(3)

Regulation

Note that the fraction of each strain in the community is a function of the AA secretion rates (pa and pb) and the amount of AA required to "construct" a cell (da and db).

Bsas2012-modeling-eq4.png (4)

For positive values of d and p the range of the function is (0; 1), consistent with what we expect for a fraction.

The fraction of each strain in the culture in equilibrium are regulated by the production and export of each AA – represented in the model through pa and pb. These fractions are independent of initial conditions! This means that the system auto-regulates itself, as intended.

In fact we only need to control the ratio between these parameters to control the culture composition:

Bsas2012-modeling-eq5revised.png(5)

The next figure illustrates how the fraction Xa varies with the variable ε.

Bsas2012-modeling-fig1revised.png

Figure 1. Fraction Xa as a fuction of the ratio of protein export for values of D=0.1,1,10 (in green, purple and   red).

To programme the percentage we need a second set of biobricks that can sense an external stimulus and transduce this signal to modify ε. The range of percentages we can control is determined by the range of ε accessible with this second device and the ratio D that is set for any two A.A.

Parameters

The parameter estimation process is detailed in parameter selection

Parameters selected
Value Units
kmax 0.4261 1/hr
K his 2.588e16 AA/ml
K trp 1.041e16 AA/ml
n his 1.243
n trp 1.636
Cc 3.0 10^7 cell/ml
death 3 to 7 days
p his p 2.16e8 AA/ cell hr
p trp p ε 2.16e8 AA/ cell hr
d his 6.348e8 AA/cell
d trp 2.630e7 AA/cell
 Table 1. Model parameters used for the simulations

Parameter selection

We estimatated values for all the parameter in the model, doing dedicated experiments or using values from the literature. This allowed to check the feasibility of the system.

  • The Kaa found in the equations are related to the concentration of AAs in the medium required to reach half the maximum rate of growth. Curves of OD600 vs [AA](t=0) were measured experimentally for each amino acid and the data fitted to a Hill equation (6) using MATLAB’s TOOLBOX: Curve Fitting Tool.

Bsas2012-modeling-eqfit.png(6)

This data from experimental results and the best fit obtained for each are shown below in Figure 2.

Bsas2012-modeling-fig aux1.png Bsas2012-modeling-fig aux2.png

Figure 2. Single strain culture density after over night growth vs the initial concentration of AA in the medium (dilutions 1:X).
The best fit to the data is also shown.

The values taken from the fits are

His-  :

      K_dil =     0.2255  
      n =          1.52
 R-square: 0.997  

Trp-  :

      K_dil =     0.0469
      n =          1.895 
 R-square: 0.998


The results were then converted to the units chosen for the simulations.

K = K_dil * [AA 1x] * Navog / (Molar mass AA)

where [AA 1x] = 0.02 mg/ml is the concentration of the AA (His or Trp) in the 1x medium and Navog is Avogadro's number. We get

  • K trp= .0469* 5.88e16 AA / ml = 2.76e15 molecules/ml
  • K his= 0.2255* 7.8e16 AA / ml = 1.76e16 molecules /ml


The competition within a strain and with the other were taken as equal. The system's general carrying capacity considered in the model is the one often used here, in a lab that works with these yeast strains:

    • Cc= 3e7 cell/ml.


  • The death rate was taken between 3 and 7 days.


  • The rate of “production and export” was given an upper bound value (P_MAX) of 1% of the maximum estimated number of protein elongation events (peptidil transferase reactions) for a yeast cell per hour. That is, if all the biosynthetic capacity of the cell was used to create the AA rich exportation peptide, the export rate would equal P_MAX. Of course this is not possible because the cell has to do many other things, therefore we considered 1% of P_MAX as a reasonable upper bound for p.

From [http://www.biomedcentral.com/1752-0509/2/87| von der Haar 2008] we get an estimate for the total number of elongation events (peptidil transferase reactions): 6e6 1 / cell sec [ ]

  • P_MAX = 2.16e8 1/ cell hour

These parameters will be regulated in the simulation by ε and p, to alter the fraction between populations.

  • p trp = p* ε*2.16e8 AA / cell hour
  • p his = p*2.16e8 AA / cell hour

where p < 1 controls the fraction of the P_MAX value and ε the ratio between the export of each amino acid.


d is the total number of amino acids in a yeast cell –same as the ones needed to create a daughter -per cell:

d = # A.A. per cell = (mass of protein per yeast cell ) * relative abundance of AA * Navog / AA's molar mass.

  • d trp= 2.630e7 AA / cell
  • d his= 6.348e8 AA / cell

Numerical Simulations

Initially we relied on numerical simulations (NS) performed with [http://www.mathworks.com/products/matlab/| MATLAB] to explore possible behaviors of the model, ODEs rarely have solutions that can be expressed in a closed form.
To run the simulations all that is left is to choose values {p, ε, initial conditions}.
  • Since this is merely a framewok we have taken values p from a log scale from -3 to 1.
  • ε was ranged from 0.0001 to 4; which varies the fraction of each population from 10% to 90%.
  • Cells initial concentration (i.c.c.): dilutions from 1:10000 to 10x of the carrying capacity (Cc).
  • Amino acids (AA) initial concentration: null, unless stated otherwise.

Die or thrive?

Some basic properties were observed by simplifying the system; we wanted to check that our model was sound. We considered a "symmetric interaction" in which both secretion rates and amino acid requirement were equivalent.

Bsas2012-modeling-eq6.png(7)

where we can define

Bsas2012-modeling-eq7.png (8)

This parameter λ has an intuitive interpretation. The ratio d/p is the time it takes a cell to export enough AA for a cell of the other strain to build a daughter. On the other hand, 1/death is the average life span of a cell. Therefore λ is the ratio between the "life span" of a cell and the time it takes to export sufficient AA to build a cell, thus λ represents the amount of cells that can be constructed with the AA secreted by a cell in its lifetime. If this value is grated than one (λ > 1) the culture grows, otherwise it dies (This is completely analogous to the "force of infection" as defined in epidemiology).
Keeping d and death constant, the total steady state number of cells in the culture depends on the secretion rate p as shown in the following figure 3. There is a threshold value p given by λ=1; with lower values the culture dies.

Bsas2012-modeling-fig3revised.png

Figure 3. Total number of cells in the mix vs the strengh p of the production and export of AAs for a set of parameters Cc,d, death. 


Another condition is related to extra-cellular concentration of amino acids. Looking at the equations for the evolution of the AA in the model, we can identify the parameter δ as follows

Bsas2012-modeling-eq8.png(9)

As before, the ratio √(da db)/√(pa pb) is the average time required for a cell to export enough AA to construct another cell. 1/kmax is the time it takes a cell to replicate in optimal growth conditions. δ can be interpret as the amount of cells that can be made with the material secreted a by a single cell before it divides (in optimal conditions). If δ > 1 the amount of nutrients produced exceeds the consumption, the medium gets saturated with AA and the regulation fails. So the system is auto-regulated only if δ < 1.
Combining these two conditions (λ > 1 and δ < 1), we conclude that the production rate has to be in a defined region for the system to work. To low and the culture dies out, to big and it gets out of control.

Parameter Space and Solutions

Uniting the conditions for λ and δ we see that in general p and ε must be bound for our solution to make sense (i.e. all concentrations > 0):

Bsas2012-modeling-eq14.png(10)


Bsas2012-modeling-fig3.png

Figure 4.  Regions in the parameter space that present different types of solutions. Only in Region III the system works as intended.


Now, let's classify these regions numerically to see whether the ideas we just presented are supported. Taking random values from each region we found that all four variables are always positive, but each region is associated with a different kind of solution. The conditions shape the behavior and not the existence of a solution.
Taking AA(t=0)=0:
I. The culture doesn't grow, it decays with varying velocities for any initial concentration of cells (i.c.c.). This is consistent with solution 1, where Nt= 0 as t → ∞. This is the λ < 1 case.


II. The extracellular AA concentration keeps growing and the medium gets saturated. The culture reaches equilibrium for every i.c.c. no matter how low, but there is no regulation of the composition of the culture. This is the δ > 1 case.
The growth rate tends to the theoretical 90 min and the AA dependence disappears from the equations for each population. Remarkably the fraction of each strain is still independent of the i.c.c.


III. The four variables reach SS. The mole fraction is solely dependent on ε according to the formula on equation (5). The total population of the community (Nt) is the same as in Region II.
This is the region where the auxotrophy leads to regulation of the community. However the culture doesn’t grow for all i.c.c. – will get back to this point shortly.

Conclusions

Some conclusions we've drawn from these numerical simulations:
  1. We had postulated that the solution doesn't hold in region II and expected a dramatic change in the system’s behavior. However there is a smooth transition between regions II and III. For a fixed ε we can go from region III to region II by increasing p. Once we cross the threshold we notice that the total population of cells (Nt) doeen't vary from one region to the next. Though the AA continues to grow over time (and accumulate) it does so slowly; more so than we'd expect for a vertical asymptote.
  2. More so, the formula (5) is a good approximation for the fraction of each strain in the region II close to the limit with region III; the difference between the two increases gradually as p gets further away from this threshold value. The formula is not true in II; thus regulation fails.
  3. There is no dramatic change caused by the asymptote, in fact the two regions can't be distinguished experimentally through a strain's fraction measurement in each side of the frontier; so small is the difference between the values across.
  4. Taking into account that the threshold is based on estimations; perhaps is safer to choose points well inside Region III to ensure proper regulation. This has a cost, for a fixed percentage the production and export of AAs is now lower and the dynamic of the culture slower.
The original purpose of this analysis was to undertand under which initial conditions do the populations thrive or die. This has to be analyzed in terms of Figure 4, not merely the conditions for λ.
We noticed that the culture's survival in region III also depends on the values of Kaa, even though is not explicit in any formula so far. Remember that Kaa is the concentration of an AA in the medium required for half maximal growth rate. For a fixed value of Kaa there is a critical initial concentration of cells below which the culture doesn't prosper. This effect is important when the percentage desired is extreme (close 100%), where the p value required for the strain gets really low.
This is reasonable considering that another time scale is introduced. Not only have the cells to produce and absorb the required amounts of AA before they die, they also need to modify the AA concentration of the medium. The cells must produce enough amino acids so that the concentration of AA approximates 'Kaa, before the original cells die out. If the initial cell density is to low, this won't happen.


Lag Time

Let's finally look at some numerical simulations to see how the system evolves in time. We present a sample simulation with parameters from Region III to the rigth. There's a new characeristic that our previous analysis didn't show: the system reaches the desired steady state, however there will be significant lag phase.

We found numerically that it's related to time it takes for the Amino Acids in the medium to reach a certain concentration; one similar to the parameters Kaa and Kbb respectively. To this effect we note that it's possible to reduce this lag time by

  • Using higher initial concentration for each population.
  • Lowering parameters Kaa and Kbb.

Given the relationship between these parameters and the AA absortion rates, we deviced a way to facilitate absortion. We decided to work with

Troyan Peptides

Timeevol.jpg

This lag time was quantified using the slowest strain's third cell division as a lag time.

The unit of the color scale is hours, it's clear than for lower initial concentrations and export rates

BsAs2012 11LagK.jpg.jpg