Playing with silicon cells; a short course

Hans V. Westerhoff (in collaboration with Jacky L. Snoep, Frank Bruggeman & Barbara M. Bakker)

Introduction

With the arrival of completely sequenced genomes, one might have thought that biology is now understood. The reality is that it is not however. Visit www.bio.vu.nl/hwconf/Lectures/ for explications of why the sequence information is a great step forward but not yet enough. What is important for cell function, is the dynamic interaction of the macromolecules, either through communicating metabolites or by more direct contacts. Although in principle partly encoded in the primary sequence of the macromolecules, we are and will always be far from being able to calculate those dynamic interactions from the genome sequence (just like we cannot ab initio from a Schroedinger equation calculate the properties of a glucose molecule).

One may determine experimentally the relevant kinetic properties of the macromolecules and then start calculating from there. The challenge then is whether we understand the functioning of living cells on the basis of the experimental data.

The ultimate aim is to put all the relevant kinetic data for a living cell together and thus make a computer replicon of the living cell. That is what we call the Silicon Cell (cf. www.siliconcell.net ) and there is an Amsterdam-Stellenbosch based initiative in this direction. It may be clear that we are not even close to the aim of achieving this for a complete living cell. Yet, the project is under way, in that silicon replica of parts of cells’ are available. These are being made accessible on the world wide web by the Stellenbosch/Amsterdam group of Jacky Snoep: www.jjj.bio.vu.nl

If you open the following website: http://www.jjj.bio.vu.nl/, you will find the database of mathematical models of metabolism in various cell types. If you click on ‘Model database’ you will find a list of the models in the database. They are organized in three groups:

1. ‘Silicon Cell’ type models, which are based on experimentally determined kinetic input and have been or will be validated experimentally. These models are the starting point of a true Silicon Cell, a computer replica of the living cell.

2. Core models, which do not contain much detail and serve mainly to demonstrate and explore a possible mechanism or principle.

3. Demonstration models, developed for teaching of Metabolic Control Analysis.

Each of the models is accompanied by some background literature under ‘more’. If you click on ‘model’, you will enter the modelling laboratory. Try for instance ‘Glycolysis in Trypanosoma brucei’. The screen is divided in two parts: at the right hand side the reaction scheme and at the left hand side the parameters of the model and the initial values of the variables as they were used in the published version of the model. If one moves the pointer over an enzyme in the reaction scheme, its kinetic equation will be displayed. At the left hand side you are able to change the parameter values. Also you can indicate if you want to run a time simulation or calculate a steady state. If you select ‘evaluate’, the program will calculate what you asked for. The other models in the database work in a similar way.

Exercises:

Even when one is not directly interested in a silicon cell, the silicon parts of cells enable one to do biochemistry in silico. In this short course, we shall do this with a number of aims in mind:

1. Demonstrate how one makes a kinetic model (see also the accompanying reprint [Bakker et al. 1999, JBC 274, 14551]).

2. Demonstrate the meaning of flux control coefficients

3. Demonstrate some principles of metabolic control

4. Show that metabolism can explode

5. Look at the race against time in the early development of Xenopus laevis

Kinetic modelling:

The T. brucei model:

1. Go to the silicon web site (www.siliconcell.net ), then to the silicon cell site, then to model database, then to the Glycolysis in Trypanosoma brucei model. Inspect the pathway diagram to understand that: (i) in this organism glycolysis is organized in part in an organelle (a ‘glycosome’), (ii) in the aerobic situation there is no glycerol export from the organelle and 2 ATP's are produced per glucose in the cytosol, (iii) in the anaerobic situation one glucose plus one pyruvate could be produced hence only 1 ATP per glucose. Run your mouse pointer over a red oval corresponding to an enzyme and observe the rate equation for that enzyme appear in the bottom window and the name of the reaction appear on the plot. Check that these are realistic rate equations, not phenomenological. On the left hand side now note the parameter values, again realistic.

The reaction numbers refer to the following processes:

1. vGLCup: glucose transport

2. vHK: glucose phosphorylation

3. vPGI: glucose phosphate to fructose phosphate

4. vPFK: phosphofructokinase

5. vALD: aldolase

6. vTPI: triosephosphate isomerase

7. vGAPdh: glyceraldehyde-3-phosphate dehydrogenase

8. vGdh: glycerol dehydrogenase

9. vGPO: glycerol phosphate dehydrogenase (‘oxidase’; mitochondrial)

10. vPYRtr: pyruvate transport

11. vPGK: phosphoglycerate kinase

12. vPK: pyruvate kinase

13. vATPase: ATP hydrolysis though all other cellular processes

14. vGlyK: glycerol kinase

2. Now click on Sim and on Evaluate Model. Observe the result of the simulation. Now click on State and then on Evaluate Model. Whenever the box in the left hand lower corner begins to alternate between a yellow and a green border, you know that the server (=computer at the website; this one is doing the calculations, not your own computer) is working for you. Wait until the message evaluation complete appears in this box and then inspect the window that came up to see the steady-state values for concentrations and fluxes. Observe the calculated steady state concentrations and fluxes. Which flux value is most relevant for the ATP production, which is the sole purpose of glycolysis for T. brucei? Note down its value[Wff1] .

3. Under anaerobic conditions the trypanosome can produce ATP by producing equimolar concentrations of pyruvate and glycerol. Check the metabolic scheme to verify that this is possible. Then carry out the experiment by setting the V_max of the glycerol phosphate oxidase reaction (Vm9) to 0.368. Go the Vm9, double click on the 368.0 next to it, paint the 368.0 with your mouse, press delete on your key board, then type .368 (the computer does not like real zeros, so always use a small number instead of zeros) and Enter. Then press State and Evaluate Model again. What is the effect on the ATP [Wff2] production in the cytosol? This may still be sufficient ATP to survive. Also note down the NADH concentration[Wff3] .

4. Related to this, it has long been assumed that the enzyme triose-phosphate isomerase (TPI), interconverting DHAP and GAP, is a non-essential enzyme and consequently not suitable as a drug target. Explain why this is so; think what the flow pattern would become in the scheme. Would there still be respiration (vGPO, process 9) at steady state possible? Now check this by first resetting the parameter values (i.e. making the system aerobic again; Press ‘Reset’), recalculate to check that indeed ATP production is back up, then set the maximum rate of TPI (Vm6) to very low levels: Go the Vm6, double click on the 842 next to it, paint the 842 with your mouse, press delete on your key board, then type .842 and Enter. Then press State and Evaluate Model again. What happens to ATP production[Wff4] ? And to the NADH concentration[Wff5] ?

5. Indeed, when the enzyme was deleted experimentally, it proved to be essential for growth and survival of the trypanosomes. Can you explain from the results why TPI is essential? (Hint: what happens when you inhibit TPI (Vm6) and the mitochondrial glycerol-3-phosphate oxidase (GPO; Vm9) simultaneously[Wff6] ? And, look at the NADH concentration[Wff7] ).

Explain now in words [Wff8] that you have found a new and unexpected drug target.

The yeast model:

1. Click Model database, then the ‘Glycolysis in S. cerevisiae’ model. Click steady state, and Evaluate Model. Note that not all fluxes are equal. Do they balance at nodes? Is there a steady state? Now click Simulation and then Evaluate. How long does it take before a steady state is attained?

2. Increase the V_max of the glucose transporter (VmGLT) to 197.264 (do not forget to press ‘Enter’ after having done this). Then Click steady state, and Evaluate. What happens? What is wrong[Wff9] ? Set metabolites to be looked at to F16P and P (i.e. high energy phosphate only; check/uncheck the boxes on the right of compound and rate names), and then evaluate for 40 minutes rather than the default 10 (change ‘End time form 10.0 to 40.0, click Simulation and then Evaluate). You will see a metabolic explosion. This is a consequence of the turbo structure of the glycolytic pathway. Explain what this might mean[Wff10] . This is also called ‘substrate accelerated death’.

3. Reduce the activity of hexokinase (VmGLK) until the metabolic explosion disappears (do not go below 125). What is this value[Wff11] ? How could the cells protect themselves against substrate accelerated death[Wff12] ?

4. Reset the parameters, set the Vmax of the transporter (VmGLT) again to 197, and now determine for what extracellular glucose concentration (GLCo) the explosion would disappear. Would you expect the cells to be able to grow at high glucose concentrations[Wff13] ? And at low glucose concentrations[Wff14] ?

Early development:

Go to the Kinetics of histone gene expression model. This describes the race against time in early embryos of Xenopus laevis. The early embryos have a very short cell cycle and a large genome. Eukaryotic DNA needs to be packaged in histone protein, at a ratio (mass/mass) of 1/1 and therefore lots of histone protein is required, but how can sufficient of that protein be made?
Plot DNA as a function of time, first without the DNA replication attenuation box ticked (Have the DNA content circle ticked then click on Evaluate Model). The line shows the calculated DNA as a function of time, calculated for a constant, short cell cycle time (some 20 minutes, an absolute record for a vertebrate!). This is comparable to the cell cycle time of E. coli and within 24 hours, X laevis would overgrow the world if this were to continue. Now click the box DNA replication attenuation, which puts into the model a gradual increase of the cell cycle time [inspect the rate equation for DNA synthesis to see this; this sets n from 0 to 1], and then leads to a fit of the experimental data (the black dots). At early times there is no correspondence to experimentally observed total DNA per embryo (the dots) and calculated nuclear DNA (the line). Would you have an explanation for this[Wff15] ?
Show that insufficient protein is made if there is inheritance of neither histone mRNA nor histone protein, stored in the unfertilized egg (this is the his RNA store value that one can adjust) nor Histone (protein) store (Put His RNA store to 5.7E1; this means 5.7 times 10to the power 1 and put Histone store to 1.33E1). Now select the plot Histone per DNA; remember that the logarithm of 1 equals 0). Is sufficient histone protein made [Wff16] (remember the histone protein to DNA ratio should exceed 1)?
Would either store alone suffice[Wff17] ?
Reset the parameter values. Then set the parameter ktranscr. Att. to 1. This means that there is no transcriptional regulation. Show that this leads to vast access of histone protein. When? Much too much energy would be wanted to make all this protein.
Double ktranscr.att to 2. Inspect the rate equation for transcription. At what time is transcription attenuated in the model[Wff18] ? What regulatory effect does this have? Is there still too much histone protein?
Vary ktranscr.att until you find a value that removes the excess histone protein problem. What value do you need for this parameter? Explain in words what you think that the system needs to function well. If this is right then you could check this by inspecting the histone mRNA content. Now plot histone mRNA (RNA content) for the value of ktranscr. att. you found to be optimal, and also for a value of 1 for this parameter. What difference does it make?
Note the black dots in this plot, which correspond to the experimental result. What do you conclude about the realism of the transcription attenuation that you propose?
Is the value for the trancription attenuation parameter you proposed O.K? If not, try to change it so that the mRNA data are fitted better. What value of ktranscr.att is best? Check whether for this value there is not too much excess histone protein and enough histone protein.
Suggest another regulation phenomenon that could have the same effect on the histone protein to DNA ratio. How could you test this possibility?
Reset the parameter values. Figure out why there are 50 histone genes instead of one per haploid genome.

MCA in T. brucei

Return to the silicon T. brucei.

1. Increase the V_max of the fourth reaction (vPFK) by 10 % (type the new value in the box, by eliminating the 780 and retyping 858 and pressing Return) and calculate the percentage increase in glycolytic flux[Wff19] . Estimate the flux control coefficient of this enzyme[Wff20] . The enzyme is phosphofructokinase; the suspected rate limiting step of the pathway. Is it rate-limiting[Wff21] ?

2. If you click on MCA, you can indicate if you want to calculate control coefficients, elasticity coefficients: Now click Reset, MCA and then again Evaluate and read the flux control coefficient by reaction 4 (vPFK; 0.00428 etc).

3. Try to find the step with the largest flux control coefficient[Wff22] and then reduce the Vmax of this step by 50 % (Press reset, half the Vm of this step, press Return, click State and Evaluate Model) and observe whether this will reduce the flux of ATP production[Wff23] . (look at vATPase [process/flux 13] for this). Would this be a good drug target[Wff24] ?

4. Click ‘Reset’ then reduce Vm9 (mitochondrial respiration) to 0.368 to simulate virtual anaerobiosis. Show that the ATP yield drops to one per glucose[Wff25] . Is the rate still sufficiently high[Wff26] ?

5. Click MCA again and Evaluate a steady state. Are the same enzymes limiting the flux as in the aerobic case[Wff27] , and to the same extent[Wff28] ? Add all control coefficients with respect to any flux. What is the result[Wff29] ? Same for the control coefficients with respect to any concentration[Wff30] ?

6. Click ‘Reset’ to return to the aerobic case. Calculate the flux control coefficients of glucose transport, hexokinase and glyceraldehyde-3-phosphate dehydrogenase over the steady-state glucose consumption flux.

7. Calculate the control coefficients of the same enzymes, but now over the steady-state pyruvate production flux. Are they the same as the flux control coefficients calculated under 6? Why, or why not?

8. For drug design we are interested in a large inhibition of the flux. The question is to which extent the flux control coefficients can predict what happens upon a strong inhibition of individual enzymes. Decrease the V_max of glucose transport (process 1) by 50 % and calculate the effect on the glucose consumption flux. Repeat this for hexokinase (process 2) and glyceraldehyde-3-phosphate dehydrogenase (process 7). Is the effect on the flux what you expected from the flux control coefficients calculated under b? Why, or why not?

Explore further as you like; lots remains to be discovered.