Evaluation and Development of Methods for Identification of Biochemical Networks

(1)

Department of Physics and Measurement Technology

Master’s Thesis

Evaluation and Development of Methods for

Identification of Biochemical Networks

Alexandra Jauhiainen

LITH-IFM-EX-05/1378-SE

Department of Physics and Measurement Technology Link¨opings universitet

(2)

(3)

Master’s Thesis LITH-IFM-EX-05/1378-SE

Evaluation and Development of Methods for

Identification of Biochemical Networks

Alexandra Jauhiainen

Supervisor: Mats Jirstrand

Fraunhofer-Chalmers Research Centre for Industrial Mathematics

Examiner: Jesper Tegn´er

IFM, Link¨opings universitet

(4)

(5)

Avdelning, Institution Division, Department

Division of Computational Biology

Department of Physics and Measurement Technology Link¨opings universitet

SE-581 83 Link¨oping, Sweden

Datum Date 2005-02-22 Spr˚ak Language Svenska/Swedish Engelska/English Rapporttyp Report category Licentiatavhandling Examensarbete C-uppsats D-uppsats ¨Ovrig rapport

URL f¨or elektronisk version http://www.ep.liu.se/exjobb/ifm/bi/ 2005/1378/ ISBN — ISRN LITH-IFM-EX-05/1378-SE Serietitel och serienummer Title of series, numbering

ISSN —

Titel Title

Evaluering och utveckling av metoder för identifiering av biokemiska nätverk Evaluation and Development of Methods for Identification of Biochemical Net-works Författare Author Alexandra Jauhiainen Sammanfattning Abstract

Systems biology is an area concerned with understanding biology on a systems level, where structure and dynamics of the system is in focus. Knowledge about structure and dynamics of biological systems is fundamental information about cells and interactions within cells and also play an increasingly important role in medical applications.

System identification deals with the problem of constructing a model of a tem from data and an extensive theory of particularly identification of linear sys-tems exists.

This is a master thesis in systems biology treating identification of biochemical systems. Methods based on both local parameter perturbation data and time series data have been tested and evaluated in silico.

The advantage of local parameter perturbation data methods proved to be that they demand less complex data, but the drawbacks are the reduced information content of this data and sensitivity to noise. Methods employing time series data are generally more robust to noise but the lack of available data limits the use of these methods.

The work has been conducted at the Fraunhofer-Chalmers Research Centre for Industrial Mathematics in G¨oteborg, and at the division of Computational Biology at the Department of Physics and Measurement Technology, Biology, and Chemistry at Link¨oping University during the autumn of 2004.

Nyckelord

(6)

(7)

Abstract

Systems biology is an area concerned with understanding biology on a systems level, where structure and dynamics of the system is in focus. Knowledge about structure and dynamics of biological systems is fundamental information about cells and interactions within cells and also play an increasingly important role in medical applications.

System identification deals with the problem of constructing a model of a system from data and an extensive theory of particularly identification of linear systems exists.

This is a master thesis in systems biology treating identification of biochemical systems. Methods based on both local parameter perturbation data and time series data have been tested and evaluated in silico.

The advantage of local parameter perturbation data methods proved to be that they demand less complex data, but the drawbacks are the reduced information content of this data and sensitivity to noise. Methods employing time series data are generally more robust to noise but the lack of available data limits the use of these methods.

The work has been conducted at the Fraunhofer-Chalmers Research Centre for In-dustrial Mathematics in G¨oteborg, and at the division of Computational Biology at the Department of Physics and Measurement Technology, Biology, and Chemistry at Link¨oping University during the autumn of 2004.

Keywords: Systems Biology, System Identification, Biochemical Networks

(8)

(9)

Acknowledgement

I would like to thank my supervisor at Fraunhofer-Chalmers Research Centre (FCC), Mats Jirstrand, for his help and enthusiasm in this project. Additional thanks to my friends and everyone who has provided comments and support on my thesis work. Final thanks to the staff at FCC.

(10)

(11)

Notation

Symbols and abbreviations used in the thesis are gathered here for clarification. All the abbreviations are also explained in the main part of the thesis, at their first occurrence.

Symbols

x, X Boldface letters are used for vectors, matrices, and sets. θ Parameter vector.

D_M Set of values over which θ ranges in a model structure.

Abbreviations

SITB System Identification ToolBox MAPK Mitogen Activated Protein Kinase PRBS Pseudo Random Binary Signal

(12)

(13)

Introduction

The issues examined in this master thesis work are introduced in this chapter, together with the structure and aim of the thesis. Some background information on systems biology is also given.

1.1 Problem Formulation

The main task of this master thesis is to examine and evaluate different methods for reconstruction and identification of biochemical networks. The methods are also investigated with an aim to possibly improve their applicability.

1.2 Systems Biology

The task of the master thesis is within the area of systems biology. This section aims to explain different interpretations of the systems biology concept as well as give some motivation on why we perform research on this subject.

1.2.1 What is Systems Biology?

Several different opinions on how to define systems biology exist. One description of the area is as a field of research concerned with understanding biology on a systems level. A systems biologist is interested in understanding the structure and dynamics of a system (Kitano, 2002). In more detail, this understanding of a system can be divided into four parts: system structure, system dynamics, control

(18)

2 Introduction

methods, and design methods. The approach on how to thoroughly apprehend biological systems is proposed by (Kitano, 2001):

System structure Networks of biochemical character is to be identified. This means that signal transduction as well as controlling mechanisms and mass flow between entities of the system need to be recognised.

System dynamics When fundamental knowledge about the structure of a system is found, it is of interest to learn more about system behaviour over time and under different conditions.

Control methods Methods on how to control the system can be the next step after identification and analysis of dynamic activities. The interest might be how to control or sustain a state for the system.

Design methods Design principles and simulations with information on strate-gies to modify the system. The aim is to achieve the ability to modify, or even create, a system to fulfill certain purposes; for example provide cures for diseases.

The systems biology research hence strives to understand the complex interactions between DNA, RNA, proteins, metabolic and informational pathways. Both inter-and intra-cellular networks are of interest.

The application of systems and control theory to biological systems is one of many definitions of systems biology (Wolkenhauer et al., 2003). This view can be illus-trated with the following example from (Yi et al., 2000).

Example 1 Robust perfect adaption

Adaption, or desensitisation, is a process of biological systems, that allows the system to return to a steady state, despite continuously being subject to stimulating signals. An example of this is the adaption in bacterial chemotaxis to changes in the levels of constituent proteins. The process of adaption is robust and a result of integral feedback control.

The process is exemplified with the model of Figure 1.1. The amount of the species Y is held at a constant level by the integral feedback loop. The amount is only dependent on the rate constants for the rates leading to and from the A-component as follows: dA dt = V3Y K3m+ Y − V4 Y = V4K3m V3− V4

(19)

1.2 Systems Biology 3 Activate Inhibit Y A v3 v4 E2 E1 (saturation)

Figure 1.1. The process of robust perfect adaption. E1 and E2 represent perturbations to the system.

The rate v4must be constant, i.e., operating at saturation as the figure indicates.

The interpretation of the rate constants and a more thorough treatment on kinetics can be found in Section 2.5. The block representation of this system is shown in Figure 1.2, where the reference signal r is the level of the component that the systems strives to maintain constant:

G(s) Σ + -1 s r e u y

Figure 1.2. Robust perfect adaption as a block representation of integral feedback

control.

In addition to systems and control theory, mathematical and computational tools are utilised in research, such as data preprocessing, statistical and informatics mining tools (Morel et al., 2004).

The area of systems biology is not new but has evolved concurrently with the development of suitable tools and increasing experimental skills and technology in producing useful data. Bioinformaticians have large pools of information owing to several genome projects producing sequence data. In contrast, a systems biologist is in need of a different kind of data, and to retrieve this type of data is not an easy task (Wolkenhauer et al., 2003).

(20)

4 Introduction

With the definition of systems biology as the application of system theory to biol-ogy, the lack of applicable data can be understood. In systems and control theory, perturbations, performed in a systematical manner, and the following dynamic response of the system are used to deduce facts and information about the system.

If this approach is copied to biology, time series data of the biological system is required. These kinds of measurements are preferably performed in vivo, to ensure that the system is in its natural environment, and that is not an easy task. In addi-tion, data need to be collected systematically and, as in all applications, preferably with small influence of noise, which is hard to achieve in biological systems. It is also difficult to perform measurements without altering the environment or the state of the system under observation. The available data is more often steady state data, which is a lot less informative since it does not contain any information about the dynamic behaviour of the system.

The systems biology approach to a problem can be illustrated with the following example of different levels of modelling of a enzyme catalysed transformation.

Example 2 Modelling of an enzymatic transformation

The three figures gathered below in Figure 1.3 illustrate the modelling of a transfor-mation catalysed by an enzyme. The concept is further explained in Section 2.5.3. The first of the figures is a simple view of the transformation and gives the infor-mation that an enzyme catalyses the transforinfor-mation of the substrate, S, into the product, P . The second figure is a more detailed description of how the transfor-mation occurs. The enzyme and substrate form a complex which dissociates into product and free enzyme.

S P

E

(a) Textbook drawing.

S E

ES

P

(21)

1.2 Systems Biology 5 S E k₁/k_-1 k₂ ES P

(c) Circuit diagram, modified after (Wolkenhauer et al., 2003).

Figure 1.3. Three different levels of modelling of an enzyme catalysed reaction. The substrate is denoted S, the product P , and the enzyme-substrate complex ES.

The last of the figures is an elaborate view of this kinetic reaction:

S + E k1 k−1 ES k2 → E + P

A mathematical model can be built from the kinetic interactions in the circuit diagram. The diagram is a mapping of the four differential equations which describe the dynamics of the reaction and the equations can be used in simulation of the system (Wolkenhauer et al., 2003).

A circuit diagram is useful in a systems biology approach. It can be used as a module to build models of larger networks with these kind of kinetics. On the other hand, this could also be the kind of information a systems biologist would like to obtain for a network with less known kinetics.

(22)

6 Introduction

1.2.2 Why Perform Research on Systems Biology?

Why do we need information about the structure and dynamics of a system? Sev-eral reasons exist. Information provided by systems biology research is a part of fundamental knowledge about cells and interactions within cells. The information can also, for example, be used in medical applications.

A vast number of diseases are not dependent on a singly mutated gene or, for example, a misfolded protein. Instead, complex multimolecular interactions explain the disease state (Morel et al., 2004). Therefore, the system of interactions needs to be understood in order to explain the disease mechanism.

In addition, knowledge about the controlling mechanisms in a system can help in identification of important nodes in the network, which are possible targets for drug action or other treatments (Kitano, 2002).

A recent example of the utility of systems biology in medicine is a research project on the possibility of replacing animal testing with computer modelling called BioSim (BioSim Network, 2005). The aim of the project is to decrease the number of, or at least improve, animal tests and speed up the drug development process. The com-putational modelling will hopefully give more detailed information on how drugs affect the patient and how the drugs are processed in the body.

1.3 The Thesis

1.3.1 Aim and Audience

The aim of this thesis is to describe system identification methods applied to bio-chemical systems and the different problems arising when doing so. The report is meant to give the necessary background information for the reader to understand how the work has been carried out and what results this has led to.

The thesis is written with a person in mind having basic knowledge in biology as well as familiarity with topics like mathematics, control theory, and signal processing.

1.3.2 Structure

The thesis is commenced with an abstract summing the work in a few lines while acknowledgements and notation clarifications are written there after. As a guide

(23)

1.3 The Thesis 7

to the thesis, a comprehensive table of contents as well as a list of figures and a list of tables are listed.

The introductory chapter is followed by a more extensive theory part meant to explain fundamental concepts needed for further reading. Ensuing the theory block is a method chapter where the different methods used in the thesis are demonstrated and analysed.

The results are presented and discussed in the consecutive chapters and the main part of the thesis is completed with a chapter containing the conclusions from the work. Finally, a reference list and appendices are given.

(24)

(25)

Chapter 2

Theory

This chapter introduces necessary theory in order for the reader to understand the different methods used in the work. The methods are further described in Chapter 3.

2.1 Systems, Modelling, and Simulation

Basic concepts in modelling, simulation, and validation are explained in this section. The view of system identification here is the one used in traditional engineering applications. The information can be valuable to compare to the system identifi-cation methods that have sofar been utilised and developed in systems biology, see for example (Wahde and Hertz, 2000, Ideker et al., 2001, Friedman et al., 2000). The information in this section is collected from (Ljung and Glad, 2004), if not otherwise stated.

2.1.1 Systems

A broad definition of a system is a group of objects, or possibly a single object, that we wish to study. Another one, not including information about our interest in the system, is that the system is a group of elements interacting and bringing about cause-and-effect relationships within the system (Close et al., 2002).

A system can be studied by performing experiments. The observable signals are usually referred to as outputs. Signals that in turn affect the system are called inputs.

(26)

10 Theory

Sometimes, it is not possible to perform all the experiments on the system one would like due to, safety, financial, time consumption, or technical reasons. A model of the system is needed if it is to be examined.

2.1.2 Models

A model is a tool, describing the system of interest, that enables us to answer ques-tions about the system without performing experiments. The model is a description of how the elements in the system relate to each other (Ljung, 1987).

How detailed a model is depends on what level we wish to examine the system. What level we choose is dependent on the purpose of the model.

The interest in modelling is not always focused on the details of the dynamics within the system. The relation between the input and output signals can instead be the main focus, and less attention is put on the interpretation of the detailed internal dynamics. A model with this purpose is called a black box model. The parameters in the model are estimated with the only purpose of connecting the output and input signals, and are not associated with physical properties (Ljung, 1987). Sets of standard models that have parameters possible to relate to physical properties exist. These mixed models are referred to as grey box models.

It is important to remember that a model is simply a model. The possibility of an exhaustive description of the behaviour of a system is non-existent. In addition, observing the system always results in different kinds of noise, stochastic in nature, bringing about unpredictable variations in the observations.

In engineering applications, the most extensively used models are of mathemat-ical nature. A mathematmathemat-ical model defines, through distinct equations, how the elements of a system relate to each other.

A mathematical model can be time discrete or continuous, deterministic or stochas-tic, linear or nonlinear, and lumped or distributed. What attributes to assign to a model is dependent on the system and what information we can retrieve from it.

A simple, time continuous, lumped model of a biological system is given in the example below.

Example 1 The Lotka-Volterra equations

d

dtN1(t) = aN1(t) − bN1(t)N2(t) d

(27)

2.2 System Identification 11

The equations describe the dynamics between the populations of prey, N1(t), and

predator, N2(t), respectively. The variables can represent population density or

biomass. Even though this model is simple it is highly nonlinear and cannot be examined analytically.

To form the equations, several assumptions have been made regarding growth, death, and predation. The assumptions, analysis of the system, and extensions of the model can be found in (Edelstein-Keshet, 1988). The equations are originally from a classical article by Volterra.

2.1.3 Simulation

A model can be used to deduce the behaviour of a system under certain conditions, which correspond to experimental tests. The deduction of results can be performed by analytical computation or by performing numerical calculations using the model. The latter is what we call simulation. If we think of a model as a set of instructions, for example equations in the case of mathematical models, a simulator obeys these instructions and generates the behaviour (Zeigler et al., 2000).

2.2 System Identification

System identification is about constructing a model of a system from observed data. The information in this section is a summary of the introductory chapters of (Ljung, 1987).

2.2.1 The Identification Process

The system identification process can be divided into three phases: recording data, obtaining a set of candidate models, and choosing the best of the candidate models.

Recording data Input and output signals are monitored and recorded during experiments. The objective is to generate data carrying as much information as possible about the system. Special identification measurements can be made or data is gathered during normal operation of the system.

Obtaining a set of candidate models This step is the most crucial and diffi-cult one in the identification procedure. A set of models have to be chosen on the basis of the current knowledge of the system. One can choose from sets

(28)

12 Theory

of basic models, or a model can be developed using physical characteristics of the system. A model from a set of standard models is a black box or a grey box model.

Choosing the best of the candidate models according to data In this phase, the actual identification occurs. Parameters in the chosen model are esti-mated from the data.

2.2.2 Validation

Assume we have obtained a model of a system from our recorded data. Can we trust the model? We have to examine if the model is useful for the purpose we intended it to. This is what is more known as validation.

A model is generally considered to be useful, if it can reproduce experimental data. We compare the behaviour of the model with the system behaviour (Ljung and Glad, 2004). Important to remember is that a model is only valid within its validated area. It is never advisory to extrapolate information from a model.

What if we are not satisfied with the result that our model produced? Then, we need to go back to our identification process and change one or more of the criteria. Going through the process again will produce a different, and perhaps better result.

2.3 Biological Network Structures

Different types of biological networks exist. In this thesis, the emphasis is on the molecular level: intra-cellular networks and components. An account of different types of intra-cellular networks are given in this section, in order for the reader to understand the kind of networks that the identification methods in the thesis may be applied to.

An arbitrary network within a cell is a biochemical network, since it consists of entities in a biological unit, the cell, and they interact with some kind of chemical reaction. The grouping into separate kinds of networks differ in various articles and textbooks. Biochemical networks are sometimes not considered to include the genetic networks of the cell. However, the networks of the cell are divided into metabolic, gene regulatory, and signalling networks in this thesis.

2.3.1 Metabolic Networks

The majority of biochemical reactions are degradative or synthetic (Zubay et al., 1995). The synthetic reactions, also called anabolic reactions, are the basis for

(29)

2.3 Biological Network Structures 13

the formation of large biomolecules from simpler units. An example is the assem-bly of proteins from nucleic acids. The catabolic reactions are degradative and involves the breakdown of larger, complex organic molecules into smaller compo-nents (Zubay et al., 1995). The β-oxidation of fatty acids is a catabolic process. The catabolic and anabolic processes comprise the metabolism of the cell.

A metabolic pathway is a complex series of reactions (Zubay et al., 1995) and the word series emphasises the directional property of the pathways; there is a net flow of mass in a specific direction or a specific end purpose of the pathway, for instance one, or several, products. A pathway can be considered as a network with a relatively small amount of interconnections. When different pathways share com-ponents, for example second messengers, networks of higher order and complexity occur.

In describing a metabolic network, one usually distinguishes biochemical connec-tions from controlling connecconnec-tions. Controlling links represent that two nodes are connected without any mass transfer. Instead one species influences the production or consumption rate of the other species.

The controlling links are not always represented in traditional pathway descriptions, although they are important. Controlling links are part of feedback and feedforward loops that are important for the stability and robustness of the pathway.

2.3.2 Gene Regulatory Networks

The central dogma of molecular biology states that DNA is replicated, transcribed into RNA, and translated into proteins. DNA interacts with a vast set of molecules in the cell, from complex proteins to simple transcription factors.

The RNA-molecules usually make up the nodes of a gene regulatory network. Mass flow is not as frequent in these networks as in a metabolic network and instead controlling mechanisms dominate.

2.3.3 Signalling Networks

Cells receive, process, and respond to information in a process called signal trans-duction (Lodish et al., 2000). The signals are mediated by signalling pathways, and since components often interact between the pathways, they form networks (Bhalla and Iyengar, 1999). Mechanisms of information transfer might be protein-protein interactions, enzyme activity regulations, or phosphorylation. The last form of transfer is exemplified below using the highly conserved kinase cascade where Mi-togen Activated Protein Kinase (MAPK) is activated (Lodish et al., 2000).

(30)

14 Theory

Example 2 MAPK Cascade

The cascade is built up by several levels where a kinase in the upstream level phos-phorolyses a kinase in the level downstream. A MAPK cascade with all controlling elements included is given in Figure 2.1.

v4 MKKK MKKK-P MKKK-PP v1 v3 v2 v8 MKK MKK-P MKK-PP v5 v7 v6 v12

MAPK MAPK-P MAPK-PP

v9 v11 v10 Ras/MKKKK Activate Inhibit Intramodular interactions

Figure 2.1. A MAPK cascade.

2.4 Identification of Biological Networks

System structure identification corresponds to the first of the four parts of under-standing a biological system described in Section 1.2.1. In this thesis, the focus is on network structure identification. Network structure identification is not as broad a concept as systems structure identification and excludes for example the identification of structural connections among cells as well as cell-cell association (Kitano, 2001).

(31)

2.5 Chemical Kinetics 15

2.4.1 Identification Approaches

The identification of a network includes finding all the components, their function, and how they interact. This is a difficult task since this kind of information can-not be inferred from experimental data based on some general rules or principles. Biological systems are stochastic in nature and not necessarily optimal (Kitano, 2001).

In addition, several network realisations might produce similar experimental data and the identification involves singling out the correct one. This corresponds to the process of finding the best model out of a set of candidate models as described in Section 2.2.1.

There are two general approaches in network structure identification: bottom-up and top-down identification.

Bottom-up Different sources of data, for example literature and experiments, are integrated in order to understand the network of interest. This data-driven approach is mostly suitable when almost all pieces of the network already are known and the quest is to find the missing parts.

Top-down A more hypothesis-driven approach where high-throughput data is utilised in trying to determine the network structure. Some information about the network is usually needed before hand, but not as extensive as in the bottom-up approach.

2.5 Chemical Kinetics

Chemical kinetics is the study of rates of chemical reactions under consideration of all the intermediates in the process (Atkins and Jones, 1999). The area also examines the details about how reactions advance and what determines the rate of a chemical reaction.

2.5.1 Rate Laws and Reaction Mechanisms

Dinitrogen pentoxide, N2O5, is an organic compound present in solid form at room

temperature. The chemical structure of the substance is depicted in Figure 2.2 with the covalent bonds linking the different atoms in the molecule. More information on the substance can, for example, be found in (Linstrom and Mallard, 2003).

(32)

16 Theory O O N O O N O Figure 2.2. Chemical structure of dinitrogen pentoxide.

nitrogen oxide and oxygen according to

2N2O5 x1 (g) → 4NO x2 (g) + 3O2 x3 (g) (2.1)

The variables under each substance denote the amount of the substance in question. A plot of data for the reaction is shown in Figure 2.3.

0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Time (min)

Amount of dinitrogen pentoxide (M)

Figure 2.3. Decomposition data for dinitrogen pentoxide.

From this graph, one can observe that the amount of N2O5 is consumed, fast

initially, but the consumption rate decreases gradually. If we make a plot of the logarithm of the data, and plot it against the time as before, we get the plot in Figure 2.4 of a straight line. This plot confirms the fact that the reaction is of first order, which means that the amount of dinitrogen pentoxide is consumed with a rate directly proportional to its amount.

(33)

2.5 Chemical Kinetics 17 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 −1.8 −1.6 −1.4 −1.2 −1 −0.8 −0.6 −0.4 −0.2 0 Time (min)

Logarithm of amount of dinitrogen pentoxide (M)

Figure 2.4. Logarithmic plot of decomposition data.

amount of each compound is as before denoted xi and the rate is denoted r.

dx1 dt =−2r =−2kx1 (2.2) dx2 dt = 4r = 4kx1 (2.3) dx3 dt = 3r = 3kx1 (2.4)

The rate constant for the decomposition step of equation (2.2) is defined as the proportionality constant (excluding the sign), and is hence 2k. It can be deduced from the experiments as the slope of the straight line in Figure 2.4. Observe that the formation parts of the reaction, equations (2.3) and (2.4), have other rate constants, since the stoichiometric ratios between the species differ.

The order of a reaction cannot in general be determined from the reaction formula; it is a property determined by experiments. The order varies from one chemical reaction to another and fractional orders also exist (Atkins and Jones, 1999).

If two species react to form a third species, the overall order is defined as the sum of the orders for each reactant. For the example below, the order is a + b.

A + B → C gives dA

dt = −k(xA)

a_(x B)b

The difficulty of extracting the rate law for a reaction directly from its equation owes to the fact that all, but the simplest, reactions occur in several steps called elementary reactions (Atkins and Jones, 1999). All steps might not be given in a

(34)

18 Theory

reaction formula. To understand how the reaction proceeds, a mechanism needs to be proposed.

Assume that we have a total reaction A + 2B → C + D. A possible reaction mechanism is given below, where X and Y represent intermediates.

A + B → X r1= k1xAxB

X → C + Y r2= k2xX

Y + B → D r3= k3xBxY

X

A + 2B → C + D

The deduction of a rate law (and the reaction order) from a mechanism can be done by employing different methods, and all include some kind of approximation of, or assumption concerning, the dynamics of the mechanism. The different methods can sometimes give the same rate law for a given mechanism. One method, called the steady state approximation, is employed in Section 2.5.3.

2.5.2 Equilibrium and Steady State

The simple reaction given by formula (2.1) is an irreversible reaction, meaning that it only proceeds in one direction. This is a simplification in the modelling of the dynamics; in fact, all reactions also have a reverse course of events compared to the forward reaction. For the dinitrogen pentoxide reaction this means that some amount of the product always decomposes back into reactants.

Reactions modelled as irreversible have a reverse reaction rate small enough to be neglected in modelling of the process. Reactions not possessing this property have forward and reverse reaction rates of the same magnitude; the reaction is reversible. This is depicted as A + B k1 k−1 C + D

If the reverse and forward reaction rates are equal, the reaction has reached chem-ical equilibrium. A chemchem-ical equilibrium is characterised by a minimum of the free energy for the reaction (Atkins and Jones, 1999). There is no inclination for change in any of the directions for the reaction.

The equilibrium constant is defined as

K =xCxD xAxB

If the reverse and forward reactions both are of simple second order, the equilibrium will correspond to k1xAxB = k−1xCxD, and K = k1/k−1. If K is large a lot of

(35)

product is produced before equilibrium is reached, since the rate constant for the forward reaction then is larger than the reverse rate constant (Atkins and Jones, 1999).

In cells, compounds can exist in concentrations far from their chemical equilibrium states (Zubay et al., 1995). These states are connected to a larger free energy, and the concentrations are not only governed by the external environment, as it is in chemical equilibria.

The rates in reaction sequences vary according to the cells requirements. The concentrations of key metabolites are held at constant levels by balancing the rates of production and consumption of reaction intermediates (Zubay et al., 1995).

A reaction with an intermediate species, X, is depicted below. The concentration of the intermediate is held at a constant level. Mathematically it corresponds to a derivative of zero, for the amount of the component; dxX/dt = 0. This is a non

chemical equilibrium situation where the intermediate resides in what is called a steady state.

A + B X C + D

2.5.3 Enzyme Kinetics

Enzymes are proteins acting as catalysts to chemical reactions in cells. Enzymes have an active site where the reaction takes place. The molecule which the en-zyme acts upon is called substrate. Enen-zymes work by decreasing the activation free energy, i.e., the energy needed for the substrate(s) to enter into a state of transformation, called the transition state (Zubay et al., 1995).

Enzyme catalysed reactions have a feature that distinguishes them from simpler, chemical reactions; they show saturation (Cornish-Bowden and Wharton, 1988). Almost all reactions of this type are of first order for small substrate concentrations, but the rate decreases with increased substrate concentration. The rate eventually becomes constant, independent of the concentration of substrate.

The behaviour can be observed if the rate of a reaction can be deduced from a set of experiments with different substrate concentrations. For each substrate concentration, a small amount of the catalysing enzyme is added, and the amount of formed product is monitored during a time span. The monitoring can be done with light spectrometry, provided that the product absorbs light. The amount of product is plotted against time, and the slope of the curve at the start of the experiment corresponds to the rate. An artificial curve of the typical behaviour is given in Figure 2.5.

The constants Kmand V in the figure are parameters in the rate law for the reaction

(36)

20 Theory [S] 0 0.5V V Km

Figure 2.5. Typical dependence of reaction rate on substrate concentration for a reaction

following Michaelis-Menten dynamics (see below). Km is the substrate concentration

giving a rate of 0.5V .

enzyme, S the substrate, and P the product.

E + S k1 k−1 ES k2 → E + P Etot− xES xS xES xP

The amount of free substrate is considered to be much larger than the amount bound to the enzyme-substrate complex, ES, and it is hence assumed that the free amount of substrate and the total amount are equal (Cornish-Bowden and Wharton, 1988). The conversion of ES to free product and free enzyme, E + P , is considered to be irreversible, if one only measures the initial rate of the reaction in the steady state (the regeneration of ES is negligible since the amount of product is small) (Zubay et al., 1995).

The steady state assumption for this mechanism is that the amount of the inter-mediate species does not change:

dxES

dt = k1(Etot− xES)xS− k−1xES− k2xES= 0 Solving the algebraic expression for xES gives that

xES=

k1EtotxS

k−1+ k2+ k1xS

and since the dissociation of the enzyme-substrate complex is of first order with respect to the intermediate, we have a rate of the form where Km= (k−1+ k2)/k1

and V = k2Etot: v = k1k2EtotxS k−1+ k2+ k1xS = V xS Km+ xS (2.5)

(37)

This equation is known as the Michaelis-Menten equation and Km consequently

as the Michaelis constant, although it was Briggs and Haldane that in 1925 pro-posed the mechanism (Cornish-Bowden and Wharton, 1988, Zubay et al., 1995). Michaelis and Menten in fact assumed that the first step was an equilibrium, which is a less general assumption than the one made above. The equation can be ex-tended to describing a mechanism with several substrate molecules.

Substances called activators and inhibitors affect enzyme activity and cause the catalysed reaction to proceed faster or slower respectively.

The most important kinds of inhibition, although several other exist, are results of the inhibitor binding to the enzyme (Cornish-Bowden and Wharton, 1988). The simplest type of inhibition is linear, meaning that terms proportional to the inhibitor concentration appear in the denominator of the rate law (Cornish-Bowden and Wharton, 1988). If the enzyme-catalysed reaction followed Michaelis-Menten dynamics in the absence of inhibitor it will still do so if inhibitor is present, with modifications of the effective parameter values.

Competitive inhibition is the most common kind (Cornish-Bowden and Wharton, 1988) and occurs when the substrate, S, and inhibitor, I, compete for the free enzyme. If the inhibitor binds, it will result in an inactive complex, EI, that does not lead to products. The Michaelis constant will have an altered effective value, while the limiting rate, V , will have the same effective value as before. The interpretation of this is that the enzyme-substrate complex is as reactive as before, but the effective affinity of the substrate to the enzyme is decreased (Cornish-Bowden and Wharton, 1988).

Uncompetitive inhibition, on the other hand, is a result of the inhibitor binding to the enzyme-substrate complex and producing an inactive complex, ESI. A pure form of this kind of inhibition is uncommon, and is most important as product inhibition (Cornish-Bowden and Wharton, 1988). The effective limiting rate and Michaelis constant are both affected, but their ratio is not.

If an inhibitor binds to both the free enzyme and to the enzyme-substrate com-plex, the inhibition is called mixed. In some books, the term non-competitive inhibition occurs, and is regarded as an additional form of inhibition. In this case, the inhibitor affects the enzyme without affecting binding of the substrate, and the inhibition is independent of substrate concentration (Cornish-Bowden and Whar-ton, 1988, Zubay et al., 1995). According to (Cornish-Bowden and WharWhar-ton, 1988), this is not a plausible mechanism in nature, and there are no examples recorded, if effects of the pH are excluded.

Activators are effectors that bind to an enzyme and increase enzyme activity with-out being changed in the reaction (Cornish-Bowden, 1995). Specific activation occurs when the enzyme, in the absence of the activator, does not have any activ-ity. In analogy with the competitive inhibition, the effective value of the Michaelis

(38)

22 Theory

constant is altered. An activation counterpart to the mixed inhibition also exists.

More complex expressions for the rate law of enzyme activation occur when the enzyme is less active, but not inactive, in the absence of activator (Cornish-Bowden, 1995).

A simple model for activation is to insert the concentration of the activator in the numerator of the rate law. The model might be used when the underlying mechanism of activation is unknown.

Most enzymes follow Michaelis-Menten kinetics, but some do not. These enzymes show a sigmoidal dependence on substrate of the rate, instead of a hyperbolic dependence as in Figure 2.5. Often, these enzymes have controlling, or regulatory, tasks in a biological network.

The sigmoidal response is linked to cooperativity of the enzyme. An enzyme ex-hibiting cooperativity will be ultra-sensitive to changes in the substrate concentra-tion, which is not the case of Michaelis-Menten kinetics (Cornish-Bowden, 1995). Positive cooperativity can for example occur if an enzyme is binding several sub-strate molecules and the binding of subsequent subsub-strates is facilitated by previous binding of substrate.

(39)

Chapter 3

Methods

Several identification methods were used in the thesis work and they are described and analysed in this chapter. Emphasis is given to the mathematical basis for each method and to the data required to apply them to biochemical systems.

3.1 Methods Based on Local Parameter

Pertur-bation Data

Local parameter perturbation data is based on steady state measurements. The amount of each interesting species does not change, which corresponds to a deriva-tive equal to zero for the variables representing these amounts.

Consider a network consisting of n components (or n groups of components), each described by a state variable, xi(t), that (in some units) represents the amount of

the component. The variables are gathered in a vector, x(t) = (x1(t), . . . , xn(t)).

The system is modelled by a set of differential equations where the rates of change of the variables are dependent, in addition to the variables xi, on a set of parameters,

p = (p1, . . . , pm). ˙ x1(t) = f1(x1(t), . . . , xn(t), p1, . . . , pm) ˙ x2(t) = f2(x1(t), . . . , xn(t), p1, . . . , pm) .. . ˙ xn(t) = fn(x1(t), . . . , xn(t), p1, . . . , pm) (3.1)

A component fi is not explicitly dependent on all the parameters in p, but the

exact dependence is not known. The steady state is given by ˙x(t) = f (ξ(¯p)) = 0,

(40)

24 Methods

where ξ(¯p) is a vector of the steady state values, associated to a specific set of parameter values ¯p, of the state variables.

We assume that it is possible to perturb the parameters in p individually. A realisation of this could be the addition of specific inhibitors or activators that affect the catalysing enzymes in the network. Another possibility is to introduce a double stranded RNA that interferes with a gene product, which can subsequently lead to a reduction of the amount of a protein. Interference RNA is often denoted RNAi. Each affected enzyme is responsible for catalysing a reaction dependent upon a parameter pj. Individual additions of activators, RNAi, or inhibitors are

repeated for all parameters.

The result will be m + 1 different steady state measurements, of all state variables, including the reference steady state, ξ_ref. All of the different measurements are associated to a specific set of parameters.

A perturbation experiment involving an alteration of the parameter pj results in

the following quotient, which approximates the sensitivities of the steady state values to changes in the parameter:

σij = ∆ξi ∆pj = ξi− ξref pj− p_{j ref} , i = 1, . . . , n (3.2)

The calculations are repeated for all parameters and the result is gathered in a ma-trix Σ = (σij). This matrix has the dimensions (n × m) and it is the approximation

of ∂ξi/∂pj for i = 1, . . . , n and j = 1, . . . , m.

If the exact changes of the parameter values are unknown, the sensitivities can be approximated in a different manner, explained in (Kholodenko and Sontag, 2002).

3.1.1 Interaction Graph Determination

An account of the top-down identification method from (Kholodenko and Sontag, 2002) is given in this section. Assuming that the components in the network are known, the method aims at determining the interaction graph of the system.

The interaction graph, also named connection graph, is a description of how each node in the network qualitatively affects the other nodes. The interaction graph is a basis for further investigation of the network and its dynamic behaviour. A simple interaction graph for a network of three components is shown in Figure 3.1.

How the different nodes in the network affect each other is displayed with the arrows - the net effect is negative or positive with respect to the amount of the component(s) represented by each node. Exactly how the network is built is not always trivial to extract from the interaction graph. The graph only describes the functional dependencies of the total rate of change of each variable with respect to

(41)

3.1 Methods Based on Local Parameter Perturbation Data 25

Figure 3.1. A simple interaction graph.

other variables in the system.

The deduction of the network wiring is almost always impossible if the nodes of the network are connected by mass flow, since several network structures can cor-respond to the same graph. The impossibility of distinguishing between network structures is shown in the example below.

The method also allows a modular approach, where several components can be represented by each node. The same method, with minor variations, is described in (Kholodenko et al., 2002) and is applied to a gene regulatory network.

The interaction graph can be deduced from the Jacobian, f_x= (fik), of the general

system (3.1) of differential equations, where

fik=

∂fi

∂xk

(¯x, ¯p)

A connection from node i to node k in the interaction graph, corresponds to a non-zero element in position (k, i) of the Jacobian.

Example 1 Network uniqueness

Two network architectures are given in Figure 3.2. The two networks have some mass flow between their nodes.

Both systems have the following quantitative Jacobian. The deduction of the net-work architecture from this Jacobian will not lead to a unique netnet-work.

    A B C A − 0 + B + − − C 0 0 −    

(42)

26 Methods Activate Inhibit A B C A B C

Figure 3.2. Two network architectures.

The intention of finding the Jacobian of the system from steady state data cannot be entirely realised, since a multiplication of each row in the Jacobian with a constant also is a solution to ˙x = 0. Hence, the rows can only be determined up to a scalar multiple. This is still enough to determine the interaction graph since we are only interested in whether an element of the Jacobian is non-zero or not.

A few assumptions are needed to make the method legitimate. There must, for each node, represented by a variable xi, exist other nodes that are not directly

connected to the current node i. This means that there is a set of parameters for which ∂fi/∂pj= 0 for each i ∈ {1, . . . , n}. The required set of parameters have to,

for each xi, correspond to n − 1 independent columns from the Σ-matrix.

For each node i, perturbations must be made to n − 1 other nodes, and the pa-rameters that are perturbed cannot be a part of any of the rates that lead to or from the node i. It is not necessary to perform n · (n − 1) perturbations, since each perturbation can be used for more than one node. To use the method, some information about how the nodes are connected must be known before hand. The information needed is, for example, that a subset of the nodes are known to be unconnected to the remaining set of nodes. If the network connections only are of controlling nature, the assumptions are easy to fulfill.

For one state, i.e., one row in the Jacobian, the following is valid: ∂ ∂pj fi(¯x, ¯p) = ∂fi ∂x(¯x, ¯p) ∂x ∂pj (¯x) + ∂fi ∂pj (¯x, ¯p)

The assumptions above and the fact that the system is in steady state produces the pursuing orthogonality criterion:

∂fi

∂x(¯x, ¯p) ∂x ∂pj

(¯x) = 0 (3.3)

for each pj that fulfills the first assumption for the particular fi.

This information is all we need to reproduce the rows of the Jacobian up to scalar multiples. The assumptions guarantee that the orthogonality criteria correspond

(43)

3.1 Methods Based on Local Parameter Perturbation Data 27

to a system of linear equations from which we can produce the values of each row. Also, by fixing the diagonal elements in the Jacobian to −1, the linear equation system is reduced and can be solved using least square methods.

3.1.2 Determination of Control Loops in Mass Flow

Net-works

A problem with the interaction graph identification method in Section 3.1.1 is that the interaction graph can correspond to several non-separable network architec-tures. It is not possible to separate mass flow connections from control loops, and the interaction graph tends to become complex when mass flow is present.

The method of this section aims to identify the control loops in a network where the mass flow is known before hand. The applicability of the method is of course limited by the assumption that the mass flow is known, but it is not an at all unlikely situation. Situations may exist where the ”back-bone” of the network is known from experiments, but nothing of the internal control mechanisms.

The set of differential equations from the system (3.1) can be represented in the following manner where the fi components are replaced by the rates. The time

dependence has been left out, as well as the dependence of the rates on the state variables, to simplify the notations:

˙ x1 = m1,1r1(p1) + m1,2r2(p2) + . . . + m1,qrq(pq) ˙ x2 = m2,1r1(p1) + m2,2r2(p2) + . . . + m2,qrq(pq) .. . ˙ xn = mn,1r1(p1) + mn,2r2(p2) + . . . + mn,qrq(pq)

The kinetics of the rates are not known but the mass flow in the network is reflected by the knowledge of the mi,l elements. If the stoichiometric ratios of the network

is one-to-one for all species, the matrix M = (mi,l) only contains 0, 1 or −1 in

appropriate places. Row i of the M -matrix has zeros for the rates not leading to or from the node represented by variable i. The differential equations are summarised as

˙

x = M r(p) (3.4)

Examining a specific row i, denoted M (i)r(p), in the relation above and differen-tiating it with respect to the parameter pj will result in:

∂ ∂pj (M (i)r(p)) = ∂ ∂pj (mi,1r1+ . . . + mi,qrq) = mi,1∇r1 ∂x ∂pj + mi,1 ∂r1 ∂pj + . . . + mi,q∇rq ∂x ∂pj + mi,q ∂rq ∂pj (3.5)

(44)

28 Methods

The differentiation of the arbitrary row from expression (3.5) can be represented in a more compact form, where rx = (∇r1∇r2. . . ∇rn)T quantifies how the rates

depend on the state variables:

M (i)rx ∂x ∂pj + M (i)(∂r1 ∂pj , . . . ,∂rq ∂pj )T (3.6)

The parameter pj in equation (3.6) can be chosen to produce a cancellation of the

last term of the equation. If the rates that lead to or from state that is represented by the row i are independent of the parameter, we have cancellation. Selecting a parameter with these properties will produce the following criterion, similar to the orthogonality criterion of equation (3.3):

M (i)rx

∂x ∂pj

= 0 (3.7)

The vector ∂x/∂pj is approximated by the corresponding column of the Σ-matrix

defined earlier in this chapter. Since the matrix M is known, the criterion results in an equation of some of the entries in the rx-matrix. Combining all different

rows and ∂x/∂pj vectors as in equation (3.7) will produce an under determined

equation system, assuming that each rate cannot be perturbed more than once. The assumption is reasonable, since we only perform small perturbations and the response should be similar whatever kind of perturbation we use on parameters in the same rate, i.e., additional perturbations will not add any independent equations. Further analysis of the method can be found in Chapter 4.

3.2 Methods Employing Time Series Data

The methods in the previous section are based on steady state data. As noted, steady state data is not as informative as time series data and it is not possible to exactly determine the Jacobian of the system.

Time series data is retrieved by recording the amounts of the involved species during a time span. The data is informative but hard to collect. Measurements in vitro of dynamic species is not an easy task, but can be achieved using elaborate methods including different kinds of spectrometry.

Sampling of the data is needed, for both continuous and discrete data records. The intention is for the data to pick up the dynamics of the system and in order to do so, a lot of thought must underlie the choice of sampling interval.

The sampling interval, T , is the time between samples. The angular frequency of the sampling is subsequently defined as ωs = 2π/T = 2πfs and naturally called

(45)

3.2 Methods Employing Time Series Data 29

A continuous signal is presumed to have a Fourier transform describing its frequency content. Sufficient, but not necessary, conditions for the transform to exist are the Dirichlet conditions (Sv¨ardstr¨om, 1999). Periodic functions as well as the unit step function do not fulfill the requirements in these conditions. For these signals, a transform is defined in the limit allowing generalised functions like the Dirac pulse. The transform is defined as

F [x(t)] = X(ω) = Z ∞

−∞

x(t)e−jωtdt

where ω is the angular frequency.

Sampling of a continuous signal is modelled by a multiplication of the signal with a sequence of impulses, δT(t) =P∞_n=−∞δ(t−nT ), where δ represents the Dirac pulse.

The sampling interval is the time between the pulses. The procedure of sampling is illustrated in Figure 3.3. The transform of the signal is repeated with a distance of ωsrad/s in the frequency domain. The signal is limited in the time domain, since

we only can sample for a finite time span. A time limited signal cannot be band limited; the transform is infinite in the frequency domain (Sv¨ardstr¨om, 1999).

t t t

∑

− ⋅ ( ) ) (t t nT x δ

∑

δ(t−nT) ) (t x T 0 0 T 0

Figure 3.3. Sampling of a continuous signal.

The unlimited transform causes problems during sampling. The frequencies above the Nyquist frequency, ωN = πfs = ωs/2 rad/s, are folded into the transform,

and are instead misinterpreted as lower frequencies. The effect is illustrated in Figure 3.4. Consequently, there are two problems in sampling; folding of frequencies above the Nyquist frequency and the necessity of limiting the signal in the process. Both effects distort the transform of the continuous signal.

(46)

30 Methods X( ) s N Xs( )

Figure 3.4. Frequencies above ωNare folded into the angular frequency range [−ωN, ωN].

The folding effect is also known as aliasing. A way to avoid the aliasing-effect is to use an anti-alias filter on the data before the sampling is made. An anti-alias filter is a regular low-pass filter with the cut-off frequency chosen a bit below the Nyquist frequency (Ljung and Glad, 2004). In this way, the frequency content above ωN is

lost, and does not disturb the remaining information by folding effects. In addition, if there is high-frequency noise in the data, it is also removed by the anti-alias filter (Ljung, 1987).

A rule of thumb is to choose the sampling frequency a ten-fold greater than the bandwidth we wish our model of the system to cover. The choice corresponds to 4-8 sample points on the flank of the step response from the system (Ljung and Glad, 2004). It might therefore be of interest to measure a step response from the system before hand.

3.2.1 Linearising Around an Operating Point

Consider a non-linear network that resides in a steady state, which is the operating point. How the network is wired in this steady state is, as before, given by the Jacobian of the system at the operating point and is what we wish to determine.

The system is described using the following setup, similar to the one given by equation (3.1): ˙ x1(t) = f1(x1(t), . . . , xn(t), u1(t), . . . , uk(t)) ˙ x2(t) = f2(x1(t), . . . , xn(t), u1(t), . . . , uk(t)) .. . ˙ xn(t) = fn(x1(t), . . . , xn(t), u1(t), . . . , uk(t)) (3.8)

Compared to equation (3.1), the dependence on the parameters is not explicitly stated in the representation and the variables ui(t) are added that correspond to

(47)

input signals to the system. The variables are gathered in an input signal vector, u(t). The representation in equation (3.8) implies that there is a way to effect the system, in a controlled manner, by changing the input signals.

The steady state of the system corresponds to the specific values (x0, u0) of the

state variables and the input signals. The system is observed in a vicinity of the steady states where small deviations, (x(t) − x0, u(t) − u0), are represented

by (_ex(t),u(t)). The deviations are often called perturbations. Substituting the_e expressions of the deviations into the differential equation system of (3.8) and omitting the time dependence will give

d

dt(x0+x) = f (xe 0+x, ue 0+u)e (3.9) By expanding the left hand side of the expression and noting that dx0/dt = 0 and

du0/dt = 0 as well as expanding the right hand side in a Taylor series around the

steady state (x0, u0) we get

d dt(x) = f (xe 0, u0) + ∂f ∂x(x0, u0)x +e ∂f ∂u(x0, u0)u + h.o.t.e

The perturbations around the steady state are small and the higher order terms can be neglected. We have a linear system of differential equations for the perturbations e

x, since f (x0, u0) = 0 by definition of a steady state:

˙ e x = ∂f ∂x(x0, u0)x +e ∂f ∂u(x0, u0)eu

The matrix ∂f /∂x is the Jacobian for the system in the steady state and describes the network wiring, as stated before.

In equation (3.8), the variables ui(t) represent input signals. The approximate

input-output modelling of the system is completed by defining a set of output signals, y(t), as linear combinations of the state variables as well as the input signals. The result is a linear state space model:

˙ e

x(t) = A(θ)x(t) + B(θ)_e u(t)_e

y(t) = C(θ)x(t) + D(θ)_e u(t)_e (3.10) The time dependence is explicitly stated to point out that θ is time independent. θ is a vector of the parameters, determining the behaviour of the system, that were omitted in equation (3.8). The true set of the parameters is what we wish to find.

Noise is inevitably added to the measurements during data collection. The state space model given by equation (3.10) does not contain any noise model. Adding a noise model is advisory, and a fairly simple model is to include additive noise sources. An extended model for the system, including noise sources w(t) and v(t) representing process noise and measurement noise respectively, is (Ljung, 1987):

˙ e

x(t) = A(θ)x(t) + B(θ)_e u(t) + w(t)_e

(48)

32 Methods

The noise sources are modelled as sequences of independent random variables with zero mean. A so called observer estimating the state variables is given by equa-tion (3.12). To avoid cumbersome notaequa-tions, the˜-accent is dropped, meaning that ˆ

x in fact is ˆx. The parameter dependence is also omitted for the same reason._e ˙ˆx(t) = Aˆx(t) + Bu(t) + K(y(t) − C ˆx(t) − Du(t)) (3.12) The quantity y(t) − C ˆx(t) − Du(t) is a measurement of how well the predictor ˆ

x(t) estimates x(t). The observer that minimises the estimation error is called the Kalman filter, represented by K. The filter is calculated from the matrices for the state space description as well as the variance and covariance matrices for the noise sources. For a more thorough treatment, see for example (Ljung and Glad, 1997) or (Ljung, 1987).

The prediction error ν(t) = y(t) − C ˆx(t) − Du(t), used as feedback quantity in the observer, is the new information in the measurement y(t) that is not available in the previous measurements. The signal ν(t) is called the innovation (Ljung and Glad, 1997, Ljung, 1987). The innovations form of the state space description of equations 3.11 is given by

˙ ˆ

x(t) = A(θ)ˆx(t) + B(θ)u(t) + K(θ)ν(t)

y(t) = C(θ)ˆx(t) + D(θ)u(t) + ν(t) (3.13)

A set of transfer functions can be deduced for the linear differential equation system on the innovations form. The mathematical expressions for the transfer functions as well as a block diagram are given in equation (3.14) and Figure 3.5, respectively. The dependence on the parameter vector θ is omitted.

Y (s) = (C(sI − A)−1B + D) | {z } G(s) U (s) + ((sI − A)−1K + I) | {z } H(s) V (s) (3.14) Σ+ u G(s) y H(s) +

Figure 3.5. A block representation of the system.

If an identification process on the system given by equation (3.13) is to be successful, some conditions need to be fulfilled.

(49)

Identifiability is central to the identification process and deals with the problem of the process being able to yield true values of the parameters of the model and/or giving a model equal to the true system (Ljung, 1987).

Assuming that our sampled data from the system is informative enough to represent the system behaviour, the question is if two parameterisations of the system can produce equal models, i.e., give non-separable output data. Another way to put the question is if the model structure is invertible or not (Ljung, 1987). Invertibility is a necessary condition for identifiability and a model,M, is identifiable with a true parametrisation θ∗if

M(θ) = M(θ∗_{) ⇒ θ = θ}∗_, _{θ ∈ D} M

where θ and θ∗ represent different parameterisations in the set of possible values, D_M, for the model (Ljung and Glad, 2004, Ljung, 1987). The condition above guarantees global identifiability. Local identifiability is defined in a similar way, but the invertibility is only valid in a vicinity of θ∗.

The input signals, u(t), must be chosen to excite the system (Ljung and Glad, 2004), giving data that is informative enough for the identification, as stated above. The input signals must have variations as fast as the smallest time constants of the system. For a sampling process to pick up the frequencies of the oscillations, it is necessary for the input signals to contain these frequencies. Equally important is then to sample fast enough not to lose the frequencies by folding or anti-alias filtering.

An input signal alternating between two levels is a good choice for a linear system, like the one given by equation (3.10). If the variations are random, this signal contains all frequencies (Ljung and Glad, 2004). Since we have a linearised system, it is important to verify that the variations in the input signal are small enough for the system to respond in a linear fashion. Crucial to the validity of the linearisation is closeness of the system to the steady state.

A common signal in identification of linear systems is a Pseudo Random Binary Signal (PRBS)(Ljung and Glad, 2004). A PRBS is a signal that alternates between two levels in an apparent stochastic manner, although it is deterministic.

Our question is now if the state space model in equation (3.13) is identifiable, provided that the inputs excite the system in a satisfactory way. In our application, the matrix D(θ) is zero, meaning that the inputs do not directly affect the outputs. By the theorem (4A.2) in (Ljung, 1987) it is stated that the system is globally and locally identifiable, given that it is parameterised according to the observer canonical form, at θ∗ if and only if {A(θ∗), [B(θ∗)K(θ∗)]} is controllable.

The observer canonical form for a multiple input - multiple output (MIMO) system is characterised by a full parametrisation of the B(θ)-matrix while the C(θ)-matrix is built up by zeros and ones. The A(θ)-matrix has the same number of fully

(50)

34 Methods

parameterised rows as the number of output channels, while the remaining rows are filled with zeros and ones in a certain pattern (Ljung, 2001).

The exact definition of the observer canonical form can be found in (Ljung, 1987) and the following is an example with three state variables, two output variables and one input variable (a SIMO system):

A(θ) =   0 1 0 × × × × × ×  , B(θ) =   × × ×   K(θ) =   × × × × × ×  , C(θ) = 1 0 0 0 0 1

To achieve a unique connection between the input and output signals, we have to limit the flexibility in the A(θ)-matrix by fixating certain elements, since we do not measure all the state variables.

The observer canonical form for the state space model is a sufficient parametrisation for some of the systems we intend to apply the method to, as stated in Chapter 4.

We have now stated that the method of this section can be used to identify lin-earised systems around an operating point. The identification itself can, for exam-ple, be made with the Systems Identification Toolbox for Matlab (SITB), which is designed for identification of linear systems.

3.2.2 A Discrete Linearisation Method

A less general and discrete method based on the same kind of linearisation as the method in the previous section is briefly explained here. The method illustrates a simple version of the estimation procedure of the SITB, and is included for basic understanding of parameter estimation.

Time series data is utilised and the availability of measurements of all state variables is a requirement, i.e., Cd(θ) is the identity matrix. The matrix Dd(θ) is assumed

to be zero. The subscript d emphasises that we have the discrete equivalents of the matrices from the continuous system of equation (3.10).

A difference equation is the basis for a discrete modelling of the system. The system is dependent on a set of parameters, but the dependence is omitted in the description:

xn+1= Adxn+ Bdun

y_n= Cdxn

(3.15)

An estimation of the matrix Ad can be used to deduce the Jacobian for the

(51)

following matrices are needed:

R = xn+1 xn . . . x1 , M =

xn xn−1 . . . x0

un un−1 . . . u0

By inserting the matrices into the system description (3.15) we will get

R =

Ad Bd M

To estimate Ad we multiply the expression from the right with the pseudoinverse

of M . The resulting matrix from this operation is of dimensions (n × n + k) where n is the number of state variables and k the number of input signals. The first n columns correspond to Ad. Taking the matrix logarithm of the relation Ad= eJ ∆T

and dividing by the sampling interval produces an estimated Jacobian.

The estimation of the SITB is more complex, and includes noise models, which is not the case here. The absence of noise models renders the method very sensitive to noise.

(52)

Evaluation and Development of Methods for Identification of Biochemical Networks

Department of Physics and Measurement Technology

Master’s Thesis

Evaluation and Development of Methods for

Identification of Biochemical Networks

Alexandra Jauhiainen

Evaluation and Development of Methods for

Identification of Biochemical Networks

Alexandra Jauhiainen

Abstract

Acknowledgement

Notation

Symbols

Abbreviations

Contents

List of Figures

List of Tables

Chapter 1

Introduction

1.1

Problem Formulation

1.2

Systems Biology

1.2.1

What is Systems Biology?

1.2.2

Why Perform Research on Systems Biology?

1.3

The Thesis

1.3.1

Aim and Audience

1.3.2

Structure

Chapter 2

Theory

2.1

Systems, Modelling, and Simulation

2.1.1

Systems

2.1.2

Models

2.1.3

Simulation

2.2

System Identification

2.2.1

The Identification Process

2.2.2

Validation

2.3

Biological Network Structures

2.3.1

Metabolic Networks

2.3.2

Gene Regulatory Networks

2.3.3

Signalling Networks

2.4

Identification of Biological Networks

2.4.1

Identification Approaches

2.5

Chemical Kinetics

2.5.1

Rate Laws and Reaction Mechanisms

2.5.2

Equilibrium and Steady State

2.5.3

Enzyme Kinetics

Chapter 3

Methods

3.1

Methods Based on Local Parameter

Pertur-bation Data

3.1.1

Interaction Graph Determination

3.1.2

Determination of Control Loops in Mass Flow

Net-works

3.2