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
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
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 ISBN — ISRN 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
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.
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.
Contents
1 Introduction 1 1.1 Problem formulation . . . 1 1.2 Objectives . . . 1 1.3 Thesis disposition . . . 2 1.4 Abbreviations . . . 22 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
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
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.
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
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
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.
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).
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
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.
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
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
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
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
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 byx9 – 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.
14 Background and theory
Figure 2.4. Postreceptor signaling
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
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 x∗6= 0 is called non-observable if, when u(t) = 0, t ≥ 0
and x(0) = x∗ the output is y(t) = 0, t ≥ 0. The system is called observable if it
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
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
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.
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
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)
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.
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.
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.
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.
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
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,
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.
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
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.
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
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
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
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.
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
′
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
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.
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
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.
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
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.