on June 21, 2011 by in Metabolic modelling, Comments (2)

# Kinetic modelling of metabolic pathways

## Authors

Natalie J Stanford1,2 , Kieran Smallbone2,3
1. Integrative Systems Biology Doctoral Training Centre; 2. Manchester Centre for Integrative Systems Biology; and 3. School of Mathematics, University of Manchester, M13 9PL, UK.

## Abstract

We show how to build a kinetic model of a metabolic pathway. We provide the example of building a model of glycerol synthesis using Copasi. However, the techniques required remain the same for any pathway, using any software.

## Introduction

A mathematical description of a kinetic metabolic model may be given in differential equation form as $latex x^\prime = N v (x, y, p)$ $latex x(0) = x_0$ and this may be used as a guide as to the data required to create and parameterise a kinetic model. First, N is the stoichiometric matrix, which may be derived from the topology of the model. x denotes metabolite concentrations. y denotes boundary metabolites, whose concentrations are not allowed to vary, but do affect the reaction rates. Initial concentrations for both x and y must be deﬁned, though note that only concentrations x will change over time. Finally, v denotes reaction rates; these are dependent on kinetic mechanisms, concentrations x and y and parameters p.

## Stoichiometry

The ﬁrst stage in kinetic modelling is deﬁning the pathway of interest and its boundaries, which are captured mathematically in the stoichiometric matrix. As an example, we shall use a model of glycerol synthesis containing only two reactions [1], though the same procedure can be used for more complex systems. The model’s two reactions are glycerol 3-phosphate dehydrogenase (GPD) – which has reactants DHAP and NADH and products glycerol 3-phosphate (G3P) and NAD – and glycerol 3-phosphatase (GPP), which has reactant G3P and product glycerol (Gly) and phosphate (Phi). We shall show here how to implement the model using Copasi [2], but the techniques we describe will be applicable to any software. To enter the stoichiometric knowledge in Copasi, navigate to
Model > Biochemical > Reactions
Double click on the equation box, which should take you to the reaction interface (see ﬁgure 2). In the chemical equation box enter
DHAP + NADH = G3P + NAD
Click commit and navigate back to Model > Biochemical > Reactions

Figure 1: Navigating to the reaction interface

Figure 2: The reaction interface.

Note: You may notice that Copasi automatically deﬁnes a rate law. This is ﬁne, and will be addressed later in the tutorial. Double click the new reaction box, and enter the following in the chemical equation box G3P -> Gly + Phi Note: You may have noticed that in reaction 1 we use = to signify the conversion between substrates and products, whilst reaction 2 uses ->. The former is used to represent a reversible chemical reaction. The latter signifies an irreversible reaction. It is important to note that Copasi will not let you assign pre-deﬁned rate laws for reversible reactions that are irreversible and vice versa. However, if you choose to deﬁne your own rate law, these constraints can be violated. This will be discussed in more detail when we address entering rate laws. Navigating now to
Model > Biochemical > Species

Figure 3: The species interface

we see that the six metabolic species have already been added to the model. However, we do not want to simulate changes in concentration of the boundary metabolites in our system. To address this, change Simulation Type in DHAP from reactions to fixed, and repeat this for the other boundary metabolites, Gly, NAD, NADH and Phi.

Figure 4: Changing the metabolite boundaries

The stoichiometry matrix N is now found at
Model > Mathematical > Matrices
As expected, it contains only G3P, which is produced by the ﬁrst reaction and consumed by the second.

Figure 5: The stoichiometric matrix.

## Metabolites

In order to build a kinetic model, concentration values must be provided for all metabolites. Typically these will come from metabolomics measurements, or databases such as HMDB [3]. In Table 1, the concentrations for all metabolites in our model are set out. These may be added by accessing the Species menu in Copasi, and entering the concentration as the initial concentration (see ﬁgure 6).

Figure 6: Setting the initial concentration.

ADP, ATP and F16BP are included in the model as modiﬁers of GPD (Reaction 1). This is represented in Copasi by changing the reaction’s deﬁnition to
DHAP + NADH = G3P + NAD; ADP ATP F16BP

metabolite concentration (mM)
ATP 2.37
DHAP 0.59
F16BP 6.01
G3P 0
Gly 15.1
Phi 1

Table 1: Metabolite concentrations used in the model. It is important to ensure that consistent units are used throughout the model. Navigating to Model, we change Time, Quantity Unit and Volume Unit to min, mmol and l, respectively. In this way, all concentrations are measured in units of mmol/l = mM, and all ﬂuxes in units mM/min.

Figure 7: Model unit interface.

## Rate equations

Kinetic rate laws v may be derived from knowledge of the mechanism underlying the enzymatic process. In many cases these are not known, and we must instead resort to approximations such as linlog [4, 5], power law [6, 7] or convenience kinetics [8]. For our model, the rate of GPD is known to be
$latex v_1 = \frac{\frac{Vf1}{K1nadh \cdot K1dhap} \left(NADH \cdot DHAP – \frac{NAD \cdot G3P}{Keq1} \right)} { \left(1 + \frac{F16BP}{K1f16bp} + \frac{ATP}{K1atp} + \frac{ADP}{K1adp} \right) \left( 1 + \frac{NADH}{K1nadh} + \frac{NAD}{K1nad} \right) \left(1+ \frac{DHAP}{K1dhap} + \frac{G3P}{K1g3p} \right)}$
whilst the rate of GPP is given by
$latex v_2 = \frac{\frac{V2 \cdot G3P}{K2g3p}} {\left( 1+ \frac{G3P}{K2g3p} \right) \left( 1+ \frac{Phi}{K2phi} \right)}$
Copasi has a wide variety of inbuilt kinetic laws that may be accessed via the Functions menu. However, the forms above must be entered manually; this is achieved by returning to the Reactions menu and selecting GPD (Reaction 1). Click New Rate Law and enter the following
Vf1/(K1nadh*K1dhap)*(NADH*DHAP-NAD*G3P/Keq1)/((1+F16BP/K1f16bp+ATP/K1atp+ADP/K1adp)*(1+NADH/K1nadh+NAD/K1nad)*(1+DHAP/K1dhap+G3P/K1g3p))
Once the equation has been entered you then need to associate each identifier as a Substrate (NADH, DHAP), Product (NAD, G3P), Modifier (F16BP, ATP, ADP) or Parameter(all other identifiers). This can be seen in ﬁgure 8

Figure 8: Entering a new rate law.

Repeat for GPP (Reaction 2)
V2*G3P/K2g3p/((1+G3P/K2g3p)*(1+Phi/K2phi))
once again identifying Substrate (G3P), Product (Gly, Phi) or Parameter (all other identiﬁers). Returning to the Reactionsmenu, these formulae need to be selected again under Rate Law for each reaction. Where there are multiple substrates, products or modiﬁers within a reaction it is important to check that the correct identiﬁer is associated with the correct metabolite. These are circled in ﬁgure 9

Figure 9: Selecting the rate law.

## Kinetic parameters

The kinetic formulae above lead to a range of kinetic parameters. Ideally, these would be measured using enzymatic assay, When this is not possible, databases such as Brenda [9] or Sabio-RK [10] could be used. In Table 2, the parameter values in our model are set out. These can be modiﬁed in the same interface highlighted in ﬁgure 9.

parameter value units
V1 147 mM/min
Keq 10000 1
K1atp 0.73 mM
K1dhap 0.54 mM
K1f16bp 4.8 mM
K1g3p 1.2 mM
V2 53 mM/min
K2g3p 3.5 mM
K2phi 1 mM

Table 2: Parameter values used in the model.

## Analysis

Having built a model of the system, we would like to use mathematical tools to further our understanding of its behaviour. Copasi provides a number of tools to carry out these tasks under the Tasks menu. In the original paper [1], the authors investigate the response of ﬂux through the pathway to changes in boundary metabolite concentrations. This analysis may be performed running the Sensitivities task, with Function: Concentration Fluxes of Reactions and Variable: Initial Concentrations. Its results are presented in Table 3.

metabolite response
ATP -0.420
DHAP 0.475
F16BP -0.162
G3P 0
Gly 0
Phi -0.0741

Table 3: Scaled response coeffcients of ﬂux to metabolite concentrations.

## Systems biology standards

Describing mathematical models as above is unwieldy and error-prone, and naturally leads to difficulties in reproduction of results. Thus researchers have developed SBML (the Systems Biology Markup Language [11]). A computer-readable format for representing models of biological processes. Exchanging our models in SBML format (rather than e.g. Copasi’s own xml format) not only allows us to describe our model in unambiguous terms. It also gives us access to over 200 software tools specialising in many elements from data integration to visualisation [12]. SBML can be combined with MIRIAM [13] to annotate the entities of those models; for example, by marking-up the molecule G3P as CHEBI:15978 allows its unambiguous identiﬁcation and automatically links to many additional sources of information. This may be achieved in Copasi under the Species submenu Annotation. BioModels.net [14] is a modelling repository containing hundreds of models marked up with the above standards. If you are having trouble reproducing the glycerol model described above, it is available from the repository at BIOMD0000000076.

## References

1. G. Cronwright and J. Rohwer, “Metabolic control analysis of glycerol synthesis in Saccharomyces cerevisiae,” Applied and Environmental Microbiology, vol. 68, pp. 4448–4456, 2002.
2. S. Hoops, S. Sahle, R. Gauges, C. Lee, J. Pahle, N. Simus, M. Singhal, L. Xu, P. Mendes, and U. Kummer, “COPASI – a COmplex PAthway SImulator,” Bioinformatics, vol. 22, pp. 3067–74, Dec. 2006.
3. D. S. Wishart, D. Tzur, C. Knox, R. Eisner, A. C. Guo, N. Young, D. Cheng, K. Jewell, D. Arndt, S. Sawhney, C. Fung, L. Nikolai, M. Lewis, M.A. Coutouly, I. Forsythe, P. Tang, S. Shrivastava, K. Jeroncic, P. Stothard, G. Amegbey, D. Block, D. D. Hau, J. Wagner, J. Miniaci, M. Clements, M. Gebremedhin, N. Guo, Y. Zhang, G. E. Duggan, G. D. Macinnis, A. M. Weljie, R. Dowlatabadi, F. Bamforth, D. Clive, R. Greiner, L. Li, T. Marrie, B. D. Sykes, H. J. Vogel, and L. Querengesser, “HMDB: the Human Metabolome Database,” Nucleic Acids Res, vol. 35, pp. D521–6, Jan. 2007.
4. D. Visser and J. J. Heijnen, “Dynamic simulation and metabolic re-design of a branched pathway using linlog kinetics,” Metabolic engineering, vol. 5, pp. 164–76, July 2003.
5. K. Smallbone, E. Simeonidis, D. S. Broomhead, and D. B. Kell, “Something from nothing: bridging the gap between constraint-based and kinetic modelling.,” The FEBS journal, vol. 274, pp. 5576–85, Nov. 2007.
6. M. A. Savageau, Biochemical Systems Analysis: A Study of Function and Design in Molecular Biology. Addison-Wesley Pub (Sd), 1976.
7. J. J. Heijnen, “Approximative kinetic formats used in metabolic network modeling,” Biotechnology and bioengineering, vol. 91, pp. 534–45, Sept. 2005.
8. W. Liebermeister and E. Klipp, “Bringing metabolic networks to life: convenience rate law and thermodynamic constraints,” Theor Biol Med Model, vol. 3, p. 41, 2006.
9. M. Scheer, A. Grote, A. Chang, I. Schomburg, C. Munaretto, M. Rother, C. Sohngen, M. Stelzer, J. Thiele, and D. Schomburg, “BRENDA, the enzyme information system in 2011,” Nucleic Acids Research, vol. 39, pp. D670–D676, Nov. 2010.
10. U. Leser, F. Naumann, B. Eckman, U. Wittig, M. Golebiewski, R. Kania, O. Krebs, S. Mir, A. Weidemann, S. Anstein, J. Saric, and I. Rojas, Data Integration in the Life Sciences, vol. 4075 of Lecture Notes in Computer Science. Berlin, Heidelberg: Springer Berlin Heidelberg, 2006.
11. M. Hucka, A. Finney, H. M. Sauro, H. Bolouri, J. C. Doyle, H. Kitano, A. P. Arkin, B. J. Bornstein, D. Bray, A. Cornish-Bowden, A. A. Cuellar, S. Dronov, E. D. Gilles, M. Ginkel, V. Gor, I. I. Goryanin, W. J. Hedley, T. C. Hodgman, J.-H. Hofmeyr, P. J. Hunter, N. S. Juty, J. L. Kasberger, A. Kremling, U. Kummer, N. Le Nov`re, L. M.e Loew, D. Lucio, P. Mendes, E. Minch, E. D. Mjolsness, Y. Nakayama, M. R. Nelson, P. F. Nielsen, T. Sakurada, J. C. Schaﬀ, B. E. Shapiro, T. S. Shimizu, H. D. Spence, J. Stelling, K. Takahashi, M. Tomita, J. Wagner, and J. Wang, “The systems biology markup language (SBML): a medium for representation and exchange of biochemical network models,” Bioinformatics, vol. 19, pp. 524–531, Mar. 2003.
12. SBML software guide
13. N. Le Novère, A. Finney, M. Hucka, U. S. Bhalla, F. Campagne, J. Collado-Vides, E. J. Crampin, M. Halstead, E. Klipp, P. Mendes, P. Nielsen, H. Sauro, B. Shapiro, J. L. Snoep, H. D. Spence, and B. L. Wanner, “Minimum information requested in the annotation of biochemical models (MIRIAM),” Nature Biotechnology, vol. 23, pp. 1509–15, Dec. 2005.
14. BioModels database – A database of annotated published models”