Team:SDU-Denmark/Tour32

Modelling

Revealing the bits behind the simulation

An introduction
We aim to optimize the conditions already present in the cell to optimize the overall production of rubber. We’ve been looking at two quite different ways of optimizing the rubber production, which should be applicable to a real-life production. These two are:

• Optimization of each cell
• Optimization of the number of cells

The optimization of each cell identifies rate-limiting steps in the pathways related to the production of the polyisoprene. Where the optimization of the number of cells is to make in vitro tests of how a cell responds to different changes in environment, and find the environment with the best influence.

Optimization of each cell
In the production of rubber, there are three primary pathways: the glycolysis, the MEP-pathway, and the rubber transferase itself. We decided to focus on the MEP-pathway.

Since no one had modelled the entire pathway before, we build the model from the ground up. We based the model on Michaelis Menten kinetics, and the kinetic terms entering the model are: ordinary uni-uni reaction, ordinary bi-bi reaction, reversible bi-bi reaction and substrate inhibitory bi-bi reaction.

Results

Figure 1 - We're interested in the largest slope (velocity) of the curve, here we see that the increase in slope stops after a certain increase in the starting velocity. A rate limiting step! Rate limiting steps
Articles have shown the reaction catalysed by dxs to be the rate limiting step, Source: Sprenger GA, Schorken U, Wiegert T, Grolle S, De Graaf AA, Taylor SV, Begley TP, Bringer-Meyer S, Sahm H (1997) Identification of a thiamin-dependent synthase in Escherichia coli required for the formation of the 1-deoxy-D-xylulose 5-phosphate precursor to isoprenoids, thiamin, and pyridoxol. Proc Natl Acad Sci USA 94:12857–12862 (Link) and we asked the question: If we were to manipulate this step, which step would then be rate limiting? Simulating reaction kinetics often leads to stiff systems, and this is no exception. Traditional methods of solving differential equations therefore did not suffice, and MATLAB’s ODE15s, which is an L-stable solver, was used.

Simulations showed that the next rate limiting step would occur after an improvement in rate of conversion of 12 percent, and that this new rate limiting step was the reaction catalysed by IspG, as seen below. This result is worth keeping in mind when deciding on the design of the plasmids. Namely, if there is time to spare, then cloning IspG into the system might further enhance functionality.

Repeating the simulation with the increased velocities (i.e. post-modifying of Dxs and IspG), it is seen that IspF catalysed the third rate limiting step. IspF would thus be an enzyme of interest, if the system was to be further improved.

Economy
Figure 2 - Here we see how much dxs reacts to B₁ vitamin. As it's reversible, it stops producing B₁ vitamin faster, than the cell can use it, which is indicated by the horizontal line. But still what we throw exponentially more and more out as a function of the speed of the pathway. The second question posed concerned the economy of the system. That is, does the increase in velocity through the pathway result in increased product, or does the substrates exit the pathway at an unintended point. Based on the literature, intermediates after the reaction catalysed by IspD pass directly through the pathway, while there exists a reversible exit to vitamin B₁ from DXP.

It has been implicitly established that the increase in velocity through the pathway does not exclusively represent an increased exit to vitamin B₁. Otherwise, there would be no increase in product, at all. So, the question is (in plain English) how much, exactly, can we increase the speed of the pathway before it starts leaking - before the amount diverted to the exit surpasses the increase in amount of product. This can be seen on the plot below.

It is noticeable that the increase in speed leads to an increase in the B₁ vitamin producing term. At increases up to 100 percent, the exit term is unproblematic, while it seems at increases above 200% that the reaction would become uneconomical, as growth of this term is exponential. However, this is not of concern, as such increases in speed through the pathway are far higher than what we hope to achieve.

In conclusion, the simulations showed the rate limiting steps to be the reactions catalysed by Dxs, IspG, and IspF. We chose to concentrate our efforts on overexpressing Dxs and IspG. The choice not to clone IspF was made in the light of the fact that an increase in speed through the pathway of more than 100 % was futile in terms of economy.

For a simpler approach on what we actually did in our model you can visit our off-iGEM-server interaktiv model, made to give a simplified glimpse of the mehanism of the model.

Sources for enzyme kinetic details of MEP pathway

• Dxs: Sprenger GA, Schorken U, Wiegert T, Grolle S, De Graaf AA, Taylor SV, Begley TP, Bringer-Meyer S, Sahm H (1997) Identification of a thiamin-dependent synthase in Escherichia coli required for the formation of the 1-deoxy-D-xylulose 5-phosphate precursor to isoprenoids, thiamin, and pyridoxol. Proc Natl Acad Sci USA 94:12857–12862 (Link)
• Dxs: Miller B, Heuser T, Zimmer W (2000) Functional involvement of a deoxy-xylulose 5-phosphate reductoisomerase gene harboring locus of Synechococcus leopoliensis in isoprenoid biosynthesis. FEBS Lett 481:221–226 (Link)
• Dxr: Koppisch AT, Fox DT, Blagg BS, Poulter CD. E. coli MEP synthase: steady-state kinetic analysis and substrate binding. Biochemistry. 2002 Jan 8;41(1):236-43. (Link)
• IspD: Richard SB, Lillo AM, Tetzlaff CN, Bowman ME, Noel JP, Cane DE. Kinetic analysis of Escherichia coli 2-C-methyl-D-erythritol-4-phosphate cytidyltransferase, wild type and mutants, reveals roles of active site amino acids. Biochemistry. 2004 Sep 28;43(38):12189-97. (Link)
• IspE: Bernal C, Mendez E, Terencio J, Boronat A, Imperial S. A spectrophotometric assay for the determination of 4-diphosphocytidyl-2-C-methyl-D-erythritol kinase activity. Anal Biochem. 2005 May 15;340(2):245-51. (Link)
• IspF: Geist JG, Lauw S, Illarionova V, Illarionov B, Fischer M, Gräwert T, Rohdich F, Eisenreich W, Kaiser J, Groll M, Scheurer C, Wittlin S, Alonso-Gómez JL, Schweizer WB, Bacher A, Diederich F. Thiazolopyrimidine inhibitors of 2-methylerythritol 2,4-cyclodiphosphate synthase (IspF) from Mycobacterium tuberculosis and Plasmodium falciparum. (Link)
• IspG: Zepeck F, Gräwert T, Kaiser J, Schramek N, Eisenreich W, Bacher A, Rohdich F. Biosynthesis of isoprenoids. purification and properties of IspG protein from Escherichia coli. J Org Chem. 2005 Nov 11;70(23):9168-74. (Link)
• B₁ and B₆: Yeh JI, Du S, Pohl E, Cane DE. Multistate binding in pyridoxine 5'-phosphate synthase: 1.96 A crystal structure in complex with 1-deoxy-D-xylulose phosphate. Biochemistry. 2002 Oct 1;41(39):11649-57. (Link)
• B₁ and B₆: D. E. Cane, S. C. Du, J. K. Robinson, Y. J. Hsiung, I. D. Spenser. Biosynthesis of Vitamin B₆: Enzymatic Conversion of 1-Deoxy-D-xylulose-5-phosphate to Pyridoxol Phosphate J. Am. Chem. Soc. 121, 7722–7723 (1999). (Link)

Modelling with the team of Edinburgh

Optimization of the number of cells
Together with the team of Edinburgh, we implemented and ran a model created by Tobias Bollenbach (2009), Source: Bollenbach, T. et al. Non-optimal microbial response to antibiotics underlies suppressive drug interactions; Cell. 2009 November 13; 139(4): P. 707–718. (Link) using measurements done by the team of Edinburgh on their system. We’re using a whole cell model, aiming to predict an entire cell's behavior and growth. The model has to restrict itself at some point, so in our model, we include the most essential features of the metabolism of the cell; resource availability, growth rate, and intracellular protein levels. Most important is growth rate, as this allows us to test different environments and see their influence on the rate. As mentioned earlier, we simulated the system created by Edinburgh using their measurements. The resulting simulation can be found on Edinburgh’s wiki page as can the results. Though we haven’t implemented our own measurements, it definitely stands as proof of concept, and when a final production stand, we’ve shown it possible to optimize the environment, so enhance the production.

Figure 3 - The four main differential equations in Tobias Bollenbach's model.

So far, it has only been possible to solve these equations by solving the dynamics, even though we’re only interested in the steady state. So it should be possible to equate the differential equations to zero (steady state), and “just” solve that system. But so far only solving the dynamics have showed to be the only option which worked.

We solved the dynamics by first using a preprogrammed function solver in python to solve the nonlinear equations. To solve the dynamics itself we first used the implicit euler, which is a Radau IIA method of order one. As a prove of concept, implicit methods were necessary because of stiffness. Using such easy method we were able to prove the concept, but it has the same structure as the SDIRK method, which is of order maximum three. So if necessary we can improve the error by easily change the solver method.