Modeling of metabolic insulin signaling in adipocytes

Full text

(1)

Institutionen för systemteknik

Department of Electrical Engineering

Examensarbete

Modeling of metabolic insulin signaling in

adipocytes

Examensarbete utfört i Reglerteknik vid Tekniska högskolan i Linköping

av

Erik Ulfhielm

LITH-ISY-EX--06/3778--SE

Linköping 2006

Department of Electrical Engineering Linköpings tekniska högskola

Linköpings universitet Linköpings universitet

(2)
(3)

Modeling of metabolic insulin signaling in

adipocytes

Examensarbete utfört i Reglerteknik

vid Tekniska högskolan i Linköping

av

Erik Ulfhielm

LITH-ISY-EX--06/3778--SE

Handledare: Henrik Tidefelt

isy, Linköpings universitet

Gunnar Cedersund

fcc, Fraunhofer Chalmers centre

Peter Strålfors

ibk, Linköpings universitet

Examinator: Jacob Roll

isy, Linköpings universitet

(4)
(5)

Avdelning, Institution

Division, Department

Division of Automatic Control Department of Electrical Engineering Linköpings universitet S-581 83 Linköping, Sweden Datum Date 2006-02-05 Språk Language  Svenska/Swedish  Engelska/English  ⊠ Rapporttyp Report category  Licentiatavhandling  Examensarbete  C-uppsats  D-uppsats  Övrig rapport  ⊠

URL för elektronisk version

http://www.control.isy.liu.se http://www.ep.liu.se/2006/3778 ISBNISRN LITH-ISY-EX--06/3778--SE

Serietitel och serienummer

Title of series, numbering ISSN

Titel

Title Modellering av metabol insulinsignalering i fettcellerModeling of metabolic insulin signaling in adipocytes

Författare

Author Erik Ulfhielm

Sammanfattning

Abstract

Active insulin receptors (IR) phosphorylate insulin receptor substrate (IRS), but it is not clear whether IRS is phosphorylated mainly by IR at the plasma membrane or by internalized IR in the cytosol. In this thesis, structural identifiability analysis and parameter sensitivity analysis is performed for models of the first steps in the metabolic insulin signaling pathway. In particular, the identifiability of the kinetic parameters governing IRS phosphorylation are investigated.

Given measurements of the relative increase in phosphorylation degree of IR and IRS, the structural identifiability analysis revealed that the parameters gov-erning IRS phosphorylation are non-identifiable, but their ratio is identifiable. This is sufficient to study whether phosphorylation of IRS proceeds more rapidly by IR at the plasma membrane or by internalized IR in the cytosol. In the exam-ined model structure, internalization of insulin receptors is shown to be necessary to reproduce the experimental data.

Sensitivity analysis of the parameters governing IRS phosphorylation showed that many parameters need to be known in order to obtain “practical identifiabil-ity” of the interesting parameters.

Nyckelord

(6)
(7)

Abstract

Active insulin receptors (IR) phosphorylate insulin receptor substrate (IRS), but it is not clear whether IRS is phosphorylated mainly by IR at the plasma membrane or by internalized IR in the cytosol. In this thesis, structural identifiability analysis and parameter sensitivity analysis is performed for models of the first steps in the metabolic insulin signaling pathway. In particular, the identifiability of the kinetic parameters governing IRS phosphorylation are investigated.

Given measurements of the relative increase in phosphorylation degree of IR and IRS, the structural identifiability analysis revealed that the parameters governing IRS phosphorylation are non-identifiable, but their ratio is identifiable. This is sufficient to study whether phosphorylation of IRS proceeds more rapidly by IR at the plasma membrane or by internalized IR in the cytosol. In the examined model structure, internalization of insulin receptors is shown to be necessary to reproduce the experimental data.

Sensitivity analysis of the parameters governing IRS phosphorylation showed that many parameters need to be known in order to obtain “practical identifiability” of the interesting parameters.

(8)
(9)

Acknowledgements

I would like to express my gratitude to Jacob Roll, Henrik Tidefelt and Gunnar Cedersund for their supervision. I would also like to thank Siri Fagerholm and Peter Strålfors for rewarding discussions about experimental and biological aspects of insulin signaling. Many thanks to Jenny Falk for doing the opposition. Finally, i would like to thank my family and Sofia Åberg for always supporting me.

(10)
(11)

Contents

1 Introduction 1 1.1 Problem formulation . . . 1 1.2 Objectives . . . 1 1.3 Thesis disposition . . . 2 1.4 Abbreviations . . . 2

2 Background and theory 5 2.1 Systems biology . . . 5

2.2 Basic cellular biology . . . 5

2.3 The biology of insulin . . . 6

2.3.1 Physiological actions of insulin . . . 7

2.3.2 Insulin signaling in adipocytes . . . 7

2.4 Systems, models and simulations . . . 9

2.4.1 Systems . . . 9

2.4.2 Mathematical models . . . 9

2.4.3 Simulations . . . 10

2.5 A mathematical model of insulin signaling . . . 11

2.5.1 Insulin receptor subsystem . . . 11

2.5.2 Postreceptor signaling . . . 13 2.6 Dimensional analysis . . . 13 2.7 Observability . . . 16 2.7.1 Parameter identifiability . . . 16 2.7.2 Linear systems . . . 16 2.7.3 Nonlinear systems . . . 17

2.8 Controllability and observability gramians . . . 18

2.8.1 Empirical gramians . . . 19

2.9 Parameter estimation . . . 20

2.10 Parameter sensitivity analysis . . . 21

2.10.1 Derivatives of the prediction . . . 21

3 Experimental data 23 3.1 Time series . . . 23

3.2 Excitation signal . . . 23

3.3 Output equations . . . 24 ix

(12)

3.4 Data sets . . . 24 3.4.1 Pre-processing of data . . . 24 4 Results 27 4.1 Model structures . . . 27 4.1.1 Empirical gramians . . . 31 4.1.2 Some observations . . . 33 4.1.3 Non-dimensionalization of M3 . . . 33 4.2 Structural observability . . . 34

4.2.1 Measurements that give observability . . . 37

4.3 Simulations . . . 37

4.3.1 Comparison with experimental data . . . 38

4.3.2 Transfer function of the IR-subsystem . . . 39

4.3.3 Recreating the overshoot of the phosphorylated IR . . . 41

4.4 Sensitivity analysis . . . 42

5 Conclusions 45 5.1 Summary of the thesis . . . 45

5.2 Discussion . . . 46

5.2.1 Model structures . . . 46

5.2.2 Structural identifiability . . . 47

5.2.3 Parameter sensitivity . . . 47

Bibliography 49

(13)

Chapter 1

Introduction

This chapter introduces the reader to the thesis. The raison d’être for doing this work is explained, the objectives of the thesis are stated, and the thesis disposition is outlined.

For a short introduction to basic cellular biology and the biology of insulin, see Sections 2.2 and 2.3.

1.1

Problem formulation

Insulin receptors (IR) that are stimulated with insulin become phosphorylated and later internalized, i.e., the receptor leaves the plasma membrane and enters the cytosol in an endosome.

Activated IR has tyrosine kinase activity and phosphorylates various intracellular substrates, one of the major being the insulin receptor substrate (IRS)-family. There is today no consensus whether IRS is phosphorylated mainly by IR residing in the plasma membrane or by endocytotic IR. The fact that this is a question difficult to answer by direct measurements makes the application of mathematical modeling and system identification techniques interesting.

1.2

Objectives

The objectives of this thesis are:

• To modify an existing mathematical model of insulin signaling by including the possibility of IRS phosphorylation by internalized IR.

(14)

2 Introduction

• To investigate the possibility of doing parameter estimation on the model by carrying out observability and sensitivity analysis of the model parameters, in particular, the parameters governing IRS phosphorylation.

• To compare simulations of the model to time series data.

1.3

Thesis disposition

After this introduction, a chapter with the background and theory needed in this thesis is presented. A reader familiar with the biology of insulin can skip the corresponding section, and a reader accustomed to the notation and methods of control theory and system identification can omit those sections.

The third chapter contains a brief description of the experiments that generated the data, the resulting data sets and the preprocessing of them.

In the fourth chapter, a few model structures are presented, followed by the results of structural observability analysis, simulations and sensitivity analysis.

The fifth and last chapter contains a discussion of the results, along with some conclusions and directions for further work.

1.4

Abbreviations

Biological abbreviations

• DNA : Deoxyribonucleic acid • ECM : Extracellular matrix • GLUT4 : Glucose transporter 4

• IDDM : Insulin dependent diabetes mellitus • IR : Insulin receptor

• IRS : Insulin receptor substrate

• NIDDM : Non-insulin-dependent diabetes mellitus • PI3-K : Phosphatidylinositol 3-kinase

• PTP : Protein tyrosine phosphatase • RTK : Receptor tyrosine kinase

(15)

1.4 Abbreviations 3 Mathematical abbreviations

• DAE : Differential algebraic equation • LTI : Linear, time invariant

• NDF : Numerical differentiation formula • ODE : Ordinary differential equation • SEM : Standard error of the mean

(16)
(17)

Chapter 2

Background and theory

This chapter describes the background and theory needed in this thesis. A reader familiar with cellular biology and insulin metabolism can skip or skim the first sec-tions, and a reader familiar with the subjects of system identification and control theory can skip those sections.

2.1

Systems biology

Systems biology aims at systems level understanding of biology and biochemical networks. The structure and dynamics of a system are studied, rather than isolated parts of a cell or an organism, and system properties such as robustness and feedback become central issues [8].

A work of this type is interdisciplinary with points of contact in both molecular biology and control theory.

2.2

Basic cellular biology

All living organisms are composed of cells. A cell is a biomolecule containing com-partment, surrounded by an envelope called the plasma membrane. The plasma membrane consists mainly of phospholipids, cholesterol and structural proteins. A cell is internally organized into several smaller, membrane enclosed organelles, e.g. the nucleus and the mitochondria. The nucleus contains the genetic code, stored in DNA (deoxyribonucleic acid). The mitochondria is the “power plant” of the cell, where most of the energy-carrying molecule ATP (adenosine triphosphate) is produced. The cytosol is the cellular content within the plasma membrane, and outside the organelles.

(18)

6 Background and theory

The carrier of the genetic code, the DNA, consists of two complementary strings of nucleotides (also called bases), attached to a sugar-phosphate backbone. There are four different nucleotides present in DNA; adenine–A, cytosine–C, guanine–G and thymine-T. The bases are complementary, A always binds to T, and C always binds to G. The genetic information is stored in genes. A gene is a section of DNA that encodes a specific functional product (a protein).

The DNA-encoded genes are transformed into RNA (ribonucleic acid) in a process called transcription. RNA is further processed into mRNA, i.e., messenger RNA, which is transported out of the nucleus. The mRNA is then used as a template for

ribosomes, which build proteins by connecting amino acids to form polypeptide chains, a process called translation. The polypeptide chain folds into a complex

three-dimensional structure, where it tries to minimize its potential energy. The biological function of a protein is closely related to its 3-D structure.

Proteins play diverse roles in cells. Some are structural proteins that provide the basic morphology of the cell. Enzymes are proteins that catalyze chemical reactions. Many proteins reside in the plasma membrane, where they for instance can receive signals from the exterior of the cell, by acting as sensors for specific substances. Such proteins are called receptors.

The transformation of energy and matter within a cell or an organism is called

metabolism. The reactions in the cellular metabolism are governed by a large

num-ber of enzymes. Two important classes of enzymes are kinases and phosphatases. A kinase attaches a phosphate group from a donor, e.g. ATP, to an acceptor, e.g. a protein. The process of attaching a phosphate group to an acceptor is called

phosphorylation. This is an important reaction, since the phosphate group often

alters the shape – and hence the activity – of proteins. A phosphatase reverses the action of a kinase, by removing a phosphate group.

In a multicellular organism, communication between the cells is vital to maintain

homeostasis, i.e., a relatively constant internal environment. The communication

between cells can be transmitted by chemical substances from secretory cells, and received by receptors on the surface of target cells. If the signal is transmitted only to the neighbors of the secreting cell, it is called paracrine signaling. If the signal transmitter is transported to the receiver via the blood, it is called endocrine signaling.

In multicellular organisms, the cells are organized in tissues. The volume between the cells in a tissue is called the extracellular matrix (ECM). The ECM contains a complex mixture of carbohydrates and proteins.

2.3

The biology of insulin

This section describes the actions of insulin, both as a part of the endocrine (hor-monal) signaling system and how the signal is propagated at a cellular and molec-ular level in adipocytes (fat cells).

(19)

2.3 The biology of insulin 7

2.3.1

Physiological actions of insulin

Insulin is a peptide hormone which is a major regulator of glucose homeostasis. It is produced in pancreatic β-cells as a response to increased glucose level in the blood. From the pancreas it is then released into the blood, which transports it to its target tissues, which are mainly adipose (fat), hepatic (liver), and muscles. In the target tissues, insulin can activate both a metabolic and a mitogenic (cell division stimulating) response. The metabolic response stimulates the uptake and storage of nutrients from the blood and inhibits the release and degradation of triglycerides and glycogen. The result of the metabolic response are seen within minutes or hours, while the mitogenic response is more long-term and affects cell growth and proliferation.

The concentration of glucose in blood plasma is kept within narrow limits, ap-proximately between 4 and 10 mM [3]. To maintain glucose homeostasis is vital to the organism; disruptions can cause the serious disease diabetes mellitus. Dia-betes is divided in two types; type I, insulin-dependent diaDia-betes mellitus (IDDM), and type II, non-insulin-dependent diabetes mellitus (NIDDM). Diabetes type I commonly affects children, and diabetes type II affects adults, especially obese adults.

Diabetes mellitus type II is the most common form. It is the most common serious metabolic disease in the western world, and it is caused by insulin resistance in target tissues and failure to produce insulin in the pancreas. Diabetes demand external supply of insulin, normally through injections [1].

2.3.2

Insulin signaling in adipocytes

Adipose tissue is one of the main targets of insulin. It consists mainly of adipocytes held together in a framework of collagen (a structural protein abundant in the ECM). The adipocytes are highly differentiated, non-dividing large cells; they can reach 120 µm in diameter. Their main constituents are a large droplet of triglycerides and a nucleus (see Figure 2.1). The cytosol is only a very thin rim between the plasma membrane and the lipid droplet [1].

In adipocytes the IR is localized in structures called caveolae. The caveolae are about 50 nm in diameter, flask-shaped invaginations in the membrane, stabilized by the protein caveolin (see Figure 2.2).

The IR consists of two α- and two β-subunits, where the α-subunits reside on the extracellular side of the plasma membrane and contain ligand binding sites. The β-subunits are transmembrane, and connected through a disulfide bond to the α-subunits on the cell surface. The β-subunits have tyrosine kinase activity and can phosphorylate a number of substrates.

The insulin receptor belongs to the family of receptor tyrosine kinases (RTKs). Propagation of a signal by an RTK is initiated by ligand binding, which causes a

(20)

8 Background and theory

Figure 2.1. An adipocyte (a fat cell). The cell is almost completely filled by a large droplet of triglycerides. The smaller structure in the top is the nucleus. The cytoplasm where the metabolic insulin signaling pathway takes place, is only a thin rim below the plasma membrane.

Figure 2.2. Caveolae are flask-shaped invaginations in the plasma membrane. In adipocytes, the insulin receptors are located in caveolae, with the insulin binding α-subunits facing the exterior of the cell. The protein at the bottle-neck of the caveolae is the structural protein caveolin.

(21)

2.4 Systems, models and simulations 9

conformational change and rapid autophosphorylation of the receptor. In the case of metabolic IR signaling, the next step in the pathway is phosphorylation of a member of the insulin receptor substrate (IRS) family. Phosphorylated IRS then form a complex with, and activates phosphatidylinositol 3-kinase (PI3-K). The complex of phosphorylated IRS and activated PI3-K then propagates the signal downstream, and the signaling cascade ends with the translocation of the glucose transporter protein GLUT4 from an intra-cellular pool (i.e., no transcription is needed) to the plasma membrane.

GLUT4 transports glucose from the blood into the cytosol, thereby lowering blood glucose level. Eventually this will lead to reduced production of insulin in the pancreas. From a systems-perspective, this is a nice example of a negative feedback loop.

2.4

Systems, models and simulations

This section describes some basic notation and concepts in mathematical modeling and control theory.

2.4.1

Systems

A system is an object or a collection of objects, whose properties we wish to study [11], e.g. an electrical circuit, a cell, or a collection of cells. External signals that can be manipulated by us are called inputs, measurable signals from the system are called outputs. External signals which can not be manipulated by us are called disturbances. They are further divided into directly measurable disturbances and indirectly measurable disturbances, i.e., they are measurable only through their effect on the output.

2.4.2

Mathematical models

In a mathematical model, relations between physical quantities such as concen-tration, potential, distance, etc., are described by mathematical relations between variables. Models describing the dynamical properties of a system are commonly written as a system of differential equations. Such a system can be written in

state-space form, which means a system of coupled first-order ordinary differential

equations (ODEs). The measurable entities are given as a function h of the states and the inputs.

˙x = f (x, u) (2.1)

y = h(x, u) (2.2)

where the vector u(t) ∈ Rm is the input, the vector y(t) ∈ Rp is the output, the

(22)

10 Background and theory

Usually the equations depend on some parameters, θ ∈ Rd. To emphasize this

dependence, the system can be written:

˙x = f (x, u, θ) (2.3)

y = h(x, u, θ) (2.4)

If the systems of equations f and h are linear, the state-space description can be written in matrix form:

˙x = Ax + Bu (2.5)

y = Cx + Du (2.6)

A linear system in state-space form can also be described by a transfer function. With p as the differentiation operator, a system in state-space form can be written

˙x = px = Ax + Bu (2.7)

y = Cx + Du (2.8)

This yields

(pI − A)x(t) = Bu(t) ⇔ x(t) = (pI − A)−1Bu(t) (2.9)

y(t) = Cx(t) + Du(t) = (C(pI − A)−1B + D)u(t) (2.10)

If the initial values x(0) = 0, the transfer function G(p) can be defined by G(p) = C(pI − A)−1B + D (2.11)

A more general system description is differential-algebraic equations, DAEs. Apart from differential equations, DAEs can contain relations between undifferentiated variables. This is useful when there are conservations laws that apply to the situation the model aims to capture. E.g. in a closed chemical system, the number of atoms of each species is constant. A mathematical DAE description looks like

F ( ˙z, z, u) = 0 (2.12)

A special case of the DAE description when ˙z can be expressed explicitly is

˙z = f (z, u) (2.13)

which is the state-space description described above.

2.4.3

Simulations

A simulation of a model is the approximate solution of the model equations with the aid of numerical algorithms. State-space models are especially desirable due to the existence of efficient numerical methods for solving them. [11]

Loosely speaking, a time constant is the typical time, τ, needed for an interesting feature of a system to manifest. E.g. the metabolism of a drug could have a time

(23)

2.5 A mathematical model of insulin signaling 11

constant in the order of hours or days, while the physiological action of epinephrin has a time constant in the order of seconds. A system of differential equations with time constants of too different magnitudes (τmax/τmin ≥ 10 − 100) is called

stiff, and may need special solvers.

In Matlab there is a special ODE-solver for stiff differential equations based on numerical differentiation formulas (NDFs), called ode15s, which is used in this thesis.

2.5

A mathematical model of insulin signaling

The models in this thesis are based on the insulin signaling model proposed by Sedaghat et al. [12]. It presents a model of the signaling pathway from extracellular insulin to the translocation of GLUT4 receptors to the plasma membrane. This thesis is concerned only with the first steps of insulin signaling, from extracellular insulin to the formation of the IRS-1-YP/PI3-K complex. The corresponding part in the Sedaghat model can be divided into an insulin receptor subsystem and a postreceptor subsystem. They are described in some detail below.

2.5.1

Insulin receptor subsystem

The IR subsystem is composed of insulin binding kinetics and receptor internal-ization/dephosphorylation (see Figure 2.3).

The states in the insulin binding kinetics are called x1 – free insulin

x2 – free IR at the plasma membrane

x3 – IR bound to one molecule of insulin

x4 – phosphorylated IR bound to two insulin molecules

x5 – phosphorylated IR bound to one insulin molecule

Once an insulin molecule has bound to IR, the receptor is quickly phosphory-lated [12] (this is due to a conformational change in the IR protein which activates the intrinsic tyrosine kinase capacity of the IR β-subunits). The parameters k1,

k−1, k2 and k−2 are the association and dissociation rate constants for the first

and the second insulin molecule binding to its receptor, respectively. k3is the rate

constant for receptor autophosphorylation and k−3 is the receptor dephosphoryla-tion rate constant. k−3is multiplied by [PTP], a factor representing the activity of protein tyrosine phosphatases (PTPases) that can dephosphorylate the IR. These

(24)

12 Background and theory

Figure 2.3. Insulin receptor subsystem

are the state equations of the states x1-x5:

x1= insulin input (2.14) dx2 dt = k−1x3+ k−3[PTP]x5− k1x1x2+ k−4x6− k4x2 (2.15) dx3 dt = k1x1x2− k−1x3− k3x3 (2.16) dx4 dt = k2x1x5− k−2x4+ k ′ −4x7− k′4x4 (2.17) dx5 dt = k3x3+ k−2x4− k2x1x5− k−3[PTP]x5+ k ′ −4x8− k′4x5 (2.18)

The IR internalization is modeled by

x6 – free internal IR

x7 – internalized twice-bound phosphorylated IR

x8 – internalized once-bound phosphorylated IR

The internalization of phosphorylated IR is governed by the rate constants k′ 4

for endocytosis, and k′

−4 for exocytosis. The parameters k4 and k−4 determine

endocytosis and exocytosis rates of free IR. k5 and k−5 are rate constants for

synthesis and degradation of IR, where the synthesis of IR is a zero order rate constant. The parameter k6 describes dephosphorylation of endosomic IR. k6 is

(25)

2.6 Dimensional analysis 13 x6-x8: dx6 dt = k5− k−5x6+ k6[PTP](x7+ x8) + k4x2− k−4x6 (2.19) dx7 dt = k ′ 4x4− k′−4x7− k6[PTP]x7 (2.20) dx8 dt = k ′ 4x5− k′−4x8− k6[PTP]x8 (2.21)

2.5.2

Postreceptor signaling

Postreceptor signaling is described by

x9 – IRS

x10 – phosphorylated IRS

x11 – PI3-K

x12 – IRS/PI3-K complex

Insulin receptors at the plasma membrane can phosphorylate IRS to become IRS-P. Activated IRS can then form a complex with PI3-K. The phosphorylation of IRS is governed by the rate constant k7. The rate constant k7 is divided by

IRp, which is the concentration of phosphorylated surface receptors at maximal

insulin stimulation. IRS-P is dephosphorylated with the rate k−7 and the factor [PTP]. The complex of activated IRS and PI3-K forms with rate constant k8and

dissociates with the rate k−8 (see Figure 2.4). State equations of states x9-x12:

dx9 dt = k−7[PTP]x10− k7x9(x4+ x5)/(IRp) (2.22) dx10 dt = k7x9(x4+ x5)/(IRp) + k−8x12− k−7[PTP]x10− k8x10x11 (2.23) dx11 dt = k−8x12− k8x10x11 (2.24) dx12 dt = k8x10x11− k−8x12 (2.25)

Figure 2.5 shows both the IR subsystem and the postreceptor signaling.

2.6

Dimensional analysis

A dimension is a measurable property like time, distance or concentration. A unit is an arbitrarily chosen standard amount of a dimension, and a measured quantity is expressed as a multiple of a unit. A duration e.g., can be expressed as 30 minutes or as half an hour.

(26)

14 Background and theory

Figure 2.4. Postreceptor signaling

(27)

2.6 Dimensional analysis 15

For a model of a physical system to be valid, the dimensionality of both sides of the equations must be equal. Furthermore, only terms with the same dimensionality can be added or subtracted, and an argument to a transcendental function such as the exponential, logarithmic or a trigonometric function must be dimensionless. All physical models must hold independently of the system of units.

If F (x1, . . . , xn) = 0 is a model of a physical system with n variables, each variable

can be written as xi= x∗iˆxi where x∗i is a scalar multiple and ˆxi is a unit carrying

the dimension.

In a state-space description of a system, ˙x = f(x, u), each equation has the form dxi

dt = fi(x1, . . . , xn, u1, . . . , um) (2.26) and can thus be written as

d(x∗ ixˆi)

d(t∗τ ) = fi(x ∗

1xˆ1, . . . , x∗nˆxn, u1∗uˆ1, . . . , u∗muˆm) (2.27)

If each side of the equations is multiplied by τ ˆ

xi and the sizes of ˆxi, ˆui and τ are

chosen to some values, usually to scales inherent to the system, the result will be a system where all variables are dimension free. Such a system is said to be in

non-dimensional form.

Example 2.1

These are the equations for a Chemostat, example taken from [4]: dN dt =  KmaxC Kn+ C  N −F N V (2.28) dC dt = −α  KmaxC Kn+ C  N −F C V + F C0 V (2.29)

As a first step in non-dimensionalization of the equations, N and C are written as scalar multiples times their units, which yields

d(N∗N )ˆ d(t∗τ ) = KmaxC∗Cˆ Kn+ C∗Cˆ ! N∗N −ˆ F VN ∗Nˆ (2.30) d(C∗C)ˆ d(t∗τ ) = −α KmaxC∗Cˆ Kn+ C∗Cˆ ! N∗N −ˆ F C∗Cˆ V + F C0 V (2.31) Choosing τ = V F, ˆC = Kn and ˆN = Kn

ατ Kmax, and naming

α1 = τ Kmax α2 = τ F C0 V ˆC = C0 Kn

(28)

16 Background and theory

renders the following non-dimensional equations, where stars are dropped for sim-ple notation: dN dt = α1  C 1 + C  N − N (2.32) dC dt = −  C 1 + C  N − C + α2 (2.33)

The non-dimensional equations contain only two parameters compared to six in the original equations.

2.7

Observability

Loosely speaking, a state is observable if its value can be deduced from the input-output behavior of the system. Apart from the states, it is often interesting to know whether the value of a parameter can be identified, in principle, from the input-output behavior of a system.

2.7.1

Parameter identifiability

Identifiability is a property of a system where the parameters, θ, can be determined unambiguously from its input-output behavior. If the parameters in a system are considered as states with time derivatives zero, they can be incorporated in a state-space description and methods to decide on state observability can be used [2]. In an article by Ljung and Glad [10], an algorithm for deciding identifiability of differential-algebraic equations (DAEs) with polynomial nonlinearities is de-scribed. The algorithm has high complexity, and is not applicable to large prob-lems. Instead, an algorithm by Alexandre Sedoglavic [14] described in Section 2.7.3, can be used.

2.7.2

Linear systems

For linear systems, the following definition of observability can be given, from [6]:

Definition 2.1 The state x6= 0 is called non-observable if, when u(t) = 0, t ≥ 0

and x(0) = xthe output is y(t) = 0, t ≥ 0. The system is called observable if it

(29)

2.7 Observability 17

For linear systems there is also an easy way of checking the observability; a linear system, Σ

Σ : ˙x = Ax + Bu

y = Cx (2.34)

is observable if and only if the matrix

O =      C CA ... CAn−1      (2.35)

has full rank [6].

2.7.3

Nonlinear systems

For nonlinear systems, no equally simple method of determining observability ex-ists.

There is a probabilistic algorithm by Alexandre Sedoglavic [14], which decides nonlinear observability of a system of rational functions in polynomial time. Pa-rameters are treated as state-variables with time-derivative zero, ˙θ = 0, which makes this algorithm useful for parameter identifiability analysis. The algorithm determines the observable states, and gives the set of non-observable states with a high probability. It also reports a transcendence degree, which can be interpreted as the number of non-observable variables which need to be known in order to obtain an observable system [14].

The Sedoglavic’ algorithm makes no use of numerical parameter values. This means that a reported observability is structural, i.e., it exists some neighborhood in the parameter space where the states are observable.

There is a Maple implementation of the algorithm, available on the world wide web at [13].

To clarify the meaning of identifiability, a simple example with non-observable parameters is shown below:

Example 2.2

The simple system

dx dt = − p1 p2 x(t) (2.36) y(t) = x(t) (2.37)

has the solution

y(t) = x(t) = x0e−(

p1

(30)

18 Background and theory

In this case it is easy to see that the ratio p1

p2 can be observed but not the

param-eters them-selves. If either one of the paramparam-eters are known, the other becomes observable as well.

If this small system is analyzed with the Maple implementation of Sedoglavic’ algorithm, it yields the set of non-observable entities {p1, p2} and a transcendence

degree of one, as expected. If either p1 or p2 is known, the transcendence degree

drops to zero and there are no non-observable entities.

2.8

Controllability and observability gramians

A linear system, ˙x = Ax + Bu (2.39) y = Cx + Du (2.40) can be represented as ˙ξ = T AT−1ξ + T Bu (2.41) y = CT−1ξ + Du (2.42) by a change of variables ξ(t) = T x(t). The controllability energy function [7],

min u∈L2(−∞,0),x(−∞)=0,x(0)=x0 0 Z −∞ u(t) 2 dt = xT0Sx−1x0 (2.43)

describes the amount of input energy needed to bring the state vector of a system to x0, the matrix Sxis called the controllability gramian.

The observability energy function [7],

∞ Z 0 y(t) 2 dt = xT 0Oxx0 (2.44)

describes how the energy in the states contribute to the output. The matrix Ox

is called the observability gramian.

If input i is an impulse; ui(t) = δ(t) and x0= 0, the state vector becomes x(t) =

eAtB

i, where Bi is the ith column of B. Then the matrix Sx can be defined:

Sx = ∞

Z

0

(31)

2.8 Controllability and observability gramians 19

If T is chosen such that Sξ becomes diagonal

Sξ=      σ1 0 . . . 0 0 σ2 . . . 0 ... ... ... ... 0 0 . . . σn      (2.46)

where σ1 ≥ . . . ≥ σn, then the elements σi measures how much the states ξi are

affected by the inputs.

If there is no direct term (no matrix D in the state-space description), and the input u(t) = 0, the output becomes y(t) = CeAtx

0. Then the observability gramian

matrix Ox can be defined via: ∞ Z 0 y(t)Ty(t)dt = xT 0 ∞ Z 0 eATtCTCeAtdtx 0= xT0Oxx0 (2.47)

If T is chosen such that Oξ becomes diagonal, the importance of the different

components in ξ0 to the output can be seen.

It is possible to make a change of variables ξ(t) = T x(t) such that Sξ and Oξ are

equal to the same diagonal matrix Σ. This is called balancing of the gramians. If the gramians are balanced and the state variables and outputs are scaled so their variations are approximately equal, then the importance of the states ξi to

the input-output behavior of the model can be concluded. A state ξicorresponding

to a small value of σi is affected only slightly by the input and has low influence

on the output. In model reduction this means that a state variable ξiwith a small

value of σi can be eliminated without loosing much of the input-output qualities

of the system [6].

2.8.1

Empirical gramians

Empirical gramians is a concept introduced by Lall et al. [9] and later improved by Hahn and Edgar [7]. It is a generalization of linear gramians to the nonlinear case. If the empirical gramians are applied to a linear system they are reduced to the gramians described in the previous section. The following definition is the one by [9].

Before the definition of empirical gramians, the following sets need to be defined: Tn = {T1, . . . , Tr; Ti∈ Rn×n, TiTTi= I, i = 1, . . . , r} (2.48)

M = {c1, . . . , cs; ci∈ R, ci> 0, i = 1, . . . , s} (2.49)

En = {e

i; i = 1, . . . , n} (2.50)

where T is an arbitrary set of r orthogonal matrices, M is a set of positive con-stants and En are the standard unit vectors in Rn.

(32)

20 Background and theory

Definition 2.2 Let Tn, Ep and M be sets defined above. Then the empirical

controllability gramian Wc, and the empirical observability gramian Woare defined

by: Wc = r X l=1 s X m=1 p X i=1 1 rsc2 m ∞ Z 0 Φilm(t)dt (2.51) Wo = r X l=1 s X m=1 1 rsc2 m ∞ Z 0 TlΨlm(t)TlTdt (2.52)

Φilm(t) = (xilm(t) − xilm0 )(xilm(t) − xilm0 )T (2.53)

Ψlm

ij (t) = (yilm(t) − yilm0 )T(yilm(t) − yilm0 ) (2.54)

The empirical gramians can be used for model reduction similar to the linear case.

2.9

Parameter estimation

In parameter estimation numerical values of the model parameters, θ, that mini-mizes some criterion are sought.

Definition 2.3 The prediction error at time t is the difference between the system output, y(t), and the prediction, ˆy(t|θ):

ε(t, θ) = y(t) − ˆy(t|θ) (2.55)

If input-output data; y(1), . . . , y(N) are available, the deviation of a model with parameters θ from the data can be measured by

VN(θ) = 1 N N X t=1 ℓ(ε(t, θ)) (2.56)

where ℓ( · ) is a scalar valued function with minimum at zero. For a multi-output system, a quadratic criterion looks like

ℓ(ε) = 1 2ε

TΛ−1ε (2.57)

where Λ−1 is a symmetric positive semidefinite matrix that weighs together the

components of ε.

To get good agreement with data, parameters are chosen to minimize VN(θ):

ˆ

θN = arg min

(33)

2.10 Parameter sensitivity analysis 21

2.10

Parameter sensitivity analysis

Parameter sensitivity analysis is a procedure that investigates how much the dif-ferent parameters influence the output. If a parameter has a large effect on the output relative to the other parameters, the output is said to be sensitive to that parameter. Sensitivity analysis can be used to identify the most important pa-rameters in a model. If the output of a model is very sensitive to a parameter, it implies that the parameter can be estimated accurately, and if the output has a low sensitivity with respect to a parameter it means that a range of parameter values can give rise to very similar model outputs.

2.10.1

Derivatives of the prediction

One way of doing sensitivity analysis is by computing the derivatives of the pre-dictor with respect to the different parameters θ.

With an unbiased output-error model; y(t) = ˆy(t|θ) + ε(t, θ), the variance error in an estimate ˆθN can be approximated by

PN = E[(ˆθN − θ0)(ˆθN− θ0)T] ≈

1 NλNR

−1

N (2.59)

where θ0 is a parameter vector such that y(t) − ˆy(t|θ0) = ε(t, θ0) = white noise

with variance λ, and where λN and RN are estimated as

ˆ λN = 1 N N X t=1 ε2(t, ˆθ N) (2.60) ˆ RN = 1 N N X t=1 ψ(t, ˆθN)ψT(t, ˆθN) (2.61) where ψ(t|θ)T = d dθy(t|ˆˆ θN) (2.62) This means that ˆR is an approximation of the covariance matrix of the gradient of the prediction with respect to the parameters θ.

If the system is written on state-space form

˙x = f (x, u, θ) (2.63)

ˆ

y = h(x, u, θ) (2.64)

the derivatives ψ(t|θ)T = d

dθy(t|ˆˆ θN) can then be computed as

ψT(t|θ) = d dθy(t|θ) =ˆ (2.65) ∂ ∂θ[h(x(t, θ), u(t), θ)] = ∂ ∂xh(x, u, θ) ∂ ∂θx(t, θ) + ∂ ∂θh(x, u, θ) (2.66)

(34)

22 Background and theory With z(t, θ) = ∂ ∂θx(t, θ) we have ∂ ∂tz(t, θ) = ∂ ∂t ∂ ∂θx(t, θ) = ∂ ∂θ ∂ ∂tx(t, θ) = ∂ ∂θf (x(t, θ), u(t), θ) (2.67) = ∂ ∂xf (x, u, θ)z(t, θ) + ∂ ∂θf (x, u, θ) (2.68) if x(t, θ) has continuous second derivatives.

(35)

Chapter 3

Experimental data

This chapter describes the data sets used in this thesis.

3.1

Time series

A time series is a collection of data points, usually measured with uniform sample time (though the time series utilized in this thesis are not uniformly sampled). The time series employed in this thesis were collected from human and rat adipo-cytes. The measured signals were the relative increase in phosphorylation de-gree of IR and IRS-1 after stimulation with a step of insulin. The experimental data from humans comes from Peter Strålfors, at the Department of Biomedicine and Surgery, Linköpings Universitet. The experimental data from rats originate from [5].

3.2

Excitation signal

In both the case of human and rat adipocytes, the cells were treated with 100 nM of insulin, i.e. approximately a 100-fold the physiological concentration. The concentration was chosen large in order to get a maximum response from the cells. Since the insulin was in huge excess, the insulin concentration could be considered constant during the experiment. Therefore the excitation signal is viewed as a step switched on at time t = 0.

(36)

24 Experimental data

3.3

Output equations

Since the phosphorylation degree of IR and IRS-1 can be measured only relatively, the output equations contain proportionality constants ky1 and ky2

h(x, θ) = ky1·

(sum of phosphorylated IR-states) IRtot

ky2·

(sum of phosphorylated IRS-1-states) IRS−1tot

!

(3.1)

where IRtot = (sum of all IR-states) and IRS-1tot = (sum of all IRS-1-states).

These are the output equations used throughout the thesis.

3.4

Data sets

Two data sets were used in this thesis. One set was collected from rat white adipocytes, and one from human adipocytes. They are referred to as DataHuman and DataRat.

The standard error of the mean (SEM) statistic is defined as

Definition 3.1 SEM = σ √

N, where σ is the standard deviation and N is the

number of samples.

DataHuman consists of three time series collected at twelve time-points, sampled at 0, 0.9, 1.15, 1.3, 2.3, 3.3, 5.3, 7, 10, 15, 20 and 30 minutes. The means and the SEM statistics are shown in Figure 3.1.

DataRat consists of four time series collected at fourteen time-points, sampled at 0, 0.5, 0.75, 1, 2, 3, 5, 7, 10, 15, 20, 30, 45 and 60 minutes. Means and SEM are given in Figure 3.2. Only samples up to t = 15 min are utilized in this work. After that time, other cellular processes not captured by the models in this thesis may begin to be important.

3.4.1

Pre-processing of data

The data sets were pre-processed before usage. Each time series was divided by its largest value, and averages were calculated at each sample. The SEM statistic was computed for each sample instant as a measure of the variability of the means.

(37)

3.4 Data sets 25 0 5 10 15 20 25 30 0 20 40 60 80 100 Time [min] % of max IR−P / IR−tot 0 5 10 15 20 25 30 0 20 40 60 80 100 Time [min] % of max IRS1−P / IRS1−tot

Figure 3.1. Relative increase in phosphorylation degree of IR and IRS1 in human adipocytes treated with 100 nM insulin.

0 10 20 30 40 50 60 0 20 40 60 80 100 Time [min] % of max IR−P / IR−tot 0 10 20 30 40 50 60 0 20 40 60 80 100 Time [min] % of max IRS1−P / IRS1−tot

Figure 3.2.Relative increase in phosphorylation degree of IR and IRS1 in rat adipocytes treated with 100 nM insulin.

(38)
(39)

Chapter 4

Results

This chapter presents the results from modeling, simulation, observability analysis and parameter sensitivity analysis.

4.1

Model structures

Starting with the model structure described in the background, Section 2.5, a first simplification was to set the parameters k5 and k−5 equal to zero. These

param-eters govern the rate of IR synthesis and degradation, and the simplification of setting them to zero is justified by the fact that the amount of IR can be considered constant on the time-scale of the experiments. Synthesis involves transcription of DNA, as well as translation of mRNA, which are relatively slow processes. Previous models do not include the phosphorylation of IRS by internalized IR. This possibility was added through the parameter k′

7 (see Figure 4.2 and cf. the

original Sedaghat model shown in Figure 4.1, where k′

7 is missing and IRS is

phosphorylated from the plasma membrane only).

An inspection of the equations in the model by Sedaghat et al. revealed that the parameters [PTP] and IRp never could be estimated from data since they always

appear in a product or a quotient with another parameter, and always with the same other parameter. The parameters k−3, k6 and k−7 are always multiplied

by [PTP], and k7 is always divided by IRp. Therefore the products k−3· [PTP],

k6· [PTP], k−7· [PTP] and k7/IRp were simply called k−3, k6, k−7 and k7 to get

easier notation.

The state representing the input, x1 was renamed u to distinguish it from the

variables, and all states x2 to x12 were renamed x1 to x11. This model consists

of 11 states and was given the name M1 (see Figure 4.2). These are the model 27

(40)

28 Results equations of M1: dx1 dt = k−1x2+ k−4x5+ k−3x4− k1x1u − k4x1 (4.1) dx2 dt = k1x1u − k−1x2− k3x2 (4.2) dx3 dt = k2x4u + k ′ −4x6− k−2x3− k′4x3 (4.3) dx4 dt = k3x2+ k−2x3+ k ′ −4x7− k2x4u − k−3x4− k4′x4 (4.4) dx5 dt = k4x1+ k6(x6+ x7) − k−4x5 (4.5) dx6 dt = k ′ 4x3− k′−4x6− k6x6 (4.6) dx7 dt = k ′ 4x4− k′−4x7− k6x7 (4.7) dx8 dt = k−7x9− k7x8(x3+ x4) − k ′ 7x8(x6+ x7) (4.8) dx9 dt = k7x8(x3+ x4) + k ′ 7x8(x6+ x7) + k−8x11− k−7x9− k8x9x10 (4.9) dx10 dt = k−8x11− k8x9x10 (4.10) dx11 dt = k8x9x10− k−8x11 (4.11)

where u is the input. Given the measurements described in Section 3.3, the output equations become: y1 = ky1· x3+ x4+ x6+ x7 x1+ x2+ x3+ x4+ x5+ x6+ x7 (4.12) y2 = ky2· x9+ x11 x8+ x9+ x11 (4.13)

The next simplification was the merging of the states in the plasma membrane rep-resenting IR bound to one insulin molecule and IR bound to two insulin molecules, i.e., the states x3 and x4. This simplification is motivated by the fact that the IR

becomes phosphorylated when the first insulin molecule is bound, and that phos-phorylated IR is responsible for transmitting the signal downstream. The merged state was called x3. In the same manner, the states in the endosome, x6 and x7,

(41)

4.1 Model structures 29

Figure 4.1. The first steps of the Sedaghat model.

Figure 4.2. Model M1: Synthesis and degradation of IR is neglected, and phosphory-lation of IRS by internalized IR is taken in account.

(42)

30 Results

Figure 4.3.Model M2: The states with IR bound to one and two insulin molecules are merged to become the new states x3 and x4.

(see Figure 4.3). Model equations of M2: dx1 dt = k−1x2+ k−3x3+ k−4x5− k1x1u − k4x1 (4.14) dx2 dt = k1x1u − k−1x2− k3x2 (4.15) dx3 dt = k3x2+ k ′ −4x4− k−3x3− k′4x3 (4.16) dx4 dt = k ′ 4x3− k−4′ x4− k6x4 (4.17) dx5 dt = k4x1+ k6x4− k−4x5 (4.18) dx6 dt = k−7x7− k7x3x6− k ′ 7x4x6 (4.19) dx7 dt = k7x3x6+ k ′ 7x4x6+ k−8x9− k−7x7− k8x7x8 (4.20) dx8 dt = k−8x9− k8x7x8 (4.21) dx9 dt = k8x7x8− k−8x9 (4.22)

The output equations become: y1 = ky1· x3+ x4 x1+ x2+ x3+ x4+ x5 (4.23) y2 = ky2· x7+ x9 x6+ x7+ x9 (4.24)

Finally, it was recognized that the magnitude of the parameter k3was much higher

(43)

4.1 Model structures 31

course of events that occurs quickly compared to events governed by other param-eters. Actually, k3 is the rate of receptor autophosphorylation, which is indeed

a fast reaction [12]. Assuming that the receptor autophosphorylation occurs in-stantly once insulin is bound, the state x2 and the parameters k3 and k−3 were

removed, and the states x3 to x9 were renamed x2 to x8. The resulting model

structure was called M3 (see Figure 4.4). These are the equations describing the model structure M3: dx1 dt = k−1x2+ k−4x4− k1x1u − k4x1 (4.25) dx2 dt = k1x1u + k ′ −4x3− k−1x2− k′4x2 (4.26) dx3 dt = k ′ 4x2− k−4′ x3− k6x3 (4.27) dx4 dt = k4x1+ k6x3− k−4x4 (4.28) dx5 dt = k−7x6− k7x2x5− k7′x3x5 (4.29) dx6 dt = k7x2x5+ k ′ 7x3x5+ k−8x8− k−7x6− k8x6x7 (4.30) dx7 dt = k−8x8− k8x6x7 (4.31) dx8 dt = k8x6x7− k−8x8 (4.32)

Since no synthesis or degradation of any species is included in the model, the amount of each species is conserved.

IRtot = x1+ x2+ x3+ x4= const (4.33)

IRStot = x5+ x6+ x8= const (4.34)

PI3Ktot = x7+ x8= const (4.35)

Given the measurements described in Section 3.3, the output equations become: y1 = ky1· x2+ x3 x1+ x2+ x3+ x4 (4.36) y2 = ky2· x6+ x8 x5+ x6+ x8 (4.37)

4.1.1

Empirical gramians

Empirical gramians for the three models above were computed, with parameter values from [12]. An inspection of the singular values supports that no state vital to the input-output behavior of the system was removed during the simplifications, since the most energy-rich states are conserved. The singular values of the models are listed in Table 4.1.

(44)

32 Results

Figure 4.4. Model M3: The fast dynamics of the receptor autophosphorylation is assumed to be a static reaction, i.e., the receptor is phosphorylated instantly when insulin is bound. M1 M2 M3 σ1 1.96 · 101 1.91 · 101 1.91 · 101 σ2 1.81 · 100 1.80 · 100 1.80 · 100 σ3 5.60 · 10−3 5.10 · 10−3 5.10 · 10−3 σ4 2.10 · 10−3 4.01 · 10−5 4.05 · 10−5 σ5 3.49 · 10−5 5.84 · 10−6 1.62 · 10−8 σ6 3.65 · 10−6 8.06 · 10−9 6.97 · 10−9 σ7 8.95 · 10−9 5.13 · 10−10 4.71 · 10−15 σ8 6.44 · 10−9 2.90 · 10−15 1.62 · 10−16 σ9 4.29 · 10−9 3.80 · 10−16 σ10 4.80 · 10−15 σ11 1.82 · 10−17

Table 4.1. Empirical gramians of M1, M2 and M3. Each σi represents the relative

(45)

4.1 Model structures 33

4.1.2

Some observations

To decide whether IRS is phosphorylated mainly by plasma membrane bound IR or by endosomic IR, the parameters k7 and k′7 are interesting, since they govern

the rate of the phosphorylation. Actually it would suffice to know if the ratio k7

k′ 7

is less than or greater than one to decide whether phosphorylation by endosomic IR occurs faster than by plasma membrane bound IR or vice versa.

This insight, that the ratio k7

k′

7 is sufficient to decide the rate of phosphorylation

was used when doing the non-dimensionalization, where the ratio appears as the parameter α8.

Another observation is that the first four states, x1, . . . , x4 form a closed

subsys-tem, i.e., none of their equations depend on any of the states x5, . . . , x8. This is

advantageous, because parameter estimation and sensitivity analysis can be done for the isolated subsystem, giving shorter computation time.

Yet another observation is that the autonomous system obtained when the input to the IR-subsystem, (x1, . . . , x4), is constant, is a linear, time invariant (LTI)

system. Thus the transfer function of the IR-subsystem can be computed, and tools from linear control theory can be used to analyze its behavior.

4.1.3

Non-dimensionalization of M3

The equations of the model M3 were non-dimensionalized and re-parameterized to obtain the ratio k7

k′

7 as one parameter. The equations of M3 were rewritten with

xi = x∗i· ˆxi, u = u∗· ˆu and t = t∗· τ where x∗i, u∗ and t∗ are dimension-free.

Since the ratio k7

k′

7 was interesting, τ was chosen as τ =

1 k′

7xˆ. All states were scaled

by the same constant; ˆxi= ˆx, i = 1, . . . , n. Finally ˆu was chosen as ˆu = k

′ 7xˆ

k1 , since

then the number of parameters could be reduced by two. The new parameters were defined as ratios of the old ones:

α1= k−1 k′ 7xˆ α6= k6 k′ 7xˆ α2= k−4 k′ 7xˆ α7= k−7 k′ 7xˆ α3= k4 k′ 7xˆ α8= k7 k′ 7 α4= k′ −4 k′ 7xˆ α9= k−8 k′ 7xˆ α5= k′ 4 k′ 7xˆ α10= k8 k′ 7

(46)

34 Results

The resulting system of equations has only 10 parameters compared to 12 before the non-dimensionalization and the interesting ratio, k7

k′

7, now appears as one

pa-rameter, α8. The problem of finding whether IRS is phosphorylated more rapidly

by plasma membrane bound or endosomic IR, can thus be restated as finding whether the parameter α8is less than or greater than one.

These are the equations describing the non-dimensional system, the stars on the state variables are dropped for notational convenience:

dx1 dt = α1x2+ α2x4− x1u − α3x1 (4.38) dx2 dt = x1u + α4x3− α1x2− α5x2 (4.39) dx3 dt = α5x2− α4x3− α6x3 (4.40) dx4 dt = α3x1+ α6x3− α2x4 (4.41) dx5 dt = α7x6− α8x2x5− x3x5 (4.42) dx6 dt = α8x2x5+ x3x5+ α9x8− α7x6− α10x6x7 (4.43) dx7 dt = α9x8− α10x6x7 (4.44) dx8 dt = α10x6x7− α9x8 (4.45) The output equations are the same as those in model M3:

y1 = ky1· x2+ x3 x1+ x2+ x3+ x4 (4.46) y2 = ky2· x6+ x8 x5+ x6+ x8 (4.47) They were already dimension-free, and thus unaffected by the non-dimension-alization procedure.

This non-dimensional model structure was called M3NonDim.

4.2

Structural observability

The structural observability of the model structures M1, M2, M3, M3NonDim and the IR-subsystem of M3 was analyzed using the Maple implementation of the algorithm by Alexandre Sedoglavic [14]. The results are collected in Table 4.2. The model structures M1, M2 and M3 all have transcendence degree 2, meaning that two of the non-observable entities need to be known in order to achieve an observable system. The non-observable parameters k7, k′7and k8are also the same

(47)

4.2 Structural observability 35

Model Tr. degree Non-observable entities M1 2 k7, k7′, k8, x1, x2, x3, x4 x5, x6, x7, x8, x9, x10, x11 M2 2 k7, k7′, k8, x1, x2, x3 x4, x5, x6, x7, x8, x9 M3 2 k7, k7′, k8, x1, x2, x3 x4, x5, x6, x7, x8 M3NonDim 1 α10, x5, x6, x7, x8 M3 IR-sub 1 x1, x2, x3, x4

Table 4.2. Results from observability analysis with Sedoglavic’s algorithm, measure-ments as in Section 3.3.

The model structure M3NonDim has a transcendence degree of one, meaning that only one of the entities a10, x5, x6, x7, x8 needs to be known in order to obtain

an observable system. In the IR-subsystem of M3, all parameters but no state is observable.

Unfortunately, the interesting parameters k7and k′7are among the non-observable

entities in the model structures M1, M2 and M3. But the analysis of the model structure M3NonDim reveals that the ratio k7

k′

7 = α8 is identifiable, which is

suffi-cient for the problem addressed in this thesis.

The model structure M3 was further investigated to reveal which combinations of the non-observable entities k7, k′7, k8, x1, x2, x3, x4, x5, x6, x7, x8 that need to be

known in order to achieve observability of the system. Different combinations of parameters and states were assumed known, and the Sedoglavic’ algorithm was run again to see the effects on the observability. The results are collected in Table 4.3.

The non-observable entities of M3 can be divided into two sets, A = {k7, k′7,

x1, x2, x3, x4} and B = {k8, x5, x6, x7, x8} (where A contains the states in the

IR-subsystem of M3). The results show that when at least one entity from each of the sets A and B are known simultaneously, observability is obtained. For M3NonDim, it suffices to know any one of the non-observable entities to achieve observability. It is possible to acquire an intuitive understanding about why the parameters k7,

k′

7 and k8 are non-observable. An inspection of the equations of M3 reveals that

k7, k7′ and k8 all occur in terms that are quadratic in some states, (e.g. k7x2x5)

and all other parameters are multiplied by only one state.

A partitioning of the states into two sets, xA= {x1, x2, x3, x4} and xB = {x5, x6,

x7, x8} can be made by observing that xA is the IR-subsystem of M3, and xB

is the post-receptor signaling of M3. The IR-subsystem is closed, i.e., none of its state equations depend on any state outside the subsystem. The structural observability results suggests a partitioning of the states of M3 that coincides with A and B above.

(48)

36 Results

Known entities Tr. degree Non-observable entities k7 1 k8, x5, x6, x7, x8 k′ 7 1 k8, x5, x6, x7, x8 k8 1 k7, k′7, x1, x2, x3, x4 k7, k′7 1 k8, x5, x6, x7, x8 k7, k8 0 k′ 7, k8 0 Any of x1, . . . , x4 1 k8, x5, x6, x7, x8 Any of x5, . . . , x8 1 k7, k7′, x1, x2, x3, x4 One of x1, . . . , x4 0 and one of x5, . . . , x8 One of x1, . . . , x4 0 and k8 One of x5, . . . , x8 1 k7, k7′, x1, x2, x3, x4 and k8 One of x5, . . . , x8 0 and one of k7, k7′ One of x1, . . . , x4 1 k8, x5, x6, x7, x8 and any of k7, k7′

Table 4.3. Results from observability analysis of M3 with Sedoglavic’s algorithm, with different combinations of parameters and states assumed known.

is a quotient of states from xAonly, and the second output is a quotient of states

from xB. This means that the absolute levels of the states are unknown, i.e., the

states in each set xA and xB can only be determined up to a constant. If the

constants, say γ1for the states in xA and γ2 for the states in xB, are included in

the state equations, and the same dynamics is required, only the parameters k7,

k′

7 and k8are affected, since the constants cancel for all other parameters.

For example, the third state equation of M3 dx3

dt = k

4x2− k′−4x3− k6x3 (4.48)

contains the observable parameters k′

4, k−4′ and k6. If the states are scaled by

γ1 (no scaling by γ2 is needed since x3is a part of the closed IR-subsystem), the

equation becomes

d(γ1x3)

dt = k

4γ1x2− k′−4γ1x3− k6γ1x3 (4.49)

where γ1 can be eliminated on both sides of the equation. Thus the dynamics of

the equation remain unchanged when the states are scaled by a constant. The sixth state equation is

dx6

dt = k7x2x5+ k

(49)

4.3 Simulations 37

and contains, besides the observable parameters k−7and k−8, all the non-observable parameters k7, k′7 and k8. When the states are scaled by γ1 and γ2 respectively,

the equation becomes d(γ2x6) dt = k7γ1γ2x2x5+ k ′ 7γ1γ2x3x5+ k−8γ2x8− k−7γ2x6− k8γ22x6x7 (4.51) ⇐⇒ dx6 dt = k7γ1x2x5+ k ′ 7γ1x3x5+ k−8x8− k−7x6− k8γ2x6x7 (4.52)

This means that k7 and k′7 must be scaled by γ11, and k8 must be scaled by

1 γ2

to maintain the same dynamics of the system. Thus k7, k7′ and k8 are affected

by the state scaling, and as long as the scalings γ1 and γ2are unknown, they are

non-observable. In the ratio k7

γ1/

k′ 7

γ1 the constant γ1 will cancel, making the ratio

k7

k′

7 unaffected by

the state scaling and hence observable.

4.2.1

Measurements that give observability

Another way of obtaining observability is by changing the measurements. If the output is of the form

y1 = (sum of phosphorylated IR-states) (4.53)

y2 = (sum of phosphorylated IRS-states) (4.54)

i.e., if the amount of phosphorylated IR and IRS can be measured on an absolute scale, the systems become observable. In the model structures M3 and M3NonDim, this would yield the output equations:

y1 = (x2+ x3) (4.55)

y2 = (x6+ x8) (4.56)

Results from Sedoglavic’s algorithm for some model structures with the measure-ments above are collected in Table 4.4.

Model Transcendence degree

M3 0

M3NonDim 0

M3 IR-sub 0

Table 4.4. Results from observability analysis with Sedoglavic’s algorithm, when IR-P and IRS1-P are measured on an absolute scale.

4.3

Simulations

The model structure M3 was implemented in Matlab. Since M3 still captures dynamics with different time constants (τmax/τmin≈ 104), the ODE-solver ode15s

(50)

38 Results 0 5 10 15 −0.2 0 0.2 0.4 0.6 0.8 1 1.2 Time [min] arb. units IR−P / IR−tot

Figure 4.5. Relative increase in IR phosphorylation in human adipocytes treated with 100 nM insulin. The circles represent the mean, and the plus signs represent the SEM statistic. The solid line is the simulated output from model M3.

for stiff equations was chosen for the integration.

The units of the concentrations were also changed from M to pM to get more “manageable numbers” in the simulations. This affected the initial concentrations of the states, the input u, and the parameters k1, k7, k′7 and k8 since they have

the unit M−1min−1.

4.3.1

Comparison with experimental data

The simulated output of IR phosphorylation was plotted together with the IR data from human and rat adipocytes. In the simulations, parameter values and initial conditions were initially chosen from the article by Sedaghat et al. [12].

The simulation of IR phosphorylation and the human IR time series are seen in Figure 4.5. The circles represent the mean of the measurement series and the plus signs represent the SEM statistic.

The output from M3 in Figure 4.5 shows little agreement with the experimental data; for instance, the output from M3 has no overshoot, but behaves like the step response of a first-order system.

(51)

4.3 Simulations 39

4.3.2

Transfer function of the IR-subsystem

To understand the origin of the overshoot expressed by the data, the model struc-ture M3 was further investigated. The equations (4.25)-(4.28), together with equa-tion (4.36) can be written as

˙x = A(u)x (4.57)

y = h(x) (4.58)

x(0) = x0 (4.59)

The system is linear in the states, and by holding the input u constant, an au-tonomous LTI system is achieved:

˙x = ˜Ax (4.60)

y = h(x) (4.61)

x(0) = x0 (4.62)

where ˜A incorporates the constant u. An examination of the output equation of the IR-subsystem, (4.36), together with equation (4.33) reveals that the denominator is constant. Thus the output equation can be rewritten as a linear combination of the phosphorylated IR-states:

y1 = ky1· x2+ x3 x1+ x2+ x3+ x4 (4.63) ⇐⇒ y1 = ky1· x2+ x3 IRtot (4.64) ⇐⇒ y1 = k′y1· (x2+ x3) (4.65) where k′ y1 = ky1/IRtot.

With the output in matrix form, the equations of the autonomous system can be organized as ˙x = A(u)x (4.66) y = Cx (4.67) x(0) = x0 (4.68) where C = (0, k′ y1, k ′

y1, 0). The system can be written equivalently as

˙x = ˜Ax + x0v(t) (4.69)

y = Cx (4.70)

v(t) = δ(t) (4.71)

x(0−) = 0 (4.72)

where δ(t) is an impulse that drives the system to begin at x0. Since the new

initial values are x(0−) = 0, the transfer function, Gu(s), of the system above can

be computed using the formula

(52)

40 Results

with A = ˜A, B = x0, C as above, and D = 0.

The resulting transfer function, Gu(s) has two zeros and four poles. Gu(s) has a

pole in the origin (i.e., contains a factor 1

s) due to moiety conservation of the IR,

i.e., x1+ x2+ x3+ x4= const. This can be seen by doing the variable substitution,

z1= x1+ x2+ x3+ x4 (4.74)

z2= x2 (4.75)

z3= x3 (4.76)

z4= x4 (4.77)

and noting that ˙z1= 0, which gives a row of zeros in the new A-matrix, implying

an eigenvalue equal to zero, i.e., a pole in the origin.

Instead of studying the impulse response of Gu(s), the step response of the system

G(s) = sGu(s) was studied. The transfer function G(s) has two zeros and three

poles, which can give rise to an overshoot in the step response due to the zeros. To see whether the internalization of the IR was necessary to obtain the overshoot, the IR-subsystem of M3 was further decomposed into insulin binding kinetics and receptor internalization/recycling. The insulin binding kinetics are represented by the states x1 and x2, with state equations

dx1

dt = k−1x2− k1ux1 (4.78) dx1

dt = k1ux1− k−1x2 (4.79) and the output equation

y1= k′y1· x2 (4.80)

A transfer function was computed according to the procedure above:

Gu(s) = k1k

′ y1x1(0)

s(k1+ k−1+ s)

(4.81)

Gu(s) also has a pole in the origin, due to the moiety conservation. The step

response of G(s) = sGu(s) was computed,

G(s) = k1k

′ y1x1(0)

k1+ k−1+ s

(4.82)

which is a transfer function that lacks zeros and has no complex poles. It cannot therefore give rise to an overshoot. Hence, internalization of the IR is needed to reproduce the overshoot of the phosphorylated IR in the experimental data in the model structure M3. Simulations which recreated the overshoot suggest that the IR internalization must take place much faster than in previous models.

(53)

4.3 Simulations 41 0 5 10 15 −0.5 0 0.5 1 1.5 2 Time [min] arb. units IR−P / IR−tot 0 5 10 15 −0.5 0 0.5 1 1.5 2 Time [min] arb. units IRS1−P / IRS1−tot

Figure 4.6. M3 with fast internalization, partly fitted to human data.

4.3.3

Recreating the overshoot of the phosphorylated IR

It was desirable to recreate the overshoot in the experimental data to find a re-alistic working point of the model, i.e., to find a set of parameters making the model behave like something that could be observed experimentally. Therefore several simulations with different parameter values were performed. The choice of parameter values was based on some intuition for the model. The most impor-tant change was the increased rate of receptor internalization and recycling which yielded an overshoot in the first output (the phosphorylation degree of the IR). The parameters of the IR-subsystem of M3 are structurally identifiable according to the analysis in Section 4.2. Therefore, a parameter estimation of M3IRsub (using the Matlab fminsearch command) was performed. This was performed for both human and rat time series. The results are seen in the Figures 4.6 and 4.7. The parameter values used to obtain Figure 4.6 and Figure 4.7 can be found in appendix A.

For DataHuman, the model could also reproduce the phosphorylation degree of IRS, as seen in Figure 4.6. Important to note is that the values of the parameters governing the postreceptor signaling were not found by any estimation algorithm due to the lack of identifiability. They are to some extent chosen arbitrarily, and only show that the model structure is capable of explaining the data.

For DataRat, no successful recreation of the second output was found. The time series of IRS phosphorylation in rat adipocytes display a different pattern

(54)

com-42 Results 0 5 10 15 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Time [min] arb. units IR−P / IR−tot

Figure 4.7. M3 with fast internalization, fitted to rat data.

pared to the human time series. Surprisingly, the first peak comes before the peak in IR phosphorylation, even though IRS is located downstream of IR in the signaling pathway.

Figure 4.8 shows how the states in the IR-subsystem of M3 change with time. The first state represents the free IR at the plasma membrane, which decreases rapidly when insulin is added. The second state represents the activated IR at the plasma membrane. It has a sharp peak when insulin is added, and then stays on a low constant level. The third state represents the activated internalized IR. It has a broader peak following the one in the previous state. The estimated parameters in the IR-subsystem resulted in a larger amount of internalized IR-P compared to previous models. The fourth state, the pool of free internal IR, rises from its initially low level and attains a steady state after about five minutes.

4.4

Sensitivity analysis

To find out whether the parameter α8 (in the model structure M3NonDim,

cor-responding to k7

k′

7) could be estimated from data, a more practical measure of

identifiability than the structural observability examined in Section 4.2 was desir-able.

To this end, the variance error matrix of the parameters, P (θ), was estimated using the methods in Section 2.10.1.

Figur

Updating...

Referenser

Updating...

Relaterade ämnen :