Dynamical Modelling of a Metabolic Reaction Network

(1)

Dynamical Modelling of a Metabolic Reaction Network

JOAN GONZALEZ HOSTA

Masters’ Degree Project Stockholm, Sweden Sep 2009

XR-EE-RT 2009:019

(2)

Dynamical modelling of a metabolic reaction network

Master Thesis in Automatic Control, KTH

Joan Gonzalez Hosta

Supervisor: Elling W. Jacobsen

External contact: Veronique Chotteau

Stockholm-Barcelona, 2009

(3)

1 I would like to specially thank my supervisor Elling W. Jacobsen for making me realise when I was wrong and making me have a better understanding of everything, the responsible person in the KTH Division of Bioprocess Veronique Chotteau for always being able to find time for my doubts in her busy timetables and for the motivation she gave me in hard times and the “guy of the experiments” Andreas Andersson for his extra-hours in order to try to have the data on time and for his kindness. I’m amazed of all the things that I have learnt.

Joan Gonzalez Barcelona, 2009-08-11

(4)

2

Abstract

In a cell culture system, a lot of compounds are involved considering the ones contained in the medium (extracellular) and the intracellular environment. Studies of the different system behaviours adopted due to variations in the medium concentrations are typically based on empirical or statistical approaches when working with such an amount of compounds. In this paper is presented a method to develop a dynamical model of a cell culture system which can predict the evolution of the medium compounds and cell growth along the time for this kind of systems. The method has been implemented as a standardised tool that can be used for any kind of biological system when a pre-defined metabolic network and data of different states are available. In parallel to this project, a medium where the exact concentrations of all the compounds are known and should let the cells live properly have been developed by the KTH Division of Bioprocess. The system corresponding to the experiments carried out to develop the medium has been modelled using the proposed method. The results have been analysed and the conclusions, requirements and drawbacks of the model have been discussed.

(5)

3

1. Introduction

The phenomenon that maintains the cells living by providing feeding and consequently also energy is the metabolism. The metabolism is the set of chemical reactions that allows these organisms to grow and reproduce, maintain their structure and adapt to the environments where they have to live. The biochemical reactions that form this metabolism consist, basically, in the transformation of chemical compounds (or metabolites) into other ones. In order to study the cells behavior or how they evolve in different conditions, the reactions that form the metabolism have been tried to model into equations. The fact is that due to the huge number of reactions and compounds that form the metabolism, this modeling has been focused in some sets of reactions that are part of the metabolism and form a system themselves. Depending on the field of study, the target compounds or the processes of interest, the selection of the reactions and the metabolites that will be involved in the model can change for the same system. The collection of these reactions is what we call a metabolic reaction network.

1.1. Basic concepts

A metabolic reaction network is a complex system formed by several metabolites which are linked by biochemical reactions. The change in the concentrations of these metabolites depends basically on the magnitudes of the fluxes of the reactions where they are involved and, at the same time, the magnitude of the fluxes depend on the concentrations of the metabolites that are involved in every reaction. A metabolic pathway is a consecution of reactions within the network that starting from one metabolite leads to another one. So, according to the theory, the metabolic pathways within the metabolic network that will be

Figure 1. 1 Scheme of a metabolic network; metabolites (nodes) and reactions (arrows) [1]

(8)

6 followed by the different compounds will change depending on the concentrations of the metabolites involved in these pathways. The magnitude of the fluxes of the different pathways will determine the metabolites production or consumption rates and the change in their concentrations.

For every reaction, the metabolites that are consumed are called substrates (their concentration decreases when the reaction takes place) and the ones that are produced are known as products (in this case, the concentration increases).

A simple example of a typical metabolic network of mammal cells is represented in figure 1.1.

The nodes represent the metabolites that are considered in the network and the arrows represent the reactions that link these metabolites. The direction of the arrows means the sign of the reaction flux. For example, the reaction v1 transforms the metabolite G into G6P;

through this reaction, the concentration of the metabolite G (the substrate) will decrease while the concentration of G6P (the product) will increase and the velocity is the flux value of v1.

In order to carry out studies in this field, cells are cultivated, usually in reactors or flasks. A cell culture is the system formed by a medium where the cells live and is the environment that they have to interact with, and the cells themselves. The medium will also contain metabolites that will be used in the biochemical reactions that form the metabolism.

Then, considering a cell culture, two different kinds of metabolites can be found; the internal metabolites, which are the ones inside the cell, and the external metabolites, that are the compounds that the medium contains. The metabolic network describes the reactions that relate the metabolites (internal and external) that are wanted to be considered in the system.

In the example in figure 1.1 the metabolites in a grey circle represent the external metabolites, which are also linked to internal metabolites through reactions. These processes are what allow the interaction of the cells with the medium.

1.2. The aim of the master thesis

The aim of the master thesis is to develop a method that produces a mathematical model that allow us to predict the behaviour of a determined cell culture system given the reactions of a pre-defined metabolic network and implement it. This means a dynamical model that can predict how the concentrations of the metabolites of interest in the medium and the cell growth will evolve.

It has to be considered that only the external metabolites concentrations can be measured, known or controlled, so the inputs in the resulting model have to be variables only referring to concentrations of external metabolites.

So concretely what the model should be able to do is, given an initial values of the considered external metabolites concentrations, to predict the consumption/production rates of these metabolites, and then be able to simulate the dynamic behaviour of the metabolites along the time using this predicted values.

(9)

7

2. Initial description

Next a description of the information that is available before starting to apply any methodology or doing any calculations is explained. This available information consists basically in the equations that reach the metabolic stoichiometry, the parameters and conditions of the system that is wanted to model, a stoichiometric model of the metabolic network and the available data and the implicit information in it.

2.1. System equations description

The main elements of the metabolic network model, as said previously, are (1) the metabolites involved and their concentrations and (2) the reactions that produced the variation of these concentrations. The stoichiometric coefficients denote the proportion of substrate and product molecules involved in every reaction. For example, considering the reaction:

→ 2

The stoichiometric coefficients would be -1 for glucose and 2 for the lactate, which means that with one mole of glucose 2 of lactate can be produced. We could also say that to get one mole of lactate 0.5 mole of glucose is needed, so the stoichiometric coefficients could also be -0.5 for glucose and 1 for the lactate, the proportionality still is maintained. The sign indicates which metabolite is consumed and which is produced.

The consumption or production of a metabolite along time can be expressed using differential equations. For the previous example the equations would be the following:

( )

= − , and ^{( )}

= 2 (2-1, 2-2)

Where is the flux (or rate) of the reaction itself, but also the rate that glucose is consumed and the half of the rate that lactate is produced.

If we consider that the metabolites are produced or consumed in many reactions in the same system, the dynamics of the metabolic network can be represented by the following equation [2]:

=

Where is the concentration of the metabolite i, is the flux in the reaction j and is the stoichiometric coefficient of the metabolite i in the reaction j. In this statement it is assumed that the only cause of the change in the concentration is the mass flow due to the reactions.

The stoichiometric coefficients can be represented in a matrix where the rows are the metabolites and the columns the reactions. This matrix A is the so called stoichiometric matrix which defines the metabolic network stoichiometry and it’s the starting point of the study. So finally our system is described by the following:

(2-3)

(10)

8

⎝

⎜⎜

⎛ ⋮

⋮

⎠

⎟⎟

⎞=

⎣⎢

⎢⎢

⎢⎡ ^, ⋯ _,

⋮ ⋱ ⋮

, ⋯ ,

, ⋯ _,

⋮ ⋱ ⋮

, ⋯ , ⎦⎥⎥⎥⎥⎤

∙ ⋮

Where J is the number of reactions of the stoichiometric matrix, M is the number of external metabolites and N is the total number of metabolites (internal and external).

The modeling of the reaction rates as functions of the metabolites concentrations values will be considered in later steps due to the fact that this can vary a lot depending on the kind of reaction and the hypotheses that can be done.

2.1.1. The quasi steady-state assumption

In order to simplify the system an approximation on timescale separation is assumed, the quasi steady-state assumption [3]. This assumption is very common in this kind of systems. The reason to adopt this approximation is that the kinetics inside the cell is very much faster than the reactions which involve external metabolites. Metabolites that are changed faster by the reactions will reach the steady-state in a short period of time (for example, seconds), and for the others that are changed slower it will take much longer to reached the steady state (for example, hours or days). For the internal metabolites, after this short period when they reached the steady-state, the change in the concentration can be considered zero,

= 0. So the differential equation referring to the change in the concentration can be replaced for an algebraic equation.

For example, considering the pathway:

→ 6 → 2 The implicit differential equations are:

( )

= − , ^{( )}

= − and ^{( )}

= 2 (2-5, 2-6, 2-7) If the kinetic rate of is very much faster than the kinetic rate of ; then ^{( )}

can be approximated by 0 and 0 = − . The new equations are:

( )

= − , and ^{( )}

= 2 (2-8, 2-9)

And, following the law of Mass Action [4], the concentration of the fast changing metabolite is related algebraically to the Glucose ( is the kinetic rate of reaction ):

= and = 6 ; (2-10) 6 = ; (2-11)

Finally, by applying this, the system reads:

(2-4)

(11)

9

⎝

⎜⎜

⎛ ⋮ 0 ⋮

0 ⎠

⎟⎟

⎞=

⎣⎢

⎢⎢

⎢⎡ ^, ⋯ _,

⋮ ⋱ ⋮

, ⋯ _,

⋮ ⋱ ⋮

, ⋯ _, ⎦⎥⎥⎥⎥⎤

∙ ⋮

Where the internal metabolites consumption/production rates have been considered 0 because they change instantaneously when a variation is produced in the rates of the external metabolites.

2.2. The system object of the study

The system which is object of this study consists in a culture of CHO cells cultivated in batch mode in flasks which contain a medium with a determined composition. This means that cells will start growing with a determined initial medium concentration and a determined initial cell density (or initial state) until another state where these parameters will be measured again.

What is wanted to be modelled in this study is the dynamics of the change in these parameters.

The medium development has been done in parallel with this study. The objective of developing this medium is to achieve a medium where the cells can live properly and also its exact composition of all compounds is known in order to allow us to carry out some experiments just varying one or some of the concentrations in every one of them. These variations should induce the system to adopt different behaviours. Also has to be considered that the only compounds that can be measured during the experiments, or its concentration can be adjusted, are some extracellular compounds which are the ones that will be modelled.

These compounds are:

• The 20 amino acids, which are the main object of the study and are the ones which their initial concentrations will be changed.

• Glucose which is the main source of feeding in conventional cell cultures.

• Lactate and Ammonia; when the system begins to evolve the cells will consume part of these substrates producing the by-products lactate and ammonia which can also be consumed in a posterior state.

• The last variable that has to be considered is the cell density itself which is an important factor when it comes to determine the dynamics of the concentrations.

Another important point is that the phase of the cell life that is wanted to be modelled is the growth phase that is when the cells are increasing their density. So if in later studies it is desired to optimize the production of any compound or the cell growth itself, this should be the phase object of the study.

All these conditions must be represented by the metabolic network model that is being taken as a reference point to start building the new mathematical model.

(2-12)

(12)

10

2.3. The available network

A stoichiometric matrix from a metabolic network model was provided as a starting point from the KTH Division of Bioprocess. The stoichiometric matrix was based on the stoichiometric model proposed in [5] but adding the reaction of the biomass production which was calculated by the KTH Division of Bioprocess from the mass balances of the production of the main macromolecules (DNA, RNA, protein, carbohydrates and oleic acid). The coefficients corresponding to these equations are always estimated values due to the fact that they are referred to standard values suggested in other studies. Regarding this reaction the biomass itself can be considered as another external metabolite as long as it is linked to the other network compounds by stoichiometric relations.

In [5] the model was used to determine, from experimental data, the fluxes of the different pathways of the network and to see which pathways are followed when the Glucose feeding concentration was changed. The cells were cultivated in continuous culture. As described in 2.2 (The system object of the study), the system object of the study has different properties (the cells are cultivated in batch mode) but the stoichiometry relations must still be fulfilled.

The characteristics of the stoichiometric model are:

- 24 external metabolites, the ones that are in the system which is object of the study:

- 14 are only substrates: Glucose, Arginine, Histidine, Isoleucine, Leucine, Lysine, Methionine, Phenylanine, Threonine, Tryptophan, Valine, Asparagine, Serine and Proline. Regarding the exposed configuration for the reactions reversibility, Proline can also be a product, not only a substrate (see below).

- 8 are substrates and products: Ammonia, Glutamine, Glutamate, Alanine, Aspartate, Cysteine, Glycine and Tyrosine.

- 2 are only products: Lactate and the Biomass.

- 43 internal metabolites, among them, the main metabolites represented in the network diagram (see Appendix I), but also secondary metabolites.

- 46 reactions which describe the stoichiometry among all the network compounds.

- The reactions were considered to occur in the direction suggested in [5] excepting reactions number 5 and number 14 (see Appendix I) which have been chosen as reversible (Glutamine -> Glutamate and Proline -> Glutamate). The criterion that has been followed is that in articles [6], [7] and [8] where most of the same pathways of the same network are described, reaction number 5 goes in the other way and in article [7] reaction number 14 also goes in the other way (in articles [6] and [8]

reaction 14 is not considered as most of the amino acids).

The system defined by the stoichiometric matrix reads:

⎝

⎜⎜

⎜⎛

( )

⋮

( )

( ) 0⋮

0 ⎠

⎟⎟

⎟⎞

· =

⎣⎢

⎢⎢

⎡ ^, ⋯ _,

⋮ ⋱ ⋮

, ⋯ _,

, ⋯ ,

⋮ ⋱ ⋮

, ⋯ , ⎦⎥⎥⎥⎥⎥⎤

∙ ⋮

· (2-13)

(13)

11 Where:

- [nmol/MVC·day] is the specific consumption/production rate of the metabolite . The consumption/production rate is the rate at which one the cells are producing or consuming the metabolite , or, what is the same, the rate of the change in the concentration of the metabolite per million viable cells (MVC); thus = · (see below the definition of and ).

- ( ) is the specific biomass production or cell growth, the dimension determined by the stoichiometric coefficients is [10³ cells/ MVC·day].

- [nmol / MVC·day]is the flux value or specific rate of the reaction . Note that in the theoretic explanation in 2.1 (System equations) the rates where considered as global rates, not specific rates. The fluxes usually are expressed as a function of the concentrations:

= ( ), … , ( ), ( ), ( ), … , ( . ) - [μM] is the concentration of the metabolite . From to are the concentrations

referring to the extracellular metabolites, so they are medium concentrations and from to are concentrations referring to the intracellular metabolites, or what is the same, concentrations inside the cell. [10³ cells/ml] is the cell density.

- [MVC/ml] is the number of viable cells. Note that = · 10 , the denomination of it’s just for notation.

The complete stoichiometric matrix with all its coefficients, a schematic diagram of the main metabolic routes in the network and all the reactions that form the model are described in Appendix I.

2.4. The experimental data

In order to build a versatile model, representative experimental data is needed. This data must be obtained in different environmental conditions which represent different system states. In order to achieve this, several experiments have to be carried out with different initial conditions, or what is the same, the values of the initial concentrations of the external metabolites must vary in each experiment. Consequently, the variations in the initial conditions will induce the metabolic reaction network to adopt different behaviours so the consumption/production rates will also vary between experiments. Then, data which represents different system behaviours will be available. It also has to be considered that the variations in the initial concentrations have to be high enough to produce a significant change in the measurements among experiments.

A set of experiments was planned in order to represent the conditions specified above. The number of experiments, the amino acids that are varied and the level of variation referred to the pre-selected medium for each experiment are resumed in table 2.1.

In total there are 37 experiments: 2 using the reference medium, 2 decreasing the essential amino acids concentrations all together in two levels (a group representing the essential amino acids was created including Arginine, Histidine, Isoleucine, Leucine, Lysine, Methionine, Phenylanine, Threonine, Tryptophan and Valine), 22 changing the concentration of each non- essential amino acid one by one also in two levels and 11 removing one of the non-essential

(14)

12 amino acids in each. The distinction between the essential amino acids and the non-essential is done because they have different characteristics. The main treat of the essential amino acids is that they cannot be produced by the cells in any conditions, so they are an essential compound to maintain the cells living.

Varying Conditions Levels

0% -50% -80% -100%

Reference Medium 2 0 0 0

Essential Aminoacids 0 1 1 0

Ala 0 1 1 1

Asn 0 1 1 1

Asp 0 1 1 1

Cys 0 1 1 1

Gln 0 1 1 1

Glu 0 1 1 1

Gly 0 1 1 1

Pro 0 1 1 1

Ser 0 1 1 1

Tyr 0 1 1 1 Total

nº Experiments 2 11 11 10 34

Table 2. 1 Set of planned experiments

Next, a model of experiment methodology is exposed. The goal of the following model is to achieve various samples of data for the same initial conditions in an easy way and also to achieve data which let us compare the different behaviours of the system when the cell density grows up. Below, a summary of the methodology is exposed:

- First let the cells grow until a reasonable level where they should have a more stable behaviour. To achieve this, cells may grow between 1 and 3 days.

- Then, the cells are removed and they are introduced in to the flask filled with medium with the initial composition (the initial composition varies as it is indicated in table 2.1). In this point, all the concentrations and the cell density in this state of the system are known. The cells remain growing and consuming/producing metabolites in the medium during one day.

- Samples are taken and the concentrations and cell density are measured. The cells are diluted in the medium with the initial composition in order to achieve the cell density of the day before so the system is again in the same state.

- The process is repeated during 2 or 3 days in order to achieve samples in the same conditions.

- The same whole process should be repeated for different levels of cell density to get samples with different cell densities.

A schematic diagram of the cell density along the time is represented in figure 2.1.

(15)

13

Figure 2. 1 Cell density along the time during every experiment

The diagram for the other metabolites would have the same shape but always starting at the same level of concentration in every experiment and, of course, the concentration level during the cultivation process can be raised or downed depending on the metabolite. For example, for Arginine which is an essential amino acid and is always consumed, the concentration along the time would have a shape similar than the one represented in figure 2.2.

Figure 2. 2 Arginine concentration along the time during every experiment

2.4.1. Estimation of the values of the consumption/production rates and concentrations in a determined system state

From every sample the values of the consumption/production rate and concentration of a determined state of the system can be estimated for every external metabolite between two set of measurements in two different time instants.

First has to be considered an increase of the concentration of the metabolite during a period of time, where the values of the concentrations are known in the starting and the ending point. In our case the starting point is when the cells are introduced in the flask with medium

(16)

14 of known concentration and the ending point is when the samples are taken and the measurements are done after one day.

The method applied to estimate these parameters (consumption/production rates and concentrations for a determined system state) is based in the assumption that the value of the derivative when the concentration of the metabolite is in the middle point of the total increase can be approximated to the change of the concentration divided by the change in time. This is represented in figure 2.3.

Figure 2. 3 Estimation method to obtain and for a determined system state

Thus the specific consumption/production rate for the point where ( ) = (2-14) can be estimated by the following equation:

( )/ = −

∆ ·( + ) 2

This method is a simple linearization of the concentration curve suggested by the KTH Division of Bioprocess and is usually used for estimating the consumption rates when doing flux analyses. It is based on the fact that the kinetics for biological systems metabolites usually have an exponential growth but very smooth.

2.4.2. The available experimental data

The experiments that have been carried out differ a little from the planned ones. In reality there are sample of 12 different centrifuged tubes (experiments) which ones have different amino acids initial concentrations. The initial values of the concentrations of the centrifuged tube taken as reference (CT12) are shown in table 2.3 and the initial concentration variations respecting to CT12 of the other centrifuged tubes summarized in table 2.2.

There are available samples of the concentrations of the amino acids, glucose, ammonia and lactate of three different days for each of the centrifuged tubes (excepting Day 1 for CT04 and CT08):

(2-15)

(17)

15 Centrifugate Tube Varying Conditions Levels

0% -100%

CT12 Reference medium

concentrations 1 0

CT01 Ala 0 1

CT02 Asn 0 1

CT03 Asp 0 1

CT04 Cys 0 1

CT05 Gln 0 1

CT06 Glu 0 1

CT07 Gly 0 1

CT08 Pro 0 1

CT09 Ser 0 1

CT10 Tyr 0 1

CT11 Trp 0 1 Total

nº Experiments 1 11 12

Table 2. 2 Set of experiments with available data

• Day 1: concentrations after one day of cultivation.

• Day 3: concentrations after two days of cultivation but with the same initial concentrations but taking as initial cell density the diluted cell density of Day 1 ( · 4/5).

• Day 4: concentrations after one day of cultivation but taking as initial density the diluted cell density of Day 3 ( · 4/5).

Figure 2. 4 Cell density along the time of the experiments that have been carried out

The time diagram for the cell density equivalent to figure 2.1 but for the available samples is shown in figure 2.4.

From the samples listed in table 2.2 and for every day, the consumption/production rates have been calculated as described in 2.4.1 (Estimation of the values of the consumption/production rates and concentrations in a determined system state). So, as many sets of system states (

(18)

16 and ) as samples are obtained. The calculated values for and can be found in Appendix III.

One problem is that in the available experimental data the values of the concentrations of Serine and Asparagine are summed in the same value. The dealing with this problem will be exposed in 5.3.1 (The problem with Serine and Asparagine).

Concentrations in CT12 (µM)

Glucose 17500

Lactate 0

NH4 0

Glutamine 4000

Glutamate 258,1

Arginine 1765,3

Histidine 168,6

Isoleucine 455,5

Leucine 682,9

Lysine 542,8

Methionine 127,5

Phenylanine 199,3

Threonine 440,7

Tryptophan 45,2

Valine 607,7

Alanine 204,7

Aspargine+Serine 779,2

Aspartate 174,10

Cysteine 694,0

Glycine 252,1

Proline 402,7

Tyrosine 315,2

Cell density (MVC/ml) 0,5

Table 2. 3 Initial concentrations of the reference medium

(19)

17

3. The Modelling Method

The method presented here is the main object of the master thesis. The purpose of this methodology is to achieve a dynamical model where only the extracellular metabolites are involved an then simplify it achieving a model as much reduced as possible. The model should be able to predict the consumption/production rates of these metabolites given a determined concentration values.

In this method each cell is considered as a black box device which transforms the extracellular metabolites into other extracellular metabolites. The inputs of the mentioned device will be the extracellular substrates and the outputs the extracellular products. This is represented in figure 3.1 (Note that an extracellular substrate can be at the same time an extracellular product).

Figure 3. 1 The cell as a black box device

In order to calculate this reduced dynamical model, information consisting in two basic elements is needed:

• A pre-defined metabolic network model which should contain the reactions that are wanted to be considered in the system which implies basically the stoichiometric relations between all the compounds in the network or, what is the same, the stoichiometric matrix of the network model.

• Experimental data with samples obtained in different initial conditions describing the different system behaviours that are wanted to be modelled. The format of the data should be the one exposed in 2.4 (The available data); and .

The method consists of two parts:

• The design of the macroscopic model, where a stoichiometric matrix of macroscopic reactions is calculated from the stoichiometric matrix of the pre-defined metabolic network.

• The design of the reduced dynamical model, where a reduced dynamical model which relates the consumption/production of the external metabolites with their concentrations is calculated from the experimental data and the previously calculated stoichiometric matrix of the macroscopic reactions.

3.1. Design of the Macroscopic Stoichiometric Model

In this first part of the method a macroscopic model of the system is calculated by using elementary flux analysis. The complexity of the pre-defined stoichiometric matrix and the quantity of extracellular metabolites will determine the size of the resulting stoichiometric

(20)

18 matrix of the macroscopic model. The elementary flux analysis methodology has been also used in other papers such [6] or [8] where they are dealing with a similar problem and theory basis are exposed in [9]. Both articles ([6] and [8]) use networks much simpler than the one that is object of the study, so the whole processes used the articles cannot be extrapolated to the case described in 2.2 (The system object of the study). In order to illustrate the different steps of the first part of the method the network used in [8] has been taken as example.

3.1.1. System description for the macroscopic design

As is defined in 2.1 (System equations) the system can be represented by the following equations:

⎝

⎜⎜

⎛ ⋮ 0 ⋮

0 ⎠

⎟⎟

⎞=

⎣⎢

⎢⎢

⎢⎡ ^, ⋯ _,

⋮ ⋱ ⋮

, ⋯ _,

, ⋯ ,

⋮ ⋱ ⋮

, ⋯ _, ⎦⎥⎥⎥⎥⎤

∙ ⋮ = ·

The system can be divided in two subsystems which their metabolites are assumed to be submitted to different dynamic conditions: the internal metabolites which fulfil the steady- state conditions = 0, and the external metabolites which concentration change are not considered instantaneous. The stoichiometric matrix can also be represented as:

=

Where

= ^,

⋯ ,

⋮ ⋱ ⋮

, ⋯ _, and = ^,

⋯ ,

⋮ ⋱ ⋮

, ⋯ _,

For example, considering the network used in [8] (figure 3.2) where the external metabolites are Glucose, Glutamine, Lactate, Ammonia, the Biomass and the antibody IgG1, the stoichiometric matrixes , and are the ones represented in table 3.1 (the Biomass and the IgG1 are not metabolites themselves but are considered as ones due to that they are related to the other metabolites through stoichiometric equations).

(3-1)

(3-2)

(3-3, 3-4)

(21)

19

Figure 3. 2 Metabolic network of the system used as example [8]

Note that CO2 is not represented in the stoichiometric matrices. This is due to the fact that it is an extracellular metabolite that is only produced and its production rate is not object of the study. As long as it’s always a product it has no effect in the dynamics of the remaining metabolites. Thus it could be included in but the resulting elementary flux modes would be the same ones because they only depend on (see 3.1.2, The Kernel Matrix and the Elementary Flux Modes) and the flux of the macroscopic reactions only depend on the substrates (see 3.2.1, Reaction kinetics modelling).

r1 r2 r3 r4 r5 r6 r7 r8 r9 r10 r11 r12 r13 r14 r15 r16

A Aext

Glc -1 0 0 0 0 0 0 0 0 0 0 0 0 0 -0,0208 0

Lac 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0

NH4 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0

Gln 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 -0,0377 -0,0104

Ala 0 0 0 0 0 1 0 0 0 0 0 0 0 0 -0,0133 -0,0112

Biomass 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0

IgG1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

Aint

G6P 1 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0

DaP 0 1 -1 0 0 0 0 0 0 0 0 0 0 0 0 0

G3P 0 1 1 -1 0 0 0 0 0 0 0 0 0 0 0 0

Pyr 0 0 0 1 -1 -1 -1 0 0 0 0 1 0 0 0 0

AcoA 0 0 0 0 0 0 1 -1 0 0 0 0 0 0 0 0

Cit 0 0 0 0 0 0 0 1 -1 0 0 0 0 0 0 0

α-keto 0 0 0 0 0 1 0 0 1 -1 0 0 1 0 0 0

Mal 0 0 0 0 0 0 0 0 0 1 -1 -1 0 0 0 0

Oxa 0 0 0 0 0 0 0 -1 0 0 1 0 0 0 0 0

Glu 0 0 0 0 0 -1 0 0 0 0 0 0 -1 1 0 0

Table 3. 1 Stoichiometric Matrixes corresponding to figure 3.2

Considering the definitions above, the dynamics of the external metabolites can be described as:

(22)

20 ⋮ = ∙ = ∙ ∙

Where [nmol/MVC·day] is the vector of the specific consumption/production rates of the external metabolites, is the vector of the reaction rates [nmol/MVC·day], [MVC/ml] is the cell density and is the stoichiometric matrix of the external metabolites.

This subsystem can be simplified by only taking into account the specific values (divided by ) which makes it more functional to work with.

= ∙

The remaining subsystem for the intracellular compounds considering the quasi steady-state conditions reads:

= ∙ = 0

This means that all flux vectors that can take place when the internal metabolites have reached the steady state must fulfil this condition. Thus we could say that, for each internal metabolite, the sum of the fluxes in all the reactions which this metabolite takes part multiplied by the corresponding stoichiometric coefficient is zero or, what is the same, what is produced is also consumed.

3.1.2. The Kernel Matrix and the Elementary Flux Modes

The first step when it comes to apply the first part of the methodology is to calculate the elementary flux modes. A flux mode is a configuration of the reaction flux rates which fulfils the conditions of steady state for the internal metabolites. The overall idea is that the configuration of the fluxes rates is one that maintains the production of the concentration of the internal metabolites constant for a determined value of the consumption/production rates of the external metabolites. Thus by following the pathways described by the reactions which the flux rate is not 0 for this configuration (or flux mode) a link between the external metabolites is created and a macroscopic reaction is defined for every flux mode. Then, an elementary flux mode is, as its name says, a flux mode which cannot be decomposed in other flux modes, this concept is explained in detail below. In figure 3.3 examples of elementary flux modes for the network in figure 3.2 are represented. For example, the elementary flux mode represented by e1 in figure 3.3 creates a link between Glucose and Lactate. The rates of reactions 1, 2 and 3 are 1 and the reactions 4 and 5 ones are 2 (Defined by the stoichiometry and by the steady-state condition). The rest of flux rates in the network are 0. This configuration of the flux rates constitutes a flux mode which implicitly creates the macroscopic reaction:

→ 2

From the stoichiometric matrix of the internal metabolites , where usually the number of columns is larger than the number of rows (there are more reactions than metabolites), a null space can be defined. The dimension of this null space D will be the number of reactions

(3-5)

(3-6)

(3-7)

(23)

21 minus the rank of ( − ( )). Then D independent vectors with length can be defined which ones will form a basis for this null space. These vectors set in columns will form the so called Kernel Matrix.

Figure 3. 3. Elementary flux modes of the network in 3.2 [8]

Every set of fluxes that fulfil the condition of quasi steady-state can be expressed as a linear combination of the columns of the Kernel matrix (the vectors of the basis of the null space).

Note that the definition of the Kernel Matrix is not unique, there exists an infinite number of vectors which can form a basis of the null space. In our running example one possible Kernel Matrix would be the one represented in table 3.2.

( ) is 10 and there are 16 reactions, thus the number of Kernel vectors that form a basis of the null space D is 6.

k1 k2 k3 k4 k5 k6

r1 1 0 0 0 0 0

r2 1 0 0 0 0 0

r3 1 0 0 0 0 0

r4 2 0 0 0 0 0

r5 2 -1 1 0 0 0

r6 0 0 -1 1 0 0

r7 0 1 0 0 0 0

r8 0 1 0 0 0 0

r9 0 1 0 0 0 0

r10 0 1 0 1 0 0

r11 0 1 0 0 0 0

r12 0 0 0 1 0 0

r13 0 0 1 0 0 0

r14 0 0 0 1 0 0

r15 0 0 0 0 1 0

r16 0 0 0 0 0 1

Table 3. 2 The Kernel Matrix

(24)

22 If the flux of one reaction is always equal to 0 in all the vectors of the Kernel matrix, this means that the reaction corresponding to this flux is an equilibrium reaction and it will never take place under the quasi steady-state conditions (property 1); the flux of this reactions will be always zero when the intracellular metabolites have reached the steady state. And if the flux values of two or more reactions are the same (or a multiple) in all the vectors of the Kernel matrix, this means that these reactions form an unbranched reaction path; the value of the fluxes for these reactions will be the same (or a multiple) in all the sets of fluxes under the quasi steady-state conditions (property 2).

In the matrix above (table 3.2) it can be seen that reactions 1, 2, 3 and 4 always have the same entries for all the Kernel vectors (or a multiple). This means that it is an unbranched path. Thus reactions 1, 2, 3 and 4 could be simplified for the reaction → 2 (see figure 3.2).

In this case there are no entries in the Kernel vectors that are always zero, thus there are no equilibrium reactions and all the reactions can occur in quasi steady-state conditions.

We can now define a flux mode. A flux mode is a set of flux vectors defined as ( is the number of reactions, 16 in the running example):

= { ∈ | = · ^∗; > 0}

Where ^∗ fulfil two conditions:

• The steady-state condition; · ^∗= 0.

• The sign restriction condition defined by the reversibility of the reactions. This means that the entries ^∗ corresponding to a reaction which is irreversible has to be positive, so the reactions only occurs in the direction predefined in the stoichiometric matrix.

A flux mode is an Elementary Flux Mode when it cannot be decomposed as a positive linear combination of other flux modes. For example, considering the flux modes FM1, FM2 and FM3 represented in table 3.3 corresponding to the described in table 3.1 where all the reactions are considered irreversible.

FM1 1 1 1 2 2 1 0 0 0 1 0 1 0 1 0 0

FM2 1 1 1 2 2 0 0 0 0 0 0 0 0 0 0 0

FM3 0 0 0 0 0 1 0 0 0 1 0 1 0 1 0 0

Table 3. 3 Flux Modes vs. Elementary Flux Modes

FM1 is a flux mode because it fulfils the conditions exposed above but is not an elementary flux mode because it can be expressed as the sum of FM2 and FM3. Thus FM1 is a positive linear combination of FM2 and FM3.

The existent elementary modes for a determined stoichiometric matrix (also taking into account the reversibility of the reactions) are unique.

In order to calculate the elementary flux modes the software Metatool 5.1 has been used.

(3-8)

(25)

23 Metatool 5.1 is a software application which can be run in Matlab or GNU Octave and allows us to calculate the elementary flux modes. Just by giving as input the internal metabolites stoichiometric matrix and a vector with length the number of reactions containing logical values which indicates the reversibility of the reactions (1 if irreversible, 0 if reversible), a matrix where the rows correspond to the reactions and the columns to the uniquely defined elementary modes can be calculated. The Metatool syntax and description can be found in [10].

3.1.3. The Metatool procedure

By calculating a Kernel matrix for , the first step executed by Metatool is to build a reduced system using the properties 1 & 2 described above in the previous point (3.1.2). First, the columns of the main stoichiometric matrix corresponding to the blocked reactions are removed from the system, and also the rows corresponding to the metabolites that only are involved in these reactions (property 1), the blocked reactions are the equilibrium reactions which their entry is always 0. Using property 2, subsets of reactions that represent the unbranched paths are created.

In the running example there are no blocked reactions and the subsets corresponding to the unbranched paths are:

• Reactions 1, 2, 3 and 4.

• Reactions 7, 8, 9 and 11.

• Reactions 12 and 14.

In the reduced system, these subsets will be represented as a linear combination of their reactions. These linear combinations reduce the number of metabolites, creating zeros in the rows corresponding to these metabolites which are in the middle of the unbranched paths.

The space transformation carried out by Metatool that transforms the original system to the reduced system can be checked in the calculated matrix “sub”. The matrix of the subsets for the network in figure 3.2 is shown in table 3.4.

r'1 1 1 1 2 0 0 0 0 0 0 0 0 0 0 0 0

r'2 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0

r'3 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0

r'4 0 0 0 0 0 0 1 1 1 0 1 0 0 0 0 0

r'5 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0

r'6 0 0 0 0 0 0 0 0 0 0 0 1 0 1 0 0

r'7 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0

r'8 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0

r'9 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1

Table 3. 4 Space transformation: Original system->Reduced system

So, for example, the reaction 1 of the reduced system r'1 is the sum of r1, r2, r3 and 2 times r4.

By looking at (table 3.1) it can be seen that when this linear combination is done the rows corresponding to G6P, DaP and G3P are zero, thus these metabolites are removed from the system in order to reduce it. The remaining reaction is → 2 . The other resulting

(26)

24 equations from the subsets are + → + 2 and + → + + . The last subset doesn’t affect to the metabolites reduction but reduce the system by one more reaction.

The metabolites remaining in the reduced system are registered in the vector “^rd_met” (table 3.5).

rd_met 4 7 8 10

Table 3. 5 rd_met; the remaining metabolites in the reduced system

So the reduced system only will have the metabolites number 4, 7, 8 and 10 in which are, respectively, Pyr, α-keto, Mal and Glu.

Applying the space transformation specified in table 3.4 and removing the rows that only contain zeros, the reduced system (table 3.6) reads:

r'1 r'2 r'3 r'4 r'5 r'6 r'7 r'8 r'9 Pyr 2 -1 -1 -1 0 1 0 0 0 α-keto 0 0 1 1 -1 0 1 0 0

Mal 0 0 0 -1 1 -1 0 0 0

Glu 0 0 -1 0 0 1 -1 0 0

Table 3. 6 The reduced system

The next step after checking the linear correlation coefficient of the matrix that Metatool executes, is to check which metabolites in are only produced or consumed and if anyone of them takes part in any reversible reaction. If it exists any metabolite that is only produced or consumed (all the stoichiometric coefficients are positive or negative) and it doesn’t take part in any reversible reaction, automatically, all the reactions where this metabolite is involved will never occur in steady state conditions because the only solution of the equation below (the steady-state condition for the metabolite ) is | = 0. Where | is a vector containing the fluxes which corresponds to the reactions where the stoichiometric coefficient of the metabolite , , is different than 0.

, = [ ⋯ ] · = 0

The same procedure is carried out after removing the blocked reactions due to property 1.

For the stoichiometric matrix of the running example in table 3.2 the result of Metatool is the following:

0 metabolites are only produced, 0 are only consumed;

0 metabolites take part in only one reversible reaction; 0 are unused.

Removing 0 blocked reactions from subsets

0 metabolites are only produced, 0 are only consumed;

0 metabolites take part in only one reversible reaction; 0 are unused.

In this case, there are no blocked reactions and no metabolites that are only produced or consumed. Thus there are no reactions that won’t take place due to the existence of metabolites that are only produced or consumed.

(3-9)

(27)

25 Based on the above, Metatool starts doing combinations with the Kernel vectors that fulfil the condition of elementary mode and the irreversibility conditions. Then the elementary flux modes of the reduced model are calculated in the matrix “rd_ems”, let us call it . The resulting matrix for the running example is represented in table 3.7. In this case there are 7 elementary flux modes.

e1 e2 e3 e4 e5 e6 e7

r'1 1 1 0 0 0 0 0

r'2 2 0 1 0 0 0 0

r'3 0 0 0 1 0 0 0

r'4 0 2 0 0 1 0 0

r'5 0 2 1 1 2 0 0

r'6 0 0 1 1 1 0 0

r'7 0 0 1 0 1 0 0

r'8 0 0 0 0 0 1 0

r'9 0 0 0 0 0 0 1

Table 3. 7 Elementary Flux Modes of the reduced system

The next step is to calculate the matrix containing the elementary flux modes for the original reactions . This step is not done automatically by Metatool, but can be calculated just by undoing the space transformation from the subsets of reactions to the original reactions (multiplying by the subsets matrix transposed; in table 3.4):

= ·

Finally the matrix (table 3.8) for the system in figure 3.2 reads:

e1 e2 e3 e4 e5 e6 e7

r1 1 1 0 0 0 0 0

r2 1 1 0 0 0 0 0

r3 1 1 0 0 0 0 0

r4 2 2 0 0 0 0 0

r5 2 0 1 0 0 0 0

r6 0 0 0 1 0 0 0

r7 0 2 0 0 1 0 0

r8 0 2 0 0 1 0 0

r9 0 2 0 0 1 0 0

r10 0 2 1 1 2 0 0

r11 0 2 0 0 1 0 0

r12 0 0 1 1 1 0 0

r13 0 0 1 0 1 0 0

r14 0 0 1 1 1 0 0

r15 0 0 0 0 0 1 0

r16 0 0 0 0 0 0 1

Table 3. 8 Elementary Flux Modes

The calculated elementary flux modes from e1 to e5 are the ones represented graphically in figure 3.3. The elementary flux modes e6 and e7 are the corresponding to the reactions that

(3-10)

(28)

26 don’t involve intracellular metabolites which are the Biomass production from Glucose and the considered amino acids and the IgG1 production from the considered amino acids.

3.1.3.1. Reactions that don’t involve internal metabolites

The reactions that don’t involve internal metabolites can be treated independently of the elementary flux modes calculation as long as they are independent of . These reactions are represented in the stoichiometric matrix of the internal metabolites as columns containing only zeros and Metatool automatically calculates an elementary flux mode for each one of these reactions which contains a 1 in the position corresponding to this reaction (see e6 and e7 in table 3.8). It has been found that when these reactions are considered reversible, Metatool don’t calculate the elementary flux mode corresponding to the reversible reaction, what means an elementary flux mode containing a -1 in the position corresponding to the reaction. In order to solve this problem the methodology that has been adopted is to treat independently these reactions just by adding them and the corresponding reversible reactions to the stoichiometric matrix of the macroscopic reactions defined in the next point 3.1.4. This has been taken into account in the implementation process (4.1).

3.1.4. Calculation of the macroscopic reactions stoichiometry

Once all the elementary flux modes are calculated we can define the matrix where the columns of it represent the elementary flux modes and it would have the following shape:

= ^,

⋯ _,

⋮ ⋱ ⋮

, ⋯ ,

Where _, is the value of the flux rate of the reaction in the elementary mode and N is the number of elementary modes which only depends on the characteristics of . For the running example, J would be equal to 16 (there are 16 reactions) and N is equal to 7. According to the definition the matrix will fulfil:

∙ = 0

Then we can create a link between the external metabolites by considering the steady-state in the system. Here, what is considered is that, in quasi steady-state, the only possible values for the fluxes of all the reactions are a linear combination of the elementary flux modes. So the stoichiometric matrix of the macroscopic reactions can be defined as:

= ∙ = ′_, ⋯ ′_,

⋮ ⋱ ⋮

′_, ⋯ ′_,

Every macroscopic reaction corresponds to one elementary flux modes. We could say that a macroscopic reaction is the flux of an elementary flux mode but seen from an extracellular point of view.

In the example exposed previously, the values of the matrix entries are represented in table 3.9.

(3-11)

(3-12)

(3-13)

Dynamical Modelling of a Metabolic Reaction Network

Dynamical Modelling of a Metabolic Reaction Network

JOAN GONZALEZ HOSTA

Masters’ Degree Project Stockholm, Sweden Sep 2009

Dynamical modelling of a metabolic reaction network

Master Thesis in Automatic Control, KTH

Joan Gonzalez Hosta

Supervisor: Elling W. Jacobsen

External contact: Veronique Chotteau

Stockholm-Barcelona, 2009

Abstract

Table of contents

1. Introduction

1.1. Basic concepts

1.2. The aim of the master thesis

2. Initial description

2.1. System equations description

2.2. The system object of the study

2.3. The available network

2.4. The experimental data

3. The Modelling Method

3.1. Design of the Macroscopic Stoichiometric Model