One shortcoming in the use of these equations is their dependency on the accuracy of the calibration curves (Figure 3, 4, and 5), as well as the OD readings, especially when OD readings are low and noisy (< 0.1, nearing the limits of the plate reader's accuracy).

For example, one of the experimental raw data obtained are shown below:

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

.

Note: As discussed before, this may lead to inaccuracies in calculation from low OD readings. All data sets are labelled correspondingly, and triplicate measurements are noted by wells 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 by focusing on the difference between the initial reading and the final reading with a base adjustment of 0.025 (the lower value used for all the calibration curves), they average roughly the same, with an average OD of 0.03847, 0.04043, and 0.04218, respectively, for ODmet 0.219, 0.546, and 0.764.

To illustrate this point more clearly, here is a comparison with another set of growth curve data run in the same plate.

Figure 6: Comparison between ΔtrpA and ΔtyrA auxotroph growth curve in the same supernatants of ΔmetA at its various optical density media.

Note that all three of the trpA^- sets of growth curve grew to a much higher OD (reaching about 0.17) at the end of 15 hours and did not reach stationary phase, whereas all tyrA^- sets grew significantly less, no significant change in OD at the end of the experiment.

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.

Of course, copying and pasting the code here is welcomed!



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 = [4E-5 4E-5 4E-5 0.05 0.05 0.05]';

% 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 a low concentrations, though there will be some left over amino acids within the living initial cells, and thus is set to 4E-5 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.05, which is the initial OD reading we obtained after we innoculated the culture, with the assumption that all the cells are still alive.

%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.23; % units = hr-1

u_maxmetA = 0.34; % units = hr-1

u_maxtyrA = 0.28; % units = hr-1

%define half-velocity constants

Ks_trpA = 3.30E-7; % units = M

Ks_metA = 5.10E-6; % units = M

Ks_tyrA = 8.50E-7; % units = M

%define the consumption rate constants

a1 = 3.05E-4; %units= M/OD trp

a2 = 1.15E-4; %units= M/OD met

a3 = 1.4E-4; %units= M/OD tyr

% 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) = 2.5E-7 + 7.2E-8 - a1*(u_maxtrpA*(y(1))/(Ks_trpA+y(1)))*y(4);

dy(2) = 1.7E-7 + 5.1E-8 - a2*(u_maxmetA*(y(2))/(Ks_metA+y(2)))*y(5);

dy(3) = -3.6E-8 + 1.7E-6 - 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

Simulation Results, Analysis and Future Work

When the Matlab Code was run, it generated a graphic result predicting the population dynamics without inducing, shown below:

Figure 7. A comparison of model predictions versus empirical measurements (n=3) of each auxotroph’s growth in a co-culture.

As shown in the figure above, this model can give reasonable predictions though still can be further modified for better accuracy. For our future work, we will refine our model with more experimental data. This can be accomplished by running the amino acid experiment data points of interest. Alternatively, we can use more direct methods of amino acid measurement, like HPLC analysis.

This code can also be adapted to predict the population dynamics of our cocultures when tuning effects are introduced. This is accomplished by by generating six new functions to describe the amino acid dynamics of each monoculture. We can then use this model to decide the amount of tuning (by adding in different amount of inducer) needed to obtain a desired population of each auxotroph.

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