Team:British Columbia/Consortia

where K_S is the limiting substrate concentration when the specific growth rate is at half maximum (K_S = S when µ_g = µ_m/2). K_S is also known as the saturation constant or half-velocity constant[1].

Equations and Constants:

In order to predict the population dynamics, with the assumption of Monod kinetics, our great interests are: the concentration of Tryptophan, Methionine, and Tyrosine in the media available for each specific auxotrophs, as well as how these concentrations will change with respect of time.

As there is very limited initial availability of each limiting amino acid in the media, the observation of the coculture of all three auxotrophs (Maybe a link here to the coculture data??? ) proved that there is a production and accumulation of all three limiting amino acids. Based on this, the hypothesis is that each auxotroph is feeding on its limiting amino acid excreted by the other two anxotrophs to grow, and assumed that this growth-related consumption is much more significant that their own cell maintenance. For example, ΔmetA is not capable in making Methionine, which is essential for cell growth, so it grew in the coculture must because of that there is Methionine available in the coculture media (aka the environment), which is expected to be produced by the other two auxotrophs (ΔtrpA, and ΔtyrA), and all the Methionine in the environment is assumed to be used only for ΔmetA growth.

Based on this observation and reasoning described above, combing with the Monod kinetics, we proposed the following fundamental equations for our consortia model, where each of the variable is annotated as what they actually are (ΔtrpA, ΔmetA, and ΔtyrA is annotated as trpA-, metA-, and tryA- respectively; [Trp], [Met], and [Tyr] represents the concentration of Tryptophan, Methionine, and Tyrosine respectively; ODtrpA-, ODmetA-, and ODtyrA- stands for Optical Density of ΔtrpA, ΔmetA, and ΔtyrA respectively; a1, a2, and a3 are the consumption constants with respect to the changes of ODtrpA-, ODmetA-, and ODtyrA- respectively; µtrpA-, µmetA-, µtyrA- are the specific growth rates of ΔtrpA, ΔmetA, and ΔtyrA respectively; µMaxtrpA-, µMaxmetA-, µMaxtyrA- are the maximum growth rates of ΔtrpA, ΔmetA, and ΔtyrA respectively; and KstrpA-, KsmetA-, KStyrA- are the half-velocity constants of ΔtrpA, ΔmetA, and ΔtyrA respectively). Each of the equation term, and key constants within will be discussed in details later in this section.

In the above equations, some of the key constants for the Monod Kinetics are obtained through analyzing the growth curves of a 96-well plate environment, with various designed amino acid concentrations. Various characteristic growth rates were observed and plotted:

Figure 1: Monoculture Growth Rates at various Limiting Amino Acid Concentrations

Figure 2: Log Scale of Monoculture Growth Rates at various Limiting Amino Acid Concentrations

note: we had lower amino acid concentrations designed as well, though we found out that when concentration is below 1E-7 M, cell OD over night does not show significant differences through out the incubation period. However, their OD do still varies, in the magnitude of 0.01, and is hard to differ from controls and is considered no growth and shown as the one of the error bars. For detailed protocol, please go to Link:.............protocol.

Based on Figure 1 and 2, key constants such as maximum growth rates, and half-velocity constants for each auxotrophs are obtained [1]:

Further more, as discussed in the beginning of this section, we assumed that the amino acid consumption is solely depended on the cell growth and not the cell self maintenance. In order to find out the remaining constants, consumption rate constants (a1, a2, and a3), we need to find out a way to measure the amino acid concentrations in the media with respect to the change of its OD for each auxotroph.

To do this, firstly, we need to find a way to measure media amino acids concentrations. With no good access to HPLC, we noticed in previous experiments that each of the knockout monocultures will reach to a certain final OD at a different known amino acid concentrations, and also according to a recently published paper: A Programmable Escherichia coli Consortium via Tunable Symbiosis, Kerner A et. al (2012) confirmed that there indeed is a relationship between the final auxotroph ODs with the environmental amino acid concentrations, though may not be very accurate[2].

Thus we generated a calibration curve to relate the limiting amino acid concentrations to their final OD readings based on our own wet lab data for our specific ΔtrpA, ΔmetA, and ΔtyrA E. coli strains, which are shown in the figures below:

Figure 3: TrpA Auxotroph Final OD and Trp Concentration Calibration Curve

Figure 4: MetA Auxotrph Final OD and Met Concentration Calibration Curve

Figure 5: TyrA Auxotroph Final OD and Tyr Concentration Calibration Curve

Note: all the data are obtained through triplicates, and detailed protocol can be found in Link:.............protocol. Also, we tested the accuracy of one of the calibration curves. The blue dot in Figure 5 is another data value we obtained later at a known Tyrosine concentration (5E-5 M), and it is reasonably close to the predicted curve. This test point has proven that this calibration curve can serve us well at this stage of the research; however, we are aware that this curve has limitation, mainly due to the big error bars at the higher OD data points. For future accuracy improvement, more data points (more of the same concentrations as well as more of different concentrations) should be collected.

Secondly, it is necessary to find out how the amino acid in the environment is consumed with respect to cell growth.

We grew all three auxotrphs: ΔtrpA, ΔmetA, and ΔtyrA in 50 ml flasks with known initial Tryptophan, Methionine, and Tyrosine concentrations (5E-5 M) respectively, and harvested 5ml of cell free supernatants at various ODs. All these supernatants served as new 96-well plate media for each knock out cultures with newly added M9-Glucose shots (2% glucose in the well), which made sure the amino acids of interest remain as the limiting substrate for all knock outs. Typical protocol can be found Link:.............protocol.

After growing the inoculated 96-well plate in humidified temperature control (37 deg) plate reader for 15 hours, their growth curves had been captured and most of the monocultures had reached their stationary phase, and the plate was then transferred to 37 deg room to obtain the remaining final OD readings.

Based on the calibration curves (Figure 3, 4, 5), the amino acid concentrations in the wells were calculated from their final OD readings, and when corrected with a dilution factor of 5/4 (when we added in the M9-Glucose shot, it was 1/5 of the total well media volume), we arrived at the corresponding supernatant amino acid concentrations at various auxotroph monoculture OD readings. Again, with the assumption that the only change in the amino acid concentration in the media is due to cell growths, and as they are specific knock outs, they were not able to produce any amino acid into the environment, each auxotrophs' consumption constants with respect to cell growth rate (a1, 2, 3) can be calculated. Their values are:

Thirdly, after founding out the consumption constants for corresponding auxotrophs, we were interested in finding out how the other two auxotrophs will contribute to this specific amino acid changes in the environment. After all, these equations (f(ODtrpA-)s, f(ODmetA-)s, and f(ODtyrA-)s in equation (1), (2), and (3)) are the key to the accuracy of this model.

In the coculture environment, for the behaviour of auxotrophs with respect to changing the environment amino acid concentrations, the closest assumption we can make is that teach auxotrph will behave similarly to its behaviour in its monoculture condition. Thus, we used the same calibration curve for relating the final OD to the limiting substrate amino acid concentration in the media. The difference is that in stead of finding out the corresponding limiting amino acid concentration, e.g., ΔmetA media for [Met] measurements, this time, we are finding out the other two non-limiting amino acid concentrations of each knocked out E. coli auxotrophs, e.g. ΔmetA media for [Trp] and {Tyr] concentrations, at various OD readings [2]. With the consideration of sample taken time and some data manipulation, we were able to derive six functions that relate the change of all amino acid concentrations with respect to time (e.g. d[Trp]/dt), with the other two auxotrophs' optical densities (e.g. ODmetA- and ODtyrA-), which is shown below:

Figure 6: Change in Trp concentration as a function of ΔmetA OD

Figure 7: Change in Trp concentration as a function of ΔtyrA OD

Figure 8: Change in Met concentration as a function of ΔtrpA OD

Figure 9: Change in Met concentration as a function of ΔtyrA OD

Figure 10: Change in Tyr concentration as a function of ΔtrpA OD

Figure 11: Change in Tyr concentration as a function of ΔmetA OD

Based on these approximated functions, we were able to obtain the equations used in equations (1), (2), and (3):

One big limitation of the accuracy of these equations is that their accuracy is extremely dependent on the accuracy of the calibration curves (Figure 3, 4, and 5) as well as precision and accuracy of the OD readers, especially especially when OD readings are low. As you can see from Figure 3, 4, and 5, the slope is extremely high for all of them at low OD values, and thus are very sensitive to inaccurate low OD readings (OD<0.1).

However, some of the equations are derived from the available data assuming that at low OD values, it is still reasonable to use the calibration curve to calculate the supernatant (cell free media) amino acid concentration.

Especially for Equation (15), which is obtained from Figure 11, its raw experimental OD readings are shown below:

Figure 12: ΔtyrA auxotroph growth curve in various ΔmetA ODs cell free supernatants in 96-well plate reader environment

Note: As discussed before, all of this data can cause inaccurate supernatant concentration calculation due to its extremely low OD readings. All data set lines are labelled correspondingly, triplicates are noted by well 1, 2, and 3. As we can see, there are some distinct differences in the OD readings among the three different supernatant data sets; however, when corrected through focusing on the difference between the initial reading and the final reading with a base line of 0.025 (which is the lower value used for all the calibration curves), they average roughly the same, with an OD average value of 0.03847, 0.04043, and 0.04218, respectively, for ODmet 0.219, 0.546, and 0.764.

Matlab Code

Now with all the key functions and constants approximated, we are able to encode our model into Matlab, where the six key ODE functions can be solved and output the our prediction of population dynamics

The code is shown below, which contains a detailed step by step instruction for future iGEM teams or related scientific researches use, as well as the clarifications of our own code. It is also written is a way that it can be directly copied and pasted into Matlab blank m.files with all the comments in Matlab commenting format already, so that the coder can read the coded instructions and explanations in the coding command window without having to switching between, for example, the matlab and browser windows.

Code begins here:

% input the following commands, which are all in bold, into MATLAB command window, and not the m.file, before you call the ODE function.

%firstly, define a time array, set the upper limit to be the time (hr) you want this model to predict. For this model to be accurate, it is advised not to use extreme values.

%in our cases, the upper limit is set to 12 as some of the modelling equations are based on 11 hour culture supernatant data. Thus accuracy is expected within 12 hours.

% >> xspan = [0 12];

%in order to set up the initial conditions, where they corespond to the initial values of y(1), y(2), y(3)..., to y(6) in our case, and y(n) if you have n number of ODE equations

% >> ynot = [1E-7 1E-7 1E-7 0.01 0.01 0.01]';

% the first three values are the initial values for y(1), y(2), and y(3) respectively, which in out cases are the initial amino acid concentrations, and as we are adding the cells directly into the M9 media without supplementing any amino acids, it is believed to be at very low concentrations, and thus is set to 1E-7 mole/L.

% the remaining three values are the initial conditions for y(4), y(5), and y(6), which are the initial cell ODs in out case, and they are all set to 0.01, which is based on the lower values of the OD plate readings after correcting with blanks.

%Then RUN THE FUNCTION

% >>[X,Y] = ode45(@UBCiGEM2012_ConsortiaModel,xspan,ynot);

% this takes the function name in your m.file, the range of your interest, and the initial conditions which is just defined, respectively

%Note: have to make sure that the name is exactly the same as the function name in your m.file for it to be able to run correctly.

%the following code is for the actual mfile.

% 1) Declare the name of the function:

function dy = UBCiGEM2012_ConsortiaModel(t,y);

% 2) define all your ODE variables clearly for your own understanding as they can get quite confusing if you don't.

% in our case:

% y(1) = [Trp];

% dy(1) = d[Trp]/dt

% y(2) = [Met];

% dy(2) = d[Met]/dt;

% y(3) = [Tyr];

% dy(3) = d[Tyr]/dt;

% y(4) = [OD_trpA-];

% dy(4) = d[OD_trpA-]/dt;

% y(5) = [OD_tyrA-];

% dy(5) = d[OD_tyrA-]/dt;

% y(6) = [OD_metA-];

% dy(6) = d[OD_metA-]/dt;

% 3) define your constants first that will not depend on variables

% in our case, it was all measured through our own wet lab data and is explained in the previous section

%define maximum growth rate constants

u_maxtrpA = 0.2442;  % units = hr-1

u_maxmetA = 0.5524;  % units = hr-1

u_maxtyrA = 0.3332;  % units = hr-1

%define half-velocity constants

Ks_trpA = 1E-7;  % units = M

Ks_metA = 5E-6;  % units = M

Ks_tyrA = 8.36E-8;  % units = M

%define the consumption rate constants

a1 = 4.08E-4;  %units= M/OD

a2 = 1.44E-4;  %units= M/OD

a3 = 1.23E-5;  %units= M/OD

% 4) write out ODE CODE, suggest to write out all equations on paper first, as they are a lot more clear than coding them directly, and will reduce the chance for you to make error why coding significantly.

dy = zeros(6,1);

dy(1) = ((1e-8)*(exp(1))^(16.586*y(4))) + ((-0.023*(y(5))^2)+0.021*y(5)-3.5e-3) - a1*(u_maxtrpA*(y(1))/(Ks_trpA+y(1)))*y(4);

dy(2) = ((9e-6)*(y(6))^2-(7e-6)*y(6)+1e-6) + (((2e-6)*(y(4))^2)-(2e-6)*y(4)+9e-7) - a2*(u_maxmetA*(y(2))/(Ks_metA+y(2)))*y(5);

dy(3) = ((3e-7)*(y(4))^2-(8e-7)*y(4)+4e-7) + (((1e-6)*(y(5))^2)-(2e-6)*y(5)+6e-7) - a3*(u_maxtyrA*(y(3))/(Ks_tyrA+y(3)))*y(6);

dy(4) = (u_maxtrpA*y(1))/(Ks_trpA+y(1))*(y(4));

dy(5) = (u_maxmetA*y(2))/(Ks_metA+y(2))*(y(5));

dy(6) = (u_maxtyrA*y(3))/(Ks_tyrA+y(3))*(y(6));

% note: For easy coding and checking purposes, all this equations are exactly the same as the one shown in equation (1), (2), and (3) for dy(1), dy(2), and dy(3) respectively, with also the same arrangements of each equation term within them.

% and Please also note: these equation has been derived with direct simple substitution of equations, where equation (4), (5), (6), (7), (8), (9), (10), (11), (12), (13), (14) and (15) have been substituted into the equation (1), (2) and (3).

%if you want, you can choose to plot diagrams and display the predicting dynamics within the m.file or manually plot the arrays of your interests in command window

%For example, for this project, we are interesting in predicting the ODs of each auxotrophs, thus, we manually plotted values in y(4), y(5), and y(6) (which again, are the OD values of each auxotrophs) agaisnt time, by typing the following in the command window after run the ODE45 functions:

%>>plot(X,Y(:,4:6))

End of the m.file

Note: I will indicate here that the direct translation of the previous discussed equations with simple substitutions at either the beginning or the end of this section

Simulation and Analysis

When the previous matlab code was run, it predicted the population dynamics without inducing as shown in the following diagram:

Figure#. Population Dynamics Prediction within Confidence Range

Figure#. Population Dynamics Prediction for 16 hours

Figure#. Population Dynamics 16-hour Prediction after eliminating one potential inaccurate function

Figure#. Population Dynamics 12-hour Prediction after eliminating one potential inaccurate function

700x700px

As mentioned at the beginning, this code can be generalized to predict the population dynamics when there is a tuning effect, therefore we can use this model to see decide the amount of tuning we need to control the desired relative population of each auxotrophs.

References

[1] Shuler, ML; Kargi, F; Bioprocess Engineering Basic Concepts, 2002 2nd edition. Prentice Hall PTR.

[2] Kerner A; Park J; Williams A; Lin XN; A Programmable Escherichia coli Consortium via Tunable Symbiosis, 2012. PLoS ONE 7(3): e34032. doi:10.1371/journal.pone.0034032