"
Team:Buenos Aires/Project/Model
From 2012.igem.org
Contents |
Modeling
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.
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 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
- Each cell has a fixed probability of dying per time unit
- Each cell releases to the medium the AA it produces at a fixed rate
- The system will reach a steady state.
We’ve produce strains which are auxotrophic for one AminoAcid (AA). At the same time each strain produces the AA that the other strains need for their growth; therefore linking the survival of each module in the system. We will release these AA into the medium in the form of proteins rich in the A.A. The production and export can be regulated. A troyan protein ensures that the AA reaches the auxotrophic strain.
The deterministic model is based on the following assumptions:
- The growth rate will depend on the concentration in the medium of the AA that they lack; as its concentration in the medium [AA] increases so does the growth rate until it settles at the known value for yeast (doubling time 90 min).
- The system will reach a steady state.
This minimal model has four variables:
- two populations [N his-] and [N trp-]
|
|
- the concentration in the medium of the two aminoacids that they produce Histidine [his] and Tryptophan [trp].
| 
|
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 τ is the required to build the daughter. One of the hypotheses built into our model is that τ 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.
Steady State Solution
The 4 equations in [] were equaled to zero and the non-linear algebraic system solved using MATEMATICA to find their equilibrium values.
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.
{N_{a} ; N_{b}} → { N_{t} , x_{a} }
Three solutions were found:
# {{N_t→0},
#{Nt→Cc+(Cc√da √db death)/(√pa √pb), xa→(√db √pb)/(-√da √pa+√db √pb),
AAa→〖(-(Kb√pa √pb)/(√da √db kmax+√pa √pb))〗^(1/m),AAb→〖(-(Kb√pa √pb)/(√da √db kmax+√pa √pb))〗^(1/n)}
#{Nt→Cc-(Cc√da √db death)/(√pa √pb), xa→(√db √pb)/(√da √pa+√db √pb),
AAa→〖((Kb√pa √pb)/(√da √db kmax-√pa √pb))〗^(1/m),AAb→〖((Kb√pa √pb)/(√da √db kmax-√pa √pb))〗^(1/n)}
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!
- The first solution is clear. The culture will not thrive under certain conditions… to be determined!
- The second one has no physical relevance since it yields negative concentrations for the Amino Acids.
- However the last one yields :
N_t→Cc-(Cc √(da ) √db death)/(√pa √pb),
x_a→1/(√((da pa)/(db pb))+1)
Regulation
Note that the percentage of each population in the community is a function of √((da pa)/(db pb))
In particular
x_a→1/(√((da pa)/(db pb))+1)
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 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 γ, where <math> \frac{p_{a}}{p_{b}}= 10^{\gamma} </math>
Figure 1.
.. 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 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.
Figure 2.
Parameters
The parameter estimation process is shown in xx section(linked)
| 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 Analysis
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.
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 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.
First we’ve checked the validity of:
Nt=Cc ( 1 - ( √da √db death)/(√pa √pb))
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.
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). 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 there is no difference between the cases with high and low Ks. Note that it reaches the ss faster as well.
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. 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.
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!
Parameter Space and Solutions
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.
Example
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; 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. 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?
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 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
However, falling that, it was done numerically.
Relaxation time
Not only are we interested in predicting the fraction of each strain in the culture, we want it within a time frame. 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.
Figure 6.
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!!
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.
Figure 8.
Parameter selection
We've estimatated values for all the parameter in the model to check is viability.
- 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.
Figure aux his and aux trp.
The values were taken from a fit to f(x) = ODmax*x^n/(x^n + K_dil) with MATLAB’s TOOLBOX: Curve Fitting Tool
His- :
K_dil = 0.3279
ODmax = 62.53
n = 1.243
R-square: 0.9376
Trp- :
K_dil = 0.1770
ODmax = 57.46
n = 1.636
R-square: 0.9905
The results were then converted to the units chosen for the simulations.
K = K_dil * [AA 1x] # Avog / (Molar mass AA), where [AA 1x] = 0.02 mg/ml
- K trp= 0.1770 * 5.88e16 AA / ml.
- K his= 0.3279 * 7.8e16 AA / ml.
- The capacity of the system taken in the model is the one often used in the 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 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]
These parameters will be regulated experimentally (ε , p) to alter the fraction between populations.
- p trp = p*6e4*3600 aa / cell hour
- p his = p* ε*6e4*3600 aa / cell hour
where p < 1 controls the fraction of the top 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 [ref]) * relative abundance of AA * # Avog / AA's molar mass.
- d trp= 6.348e8 AA / cell
- d his= 2.630e7 AA / cell
