U+003F

Let’s return to our ODE system one last time.

In previous posts we’ve described how to determine the stoichiometric matrix N and the kinetic formulae v. Remaining is parameterisation of the model, which is heavily dependent on locating relevant data.

x denotes metabolite concentrations, which for our pathway are ADP, ATP, glucose 6-phosphate, glucose, lactate, NAD, NADH, oxygen and pyruvate. y denotes boundary metabolites, whose concentrations are not allowed to vary but do affect the reaction rates; here these are extracellular glucose, lactate and oxygen. Initial concentrations for both x and y must be defined, but only concentrations x will change over time. For human models like this one, HMDB [1] is a great resource, but we’ll continue to take our numbers from the BrainCirc model [2] wherever we can.

p denotes the kinetic parameters, here and typically Vmax and Km. A number of databases for Michaelis-Menten constants exist, including Brenda [3] and Sabio-RK [4]. Nonetheless, there is a real paucity of this sort of data available; the “data deluge” hasn’t arrived. Rather than using 14 parameters to model enzymes like PFK [5], we look to simpler kinetic forms or parameter-free methods such as constraint-based analysis.

Vmax parameters prove even more difficult to find. They typically vary wildly, even within a single cell-type, and are highly dependent upon environmental conditions. However, fluxes of glucose (7×10^-5 mM/s) and oxygen uptake (4×10^-4 mM/s) are known [2]. Assuming a steady state, this is enough to infer fluxes throughout the system. Using these fluxes and other known parameter values and concentrations, Vmax values may be inferred.

I’m sure you readers are more than capable of putting together a model from the info above and their favourite modelling software. But save your sticky-back plastic, here‘s one I made earlier [6]. Have a play!

References

1. HMDB: the Human Metabolome Database. www.hmdb.ca
2. Banaji M, Tachtsidis I, Delpy D, & Baigent S (2005). A physiological model of cerebral blood flow control. Mathematical biosciences, 194 (2), 125-73 PMID: 15854674
3. Brenda: The Comprehensive Enzyme Information System. www.brenda-enzymes.org
4. Sabio-RK: The Sabio Reaction Kinetics Database. sabio.villa-bosch.de
5. Teusink B, Passarge J, Reijenga CA, Esgalhado E, van der Weijden CC, Schepper M, Walsh MC, Bakker BM, van Dam K, Westerhoff HV, & Snoep JL (2000). Can yeast glycolysis be understood in terms of in vitro kinetics of the constituent enzymes? Testing biochemistry. European journal of biochemistry / FEBS, 267 (17), 5313-29 PMID: 10951190
6. SBML model: normal_cell.xml

Let’s return to the differential equation description of our model.

v denotes reaction rates; these are dependent on kinetic mechanisms, concentrations (x and y), and parameters (p). A number of flavours of kinetic mechanisms can be used, varying in their level of detail. Let’s take the hexokinase reaction as an example:

glc + atp ---> g6p + adp

One of the least complex descriptions of this reaction is given by linlog kinetics [1]

whereby reaction rates are assumed to vary logarithmically with metabolite concentration. Here ε is the elasticity as defined by metabolic control analysis (MCA) and linlog kinetics can be thought of as a dynamical interpretation of MCA. Linlog is really just a linearisation of the kinetics about a reference state; on the one hand this means it has some nice analytical properties [2], but on the other the representation is only valid in a region near this state.

Michaelis-Menten kinetics (o) are fairly linear when plotted in log-space, hence linlog kinetics (-) provide a reasonable approximation

If we have fully experimentally-characterised hexokinase and know that G6P is a competitive inhibitor of glucose, and ADP a competitive inhibitor of ATP, we could use a rate law such as [3]

Sitting comfortably in between linlog and full kinetics are various simplified Michaelis-Menten-like equations, such as convenience kinetics [4], or the dynamics of Murad Banaji‘s BrainCirc model [5]. We’ll follow Murad’s lead, defining

Whilst this form requires less experimental data than a full description of the reaction, it is saturative in its reactants and activators, and thus will extrapolate better than linlog kinetics.

The remaining kinetic formulae are given by

2. glycolysis:
3. pyruvate to lactate:
4. krebs:
5. oxidative phosphorylation:
6. atp use:

Next up: parameterisation.

References

1. Visser D, & Heijnen JJ (2003). Dynamic simulation and metabolic re-design of a branched pathway using linlog kinetics. Metabolic engineering, 5 (3), 164-76 PMID: 12948750
2. Smallbone K, Simeonidis E, Broomhead DS, & Kell DB (2007). Something from nothing: bridging the gap between constraint-based and kinetic modelling. The FEBS journal, 274 (21), 5576-85 PMID: 17922843
3. Teusink B, Passarge J, Reijenga CA, Esgalhado E, van der Weijden CC, Schepper M, Walsh MC, Bakker BM, van Dam K, Westerhoff HV, & Snoep JL (2000). Can yeast glycolysis be understood in terms of in vitro kinetics of the constituent enzymes? Testing biochemistry. European journal of biochemistry / FEBS, 267 (17), 5313-29 PMID: 10951190
4. Liebermeister W, & Klipp E (2006). Bringing metabolic networks to life: convenience rate law and thermodynamic constraints. Theoretical biology & medical modelling, 3 PMID: 17173669
5. Banaji M, Tachtsidis I, Delpy D, & Baigent S (2005). A physiological model of cerebral blood flow control. Mathematical biosciences, 194 (2), 125-73 PMID: 15854674

Creating a kinetic model is not as hard as they’ve had you believe. Over the next few posts, we’ll walk through the steps needed to build and analyze a small model of cancer metabolism.

A mathematical description of a kinetic metabolic model may be given in differential equation form as

and this may be used as a guide as to the data required to create and parameterize such a model.

First, N is the stoichiometric matrix, which may be derived from the network structure. Consider, for example, the representation of metabolism in mammalian cells, taken from Gatenby & Gillies (2004) [1].

We can use this as a guide to the important processes we’d like to include in our model:

1. hexokinase
   glc + atp ---> g6p + adp
2. glycolysis
   g6p + 3 adp + 2 nad ---> 2 pyr + 3 atp + 2 nadh
3. pyruvate to lactate
   pyr + nadh ---> lac + nad
4. krebs
   pyr + adp + 5 nad ---> atp + 5 nadh
5. oxidative phosphorylation
   o2 + nA ADP + 2 nadh ---> nA ATP + 2 nad
6. atp use
   atp ---> adp
7. glucose diffusion
   glcX ---> glc
8. oxygen diffusion
   o2X ---> o2
9. lactate diffusion
   lac ---> lacX

The parameter nA is included because the total ATP produced during oxidative phosphorylation can vary from cell to cell. Whilst there is a theoretical yield of 38 ATP from full oxidation of glucose, various inefficiencies mean that in reality the maximum is more like 30 ATP [2]. Here, full oxidation of glucose is achieved by

glc + 6 o2 + (4+6 nA) adp ---> (4+6 nA) atp

so we may estimate nA = 4.3.

References

1. Gatenby RA, & Gillies RJ (2004). Why do cancers have high aerobic glycolysis? Nature Reviews Cancer, 4, 891-899 PMID: 15516961
2. Rich PR (2003). The molecular machinery of Keilin’s respiratory chain. Biochemical Society Transactions, 31, 1095-1105 PMID: 14641005