Computational Modeling and Validation of Biopathways Koh Yeow Nam, - - PDF document

Computational Modeling and Validation of Biopathways

Koh Yeow Nam, Geoffrey A dissertation proposal for the degree of Doctor of Philosophy 24th November 2005 NUS Graduate School for Integrative Sciences and Engineering National University of Singapore

ABSTRACT

The proposed research centers on the development of a hybrid modeling framework for the representation, simulation and validation of biopathways. It aims to provide not only the functional view of the cell components but also their dynamics. Most studies on cell functions focus only on either signaling pathways, gene regulatory networks or metabolic

• pathways. The distinguishing feature of this work will then be to allow these different types
• f pathways to communicate with one another so that they can function as a whole.

Our approach will encompass modeling, data interpretation, parameter estimation and model validation. In the modeling of biopathways, we will not restrict ourselves to a single type of representation. Instead, a combination of different computational models may be used, depending on the type of pathway that is to be modeled. Apart from developing the modeling and simulation framework, proper metrics and validation techniques will also be used to assert the correctness of the models. The developed techniques will then be tested for feasibility using examples from biological experiments. As part of our preliminary studies, we model two well-studied pathways - the Wnt and the Akt signaling pathways. Using these two pathways, we will identify the key issues that arise and our strategy for addressing them.

CONTENTS

Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ii List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iv List of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . v

• 1. Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

1 1.1 Objectives . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.2 Rationale . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.3 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3

• 2. Background - Concepts in Biology

. . . . . . . . . . . . . . . . . . . . . . . . . . 4 2.1 The Cell . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 2.1.1 The Central Dogma . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 2.1.2 Transcription and Translation . . . . . . . . . . . . . . . . . . . . . . 6 2.1.3 Chemical Reactions . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 2.2 Biological Pathways . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 2.2.1 Metabolic Pathway . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 2.2.2 Signaling Pathway . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 2.2.3 Gene Regulatory Network . . . . . . . . . . . . . . . . . . . . . . . . 10 2.3 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11

• 3. A Review of Current Techniques

. . . . . . . . . . . . . . . . . . . . . . . . . . . 12 3.1 Quantitative or Qualitative Modeling? . . . . . . . . . . . . . . . . . . . . . 12 3.2 Kinetic Modeling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 3.2.1 Michaelis-Menten Kinetics . . . . . . . . . . . . . . . . . . . . . . . . 13 3.2.2 Canonical Modeling and S-Systems . . . . . . . . . . . . . . . . . . . 15 3.3 Petri Nets . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 3.4 Hybrid Automata . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 3.5 Process Calculi . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 3.6 XS-System . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 3.7 Piecewise-affine differential equations . . . . . . . . . . . . . . . . . . . . . . 22 3.8 Other Models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25

Contents ii

3.8.1 Bayesian Networks . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 3.8.2 Boolean Networks . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 3.8.3 Association Networks . . . . . . . . . . . . . . . . . . . . . . . . . . 27 3.9 Standards and Tools . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 3.9.1 Systems Biology Markup Language . . . . . . . . . . . . . . . . . . . 28 3.9.2 Systems Biology Workbench . . . . . . . . . . . . . . . . . . . . . . . 30 3.9.3 Cell Illustrator . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 3.10 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34

• 4. Problem Formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

35 4.1 Framework for Cell Modeling . . . . . . . . . . . . . . . . . . . . . . . . . . 35 4.2 Data Collection and Interpretation . . . . . . . . . . . . . . . . . . . . . . . 35 4.3 Validation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 4.3.1 Metrics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 4.3.2 Methods and Paradigms for Validation . . . . . . . . . . . . . . . . . 37 4.4 Model Representation and Integration . . . . . . . . . . . . . . . . . . . . . 38 4.5 Issues to be addressed by proposed research . . . . . . . . . . . . . . . . . . 39

• 5. Preliminary Studies

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 5.1 The Wnt Pathways . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 5.1.1 The Wnt Canonical Pathway . . . . . . . . . . . . . . . . . . . . . . 40 5.1.2 The Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 5.1.3 Downstream Gene Regulatory Network . . . . . . . . . . . . . . . . 41 5.1.4 Stem Cells and Differentiation . . . . . . . . . . . . . . . . . . . . . 41 5.1.5 Issues and Non canonical pathways . . . . . . . . . . . . . . . . . . . 43 5.2 AKT and ERK Pathways . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 5.2.1 Parameter Estimation . . . . . . . . . . . . . . . . . . . . . . . . . . 47 5.2.2 Pathway Topology . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48 5.2.3 Topological Ordering . . . . . . . . . . . . . . . . . . . . . . . . . . . 49 5.2.4 Rank . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 5.2.5 Evolutionary Strategy . . . . . . . . . . . . . . . . . . . . . . . . . . 51 5.2.6 Experimental Data . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52 5.2.7 Simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 5.2.8 Model Validation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58 5.3 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61

• 6. Future Plans

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62

LIST OF FIGURES

2.1 Complexity spectrum of biological system . . . . . . . . . . . . . . . . . . . 4 2.2 The Central Dogma of Molecular Biology . . . . . . . . . . . . . . . . . . . 5 2.3 The types of biological pathways . . . . . . . . . . . . . . . . . . . . . . . . 8 2.4 The glycolysis metabolic pathway . . . . . . . . . . . . . . . . . . . . . . . . 9 2.5 Signaling pathway example . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 2.6 A diagram of a gene regulatory network . . . . . . . . . . . . . . . . . . . . 11 3.1 The MAPK pathway scheme . . . . . . . . . . . . . . . . . . . . . . . . . . 15 3.2 The S-System equations for the MAPK pathway . . . . . . . . . . . . . . . 16 3.3 The components of a Hybrid Functional Petri Net . . . . . . . . . . . . . . 17 3.4 CHARON structure of the bacterium Vibrio fischeri . . . . . . . . . . . . . 19 3.5 An example of an enzyme-catalyzed reaction described using π-calculus . . 20 3.6 Obtaining the automaton from the simulation profile . . . . . . . . . . . . . 21 3.7 A simple regulatory network and its PADE representation . . . . . . . . . . 23 3.8 The phase space box for the simple network . . . . . . . . . . . . . . . . . . 24 3.9 Transition system derived from the differential equations . . . . . . . . . . . 24 3.10 A Bayesian network . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 3.11 A Boolean network . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 3.12 Two possible trajectories of a Boolean network . . . . . . . . . . . . . . . . 27 3.13 Inheritance hierarchy of the SBML classes . . . . . . . . . . . . . . . . . . . 29 3.14 The structure of the SBW framework . . . . . . . . . . . . . . . . . . . . . . 31 3.15 SBW Broker . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 3.16 Cell Illustrator . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 3.17 Concentration profiles of a simulation run . . . . . . . . . . . . . . . . . . . 33 4.1 The validation paradigm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 4.2 Issues to be addressed . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 5.1 The Hybrid Functional Petri Net model of the Wnt pathway . . . . . . . . 42 5.2 Profiles of β-catenin . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 5.3 Cell lineage of stem cells . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44 5.4 The scheme of the Wnt/Ca2+pathway . . . . . . . . . . . . . . . . . . . . . 44 5.5 The Hybrid Functional Petri Net model of the Akt and ERK pathways . . . 46

List of Figures iv

5.6 Directed acyclic graph of the Akt pathway . . . . . . . . . . . . . . . . . . . 49 5.7 Akt and PDK1 knockdown experiments . . . . . . . . . . . . . . . . . . . . 54 5.8 Knockdown experiments on the ERK pathway . . . . . . . . . . . . . . . . 55 5.9 Akt and PDK1 knockdown effects on P90RSK and Bad . . . . . . . . . . . 55 5.10 Effects of treatment on Akt . . . . . . . . . . . . . . . . . . . . . . . . . . . 56 5.11 Effects of treatments on MEK and ERK activity . . . . . . . . . . . . . . . 57 5.12 Profiles of MEK and ERK activation levels compared to experimental data 58 5.13 Profiles of protein activity under different treatment . . . . . . . . . . . . . 59 5.14 Profiles of model simulation compared to validation set . . . . . . . . . . . 60

LIST OF TABLES

5.1 Rate reactions of the Wnt pathway . . . . . . . . . . . . . . . . . . . . . . . 43 5.2 Initial concentration of the various reactants . . . . . . . . . . . . . . . . . . 43 5.3 Rate reactions and their parameters . . . . . . . . . . . . . . . . . . . . . . 53 5.4 Initial concentration of the components . . . . . . . . . . . . . . . . . . . . 54

• 1. INTRODUCTION

Despite the advances in biotechnology, researchers are still unable to fully grasp the com- plexities of cellular functions. Just as we cannot understand how an automobile works by studying its individual components such as the seat belt, tail light and engine, we cannot fully explain how a cell behaves by looking only at its gene sequences and function. Re- searchers have realized this and are now actively trying to study the cell as a complex system containing several interacting components. This approach of studying the cell through the interactions between its components is broadly termed as Systems Biology. An important sub-area that we wish to term as Computational Systems Biology emphasizes on the use of mathematical and computational modeling techniques. Computational Systems Biology has been attracting the interests of professionals from various disciplines - mathematicians, chemists, physicists and computer scientists. The approaches and methods vary from group to group. The Davidson Laboratory at the Cali- fornia Institute of Technology focuses on building cis-regulatory gene networks that control endomesoderm specification in the sea urchin embryo [BD02, DRO+02]. They then vi- sualize the network evolution using software such as BioTapestry [LDB05]. Some groups try to infer the entire gene network from microarray data using probabilistic graphical models or information-theoretic approaches such as Bayesian Networks and Relevance Net- works [Fri04, KIM03, Hus03, MNB+04]. From the control theorists’ and computer scien- tists’ perspective, understanding the control mechanisms as well as techniques for qualitative analyses of the cellular functions are of more interest than identifying protein or gene inter-

• actions. They use models of computation such as Boolean Networks [Rae02, SDZ02], Petri

Nets [GKV01, MDNM00, ZOS03] and Hybrid Automata [GT04a, MP03] to represent and reason about biological pathways. Yet others concentrate on modeling biological systems using ordinary differential equations [FBV00, Voi00]. However, with the current techniques, it is difficult to accurately model the entire cell in silico such that it possesses predictive capabilities. This difficulty lies not in the lack of knowledge on the genes present but arises from the numerous possible interactions between the genes, proteins and other components in the cell. These interactions are often studied in the context of the pathways which they play a role in. Studying such interactions is

• 1. Introduction

2

already a non-trivial issue. Cross interaction between the components of different pathways make this endeavor even more challenging.

1.1 Objectives

Even though much work has been done to elucidate pathway structures such as gene regulatory networks, we believe that the key to further enhance understanding of cellular functions is to fuse the different types of pathways together. For our research, we are interested mainly in signaling pathways and gene regulatory networks. To this end, the following objectives have been identified:

• a. To formulate a set of principles for interpreting information from laboratory experi-

ments and hypotheses so as to be able to synthesize the structure of the pathways. Equally important is the development of techniques to validate the model.

• b. To develop a framework, possibly by extending existing models, to represent biopath-
• ways. The framework should take into account discrete and continuous entities, and

allows for multiple levels of time-scale and granularity

• c. Provide visualization tools and techniques so that both models and their simulation

results can be presented to the biologists in a manner that they can understand and reason about.

1.2 Rationale

The rationale for our objectives is two-fold. With all the model building techniques that are currently available, only a few have attempted to formalize the process of validating the model [YMM+05]. Most validation techniques in pathway modeling and simulation involves visual comparison with limited sets of experimental data [FDG+03]. For a model to be useful and credible, some form of “level of confidence” should be assigned to a model,

• r even a class of them with respect to experimental data, instead of using terms such as

“fit the experimental data quite well”. Secondly, despite trying to model large scale networks representing molecular interactions, few have tried piecing together networks that represent the different aspects of cellular path- ways such as signaling pathways and transcriptional networks. A cell is a reactive system, continuously receiving inputs from the external environment. The signals are transduced in and out of the cell via signaling pathways and physiological changes in the cell occur due to activity of transcriptional networks in response to such signals. Hence there is a need for such integration, in order to further enhance our understanding of cell functions.

• 1. Introduction

3

1.3 Overview

The rest of this document is organized into the following chapters. Chapter 2 presents some basic biological concepts. Therein we look at the three main classifications of biological pathways that researchers are interested in (although we will be only be focusing on two

• f them). Chapter 3 presents a review of recent works in Computational Systems Biology

(thereafter we will refer to it simply by the term Systems Biology). Chapter 4 states out the problems that needs to be resolved. Chapter 5 introduces two of the pathways that we have been working on - the Wnt pathways and the Akt signaling pathway. Here we present some preliminary results from the modeling of those two pathways and show some of the issues involved in the modeling process itself. The concluding chapter then summarizes the issues and outlines the plans for the remaining part of this work in order to fulfil our research objectives.

• 2. BACKGROUND - CONCEPTS IN BIOLOGY

In this chapter, we will briefly cover the basic biological concepts and the classifications

• f the various forms of chemical reactions and biopathways. This chapter is necessary to

identify the way researchers view cellular systems. Modeling in biology occurs over a very wide spectrum, as shown in Figure 2.1.

Figure 2.1: Complexity spectrum of biological system.

At the lowest level are the atoms that will influence one another to determine the structure

• f proteins and other organic molecules. Higher up will be the molecular level, where gene

and amino acid sequences are being identified and studied. The interactions between the genes and the proteins will then form the functional level. It is at this level where researchers use biopathways to visually lay out such interactions to determine how they affect cell

• functions. Several such functional pathways will then form the cellular model. This could

be further abstracted at the phenotype level, where whole organs such as the heart can be modeled [Nob02]. Typical reaction rates across this spectrum can range from 10−12 to 103 seconds. Systems Biology, as we see it, roughly spans from the functional level to the cellular level. Of course, it will be further extended in the future to encompass biological systems at its phenotype level.

• 2. Background - Concepts in Biology

5

2.1 The Cell

The cell is the basic building block of life. Except for some smaller organisms such as retroviruses which consist of just an RNA (Ribonucleic acid) molecule enclosed in a protein casing, all living things are made up of cells. They can either be unicellular (mostly prokaryotes and protists) or multicellular (mostly eukaryotes). Other than red blood cells, every cell contains the entire genetic information or the genome for that organism in the form of DNA molecules. For eukaryotes, these genetic materials are kept in a membranous nucleus within the cell. The means by which a cell utilizes the information encoded in the genes is by proteins, which are synthesized as a result of interpreting the gene sequences. 2.1.1 The Central Dogma The basic tenet of cellular activity is the Central Dogma for Molecular Biology. The central dogma states that DNA encodes genetic information and passes it down by replica- tion to form another DNA molecule. The cell uses the information by making RNA, which encodes information to make proteins. RNA is created by a process call transcription (as it is a direct one-to-one complementary copying) while the RNA leads to the synthesis of proteins by translation (the amino acid sequences are determined by interpreting the RNA sequence against a genetic code which is conserved over all organisms). There are of course exceptions to the central dogma such as reverse transcription, but we need not delve deep into that.

Figure 2.2: The Central Dogma of Molecular Biology

• 2. Background - Concepts in Biology

6

2.1.2 Transcription and Translation During transcription, RNA is being synthesized by a molecule known as the RNA poly-

• merase. The DNA molecule contains coding regions (nucleotide sequences that will encode

for the protein) and non-coding regions, which may be regulatory or junk DNA. The RNA polymerase will bind to the promoter, which is a control sequence on the non-coding part of the DNA. This promoter is usually close to the coding region, hence are also known as cis- Regulatory elements. Thereafter, the polymerase will open up the double-helix and move along the coding region of the gene, adding nucleotides to the synthesized RNA molecule according to the sequence of the gene. This will carry on until it reaches a stop sequence

• n the DNA molecule. Some proteins, called transcription factors, regulate this process by

facilitating the binding of the polymerase to the promotors. Others, such as inhibitors or repressors, work by preventing the binding of polymerase. Most genes will encode for a type of RNA called messenger RNAs (mRNA). These mRNA,

• r transcripts, are then exported to the cytoplasm (for eukaryotes). A protein complex called

the ribosome will then translate the sequences in these mRNA to produce proteins, which will then fold into their respective structures to perform their functions. Some of these proteins are transcription factors. This is how genes regulate one another - by encoding for transcription factors (or repressors) which will either up-regulate (or down-regulate) the expression of other genes. 2.1.3 Chemical Reactions The main effectors of the information encoded in the genes are the protein molecules that are being synthesized. They perform all the main functions in the cell, from motility to producing energy and they perform most of those tasks by participating in biochemical re-

• actions. Other than translation and transcription, there are several other types of reactions.

An extensive set of chemical and molecular processes has been elucidated and classified by the Gene Ontology Consortium [Gen00]. Under the gene ontology classification, all the chemical reactions are classified under the category Molecular function. There are more than 100,000 different types of chemical reactions but for our purpose, we can view them as being the following types (following the reaction types in the Mammalian Interaction Map [Koh99]).

• a. Non covalent binding and dissociation - Association reaction between two or

more proteins. This binding will form dimers or multimolecular complexes. The re- verse of this reaction is dissociation. Usually in biological systems, these two reactions

• ccur simultaneously. Stability is achieved when the binding reactions balance out

the dissociation reactions, leaving no net change in the protein concentration

• 2. Background - Concepts in Biology

7

• b. Covalent modification - This is usually an enzyme-catalyzed reaction. An example

is phosphorylation, where a phosphate group is covalently bonded to an amino acid residue in a protein molecule. Proteins such as transcription factors are only active when they are phosphorylated

• c. Stoichiometric conversion - This reaction is similar to covalent modification, except

that the product is often different from the substrate on a much larger scale

• d. Synthesis - In certain models, it is infeasible to mathematically define all the tran-

scription and translation reactions that synthesize the proteins. Hence they are usually abstracted into single synthesis reactions, which usually have constant rates

• e. Degradation - All molecules have a certain life-span and will eventually degrade,

either naturally or by some other degradation machinery such as proteosome

• f. Translocation - This is not exactly a biochemical reaction per se. However, some

proteins are required to be present in certain compartments before they can affect cellular functions, such as the need for transcription factors to be physically present in the nucleus. Hence the movement of molecules from one part of the cell to another, either by diffusion or with the aid of transport proteins, forms an important type of biochemical reaction itself The reason why these reactions are classified in such a way is due to the way their reactions are modeled mathematically. For example, in enzyme-catalyzed reactions such as most covalent modifications, the limiting factor is the concentration of the enzymes. No matter how high concentration a substrate is, the reaction it is involved in can never go faster than its maximal velocity, given a fixed concentration of the enzyme that catalyzes it.

2.2 Biological Pathways

At the heart of the cell circuitry is the biological pathways. Unlike the circuits we see in electronics, these pathways are what constitute the functional view of the cell. A cell oper- ates by the above mentioned chemical reactions, with its molecules continuously interacting with one another. If we were to denote each type of reaction by an edge in a graph and all the molecular types and genes as the nodes, a network-like structure will be formed, i.e. the biological pathway. However, there are several types of chemical reactions involving different types of molecules, i.e. proteins, mRNAs, DNA and ions, just to name a few. We can roughly classify all forms

• 2. Background - Concepts in Biology

8

• f biological pathways into the three main types [TYH+02], each concentrating on a par-

ticular aspect of cell function:

• a. Metabolic pathways
• b. Signaling pathways
• c. Gene regulatory networks

Figure 2.3: The types of biological pathways in the cell

2.2.1 Metabolic Pathway In molecular biology, metabolism is the chemical process that can be dichotomized into the synthesis of organic molecules (anabolism) and their breakdown (catabolism). These sets of chemical reactions are usually concerned with the maintenance and housekeeping

• perations in the living cell. Catabolic pathways depict how foodstuffs are being broken

down into smaller molecules. One such pathway is Glycolysis (Figure 2.4), where glucose is being broken down into pyruvate via a multistage process. In this process, energy is released in the form of ATP (Adenosine triphosphate) which will be used to drive other energy-requiring reactions. Anabolic pathways, on the other hand, shows how ATP is used to drive the synthesis of other molecules, such as the amino acids that will later be translated into proteins. All the reaction steps together will then form the metabolic pathway. There are several publicly available databases that contain information on these pathways [KEG] as well as charts depicting several metabolic pathways such as the Roche Applied Science Biochemical Pathway Chart and Nicholson Metabolic Pathway Chart [Mic, Nic]. However as we shall

• 2. Background - Concepts in Biology

9 Figure 2.4: The glycolysis metabolic pathway. The individual nodes are the molecule types. Arrows depict chemical reaction. They are labelled with the enzymes that catalyze them. For some of the reactions, ADP is consumed and ATP is produced. Overall, two ATP molecules are being released in this reaction scheme

see later, just knowing all the interactions is not sufficient to understand the dynamics of the cell. 2.2.2 Signaling Pathway A cell is a complex reactive system. Unlike metabolic pathways which show how the cell creates energy and organic molecules, signaling pathways depict the chemical reactions that occur in response to external signals. External signals exist in the form of molecules

• r ions, collectively termed ligands. The ligands will latch onto membrane-bound proteins
• called receptors. The receptors are usually transmembrane, the other half being in the
• cytosol. The binding of ligands to the receptors will then cause the receptors to undergo

conformational changes which will either activate proteins in the cytosol, or act as scaffolds so that internal proteins can latch on and interact with one another. The cell transmits signals through a second messenger system. Using this system, signals are sent into the nucleus not via the signaling molecules but through a series of protein modifications - noticeably phosphorylation, in which a phosphate group is covalently bonded to an amino acid of a protein. The reverse of this reaction is dephosphorylation where the phosphate group is removed. Some of the modified proteins can then translocate into the nucleus where it may up-regulate or down-regulate gene expression. Other than this mode

• f signal transduction, small diffusible molecules such as calcium ions can also be used to

transmit the signals. This entire network of events would then form what is known as

• 2. Background - Concepts in Biology

10

signaling pathways or signal transduction pathways. Several pathways are already known and are being actively studied, such as the Wnt canonical pathways, the Akt pathway and the MAPK cascade. Although the pathways are usually being studied individually, it is known that they do not always operate in isolation. Almost every one is amenable to cross interaction with one another [AL04, HKN+03, NP00] and it is necessary to study such interactions to fully understand how the cell handles external events. An example of a diagrammatic view of the signaling pathway is shown in Figure 2.5.

Figure 2.5: An diagram of a signaling pathway - the Akt signaling pathway

2.2.3 Gene Regulatory Network A gene regulatory network is a collection of genes in a cell which interact with one another, governing the rate of transcription. It is usually represented as a graph, where the genes would form the nodes, and the edges will represent their interactions. The edges do not imply direct interaction. Instead, it implicitly means that the gene at the source end of the edge will most likely encode for a transcription factor that, upon activation, will regulate the expression of the gene at the destination end. On diagram, a gene regulatory network may look like a graph (Figure 2.6). The gene regulatory network has generated much interest because the cells in an organism differ only in the states and the expression activity of their genes. For example in the nematode Caenorhabditis elegans, its vulval cell types differ in expression, with vulval types vulE and vulF expressing the gene egl-17 during the late third larval stage while the other vulval cell types do not express this gene. To further complicate this issue, even within the same cell lineage, gene expression differs in a temporal manner, with vulE and vulF shutting

• 2. Background - Concepts in Biology

11 Figure 2.6: A diagram of a gene regulatory network [SI04]. The nodes would represent the genes and the edges represent the influences genes have on one another

• ff their egl-17 expression in the fourth larval stage [IWR+05]. Although in a gene regulatory

network, such spatial and temporal features are not shown, but it serves as a foundation for modeling and further understanding of such a phenomena, i.e. Epigenetics [PP05].

2.3 Summary

Generally a cell can be studied from several points of view. Advances in Systems Biology, however, are more concerned with the functional aspects of the cell, i.e. the biological

• pathways. It has been recognized that there are three major classifications of pathways.

The metabolic pathways concentrate on chemical reactions and fluxes that maintain the internal state of the cell. They perform the ‘housekeeping’ functions of the cell and can

• ccur throughout the cell. Signaling pathways form the ‘data transfer’ mechanism of the

cell, guiding how signals, in the form of ligands, trigger off a series of activities that gets transmitted into the nucleus. Some signaling pathways also involve cell-cell interactions, such as the Delta-Notch signaling pathway [MDMK00]. Finally, gene regulatory networks form the ‘processor’ of the cell, receiving proteins as input and producing mRNAs as output. This view is a simplification of how the cell functions. There are other exceptions such as the presence of DNA molecules in the mitochondria, an organelle that can be found in the cell. Even though the aim of Systems Biology is to study all the parts together, some difficulties remain as to how best we can integrate the three different kinds of pathways together. As the pathways themselves are studied separately, there are several ways to model each of

• them. For this research, we will only be focusing on signaling pathways and gene regulatory
• networks. Even so, differences in the operation semantics of the models representing them

must be reconciled before the cell can be studied as a whole.

• 3. A REVIEW OF CURRENT TECHNIQUES

Systems Biology has attracted the attention of people of various fields. As such, there are various approaches to modeling biological systems. Previously, kinetic modeling was the main method for computationally simulating and analyzing these pathways. But in the recent years, several other models of computation and analysis techniques have been used to model and study biopathways. In this chapter, we shall review some of these methods.

3.1 Quantitative or Qualitative Modeling?

Biology is starting to change from a qualitative, descriptive science to a quantitative, predictive science [WBI99]. Quantitative methods use numerical representation, such as protein concentrations and the rates of chemical reaction, to describe the system. The models used in quantitative analysis can be simulated to trace the changes in concentration levels of the components in the system. Qualitative methods, on the other hand, use non- numerical examination or interpretation of observations to discover relationships between the components of the pathways. Although informal, the diagrams used by biologists are examples of models that are qualitative in nature. They depict only the interactions between the proteins and not their concentration levels and associated dynamics. Despite huge advances in biotechnology, most of the data that can be collected from experiments are qualitative in nature. For example, Western blot images can only show the absence or varying degrees of presence of the proteins of interest. Although there are methods to quantify these data, such as densitometric analysis, their accuracy are still subjected to various source of noise. As such, in order to be able to construct accurate quantitative models, Systems Biology will require further technological breakthroughs. Although quantitative, predictive models are desirable in modeling biopathways, the im- portance of qualitative models should not be undermined as they allow researchers to dis- cover the underlying principles of cell functions, such as the phenomena of having kinase cascades to achieve signal amplification [HJ96]. Interestingly, there are some models of computation that have both quantitative and qualitative features. These models are known as Hybrid Models [LT04]. At present, there is no best way to model, simulate and analyze biological systems, but the use of hybrid models do seem to be an attractive method.

• 3. A Review of Current Techniques

13

The remaining part of this chapter will briefly introduce some of the models that are being used in Systems Biology. We start off first with the basic kinetic modeling using

• rdinary differential equations.

3.2 Kinetic Modeling

Before the term “Systems Biology” has been coined, kinetic modeling [RS04b] has been the predominant method in representing and simulating biopathways. The attributes in such models are the concentration of the molecules that make up the system. A basic assumption is that the molecules are uniformly distributed within their compartments (e.g. cytoplasm), such that the system can be represented by a vector X where each xi ∈ X represents the concentration of a molecule type. As the molecules in the system react with one another, their concentration will change with time. For signaling and metabolic pathways, the reaction rates are denoted by the mass action law, which states that the speed of a chemical reaction v is proportional to the concentration of its reactants. Suppose we have n molecular species and m chemical reactions, the rate vj of reaction j is specified by the following equation vj = k

n

• i=1

xci

i

(3.1) where k is the rate constant and ci is the stoichiometric coefficient for the molecule xi. The stoichiometric coefficient refers to the number of molecules of xi that are involved in a single instance of reaction vj. Thus the rate of change for a molecular species xi will be represented by the mass balance equation dxi dt =

m

• j=1

cijvj (3.2) where cij is the stoichiometric coefficient for the reaction vj. Hence the system can be represented as a set of ordinary differential equations, each denoting the rate of change for a molecular species. Except for simple systems, there are no analytic solutions to these equations and they have to be integrated numerically to observe their dynamics. 3.2.1 Michaelis-Menten Kinetics Most biopathways have parts that are specified by enzyme-catalyzed reactions. Although the kinetics of such reactions can be modeled using the mass action law, some of the interme- diate molecules are unstable, making measurements and derivation of parameters difficult.

• 3. A Review of Current Techniques

14

The Michaelis-Menten model is a simplification of such kinetics in these pathways. It is derived from the mass action formulation with certain assumptions. One of the reasons why the Michaelis-Menten model is preferable over the mass action form is that the parameters needed are easier to obtain experimentally. The Michaelis-Menten equation is derived from the following chemical scheme for enzyme catalyzed reactions S + E

k1

⇋

k−1 S.E kcat

− → S + P where S is the substrate, E is the enzyme, P is the product, and S.E is the intermediate complex formed when the enzyme E binds to substrate S. k1 and k−1 are the association and dissociation rate constants while kcat is the catalytic constant. Based on the above scheme, it then makes the following assumptions:

• a. The concentration of substrate S is much larger than enzyme E
• b. The degradation of product P back to the complex S.E is negligible
• c. This scheme reaches steady-state almost instantaneously - quasi steady state assump-

tion

• d. There is no cooperativity. The binding of the substrate to one enzyme binding site

does not affect the activity of any adjacent sites With these assumptions, the rate of this catalyzed equation v can then be worked out to be v = Vmax[S] KM + [S] (3.3) where Vmax is the maximum possible rate, given the enzyme concentration (sometimes also expressed as kcat[E]) and KM = k−1+kcat

k1

is the Michaelis constant. KM is obtained from experiments as being the concentration of the substrate at which the rate of reaction is half the maximal rate Vmax. It also denotes the enzyme’s affinity to the substrate. The lower the KM, the higher the affinity and less substrate is needed to reach half the maximal rate. An example of the MAPK cascade scheme and the rate reactions are shown in Figure 3.1. Several pathways have been constructed based on the Michaelis-Menten model [BF00, Kho00, SHNI02]. There are also variants that take into account other effects such as al- losteric effects (by introducing the Hill coefficient) where the binding of some molecules to an enzyme changes its affinity to other molecules, competitive binding and reversible reactions [RS04b].

• 3. A Review of Current Techniques

15

Rate Reactions Differential Equations v1 =

kcat1 [Ras][MKKK] KM1 +[MKKK] d[MKKK] dt

= v2 − v1 v2 =

Vmax2 [MKKK−P ] KM2 +[MKKK−P ] d[MKKK−P ] dt

= v1 − v2 v3 =

kcat3 [MKKK−P ][MKK] KM3 +[MKK] d[MKK] dt

= v5 − v3 v4 =

kcat4 [MKKK−P ][MKK−P ] KM4 +[MKK−P ] d[MKK−P ] dt

= v3 + v6 − v4 − v5 v5 =

Vmax5 [MKK−P ] KM5 +[MKK−P ] d[MKK−P P ] dt

= v4 − v6 v6 =

Vmax6 [MKK−P P ] KM6 +[MKK−P P ] d[MAP K] dt

= v9 − v7 v7 =

kcat7 [MKK−P P ][MAP K] KM7 +[MAP K] d[MAP K−P ] dt

= v7 + v10 − v8 − v9 v8 =

kcat8 [MKK−P P ][MAP K−P ] KM8 +[MAP K−P ] d[MAP K−P P ] dt

= v8 − v10 v9 =

Vmax9 [MAP K−P ] KM9 +[MAP K−P ]

v10 =

Vmax10 [MAP K−P P ] KM10 +[MAP K−P P ]

Figure 3.1: The classic MAPK scheme and the equations that describe the kinetics

3.2.2 Canonical Modeling and S-Systems The canonical approach, also known as Biochemical Systems Theory, to model biological systems was developed over thirty years ago [Sav69a, Sav69b, Sav70]. The basis of this approach lies in the premise that mathematical representation of nonlinear systems is based

• n the power laws. The S-System is one such example of a methodology based on canonical

modeling [Voi00]. Its basic form is a highly structured set of differential equations, where each equation denoting the rate of change for a molecular type xi is represented as the difference between the production reactions and consumption reactions. dxi dt = αi

n

• j=1

xgij

j

− βi

n

• j=1

xhij

j

(3.4) In this equation, αi and βi are the rate constants, gij and hij are the kinetic orders (which may not necessarily be an integer) and x1 . . . xn are the variables (dependent and indepen- dent). Using the MAPK scheme in Figure 3.1 as an example, the S-System equivalent is shown in Figure 3.2

• 3. A Review of Current Techniques

16

Using such an equation to represent all the molecular types, the entire system can be modeled in a systematic way. The advantages of using the S-System formalism is that we do not require the quasi-steady state assumption and the assumption that the concentration

• f the substrate is more than that of the enzyme.

It also takes into account allosteric relationships with the inclusion of kinetic orders. It has also been shown that the Michaelis- Menten model may not be totally compatible with biological systems [SS92]. On the other hand, one of the reasons why this formalism is not as widely used as compared to its Michaelis-Menten counterpart is that its parameters are not intuitive. The parameters gij and hij represent non-integer kinetic orders, not stoichiometric coefficients that are so common to biochemists. Also, it fails to describe some biochemical effects such as saturation and sigmoidicity [ESM04]. These two models, Michaelis-Menton model and S-Systems will form the basis for most quantitative models that we will come across later in this chapter.

3.3 Petri Nets

The Petri Net model is an attractive formalism for the modeling of biological net- works [PWM03]. It is a graph-based model, making it visually intuitive for biologists when reasoning about the pathways [MTA+03]. At the same time, it is based on sound mathematical foundation, having its own operational semantics. Thus it can be executed to simulate the dynamics of cellular pathways. The basic Petri Net is a bipartite digraph, consisting of two types of vertices - Places and

• Transitions. Places, usually representing passive entities such as molecules or experimental

conditions, are drawn as circles. Transitions, on the other hand, represent active components like chemical reactions or operations, and are depicted as rectangles. The state of the system is then represented by a vector, consisting of the markings, or tokens, in each place. A transition is enabled if all its incoming places contain a specified number of tokens, and it may then ‘fire’, consuming the tokens in its incoming places and creating new tokens

Differential Equations

d[x2] dt

= α2[x3]g23 − β2[x1]h21[x2]h22

d[x3] dt

= α3[x1]g31[x2]g32 − β3[x3]h33

d[x4] dt

= α4[x5]g45 − β4[x3]h43[x4]h44

d[x5] dt

= α5[x3]g53[x4]g54[x6]g56 − β5[x3]h53[x5]h55

d[x6] dt

= α6[x3]g63[x5]g65 − β6[x6]h66

d[x7] dt

= α7[x8]g78 − β7[x6]h76[x7]h77

d[x8] dt

= α8[x6]g86[x7]g87[x9]g89 − β8[x6]h86[x8]h88

d[x9] dt

= α9[x6]g96[x8]g98 − β9[x9]h99

Figure 3.2: The S-System equations for the MAPK pathway. Note that the variables have been substituted such that x1 = Ras, x2 = MKKK, x3 = MKKK − P, x4 = MKK, x5 = MKK − P, x6 = MKK − PP, x7 = MAPK, x8 = MAPK − P and x9 = MAPK − PP

• 3. A Review of Current Techniques

17

in its outgoing places. The entire model is then executed by the systematic firing of its

• transitions. Model characteristics and properties such as reachability and liveness can then

be deduced by observing the execution patterns of the Petri Net model. There are several variants and extensions to the Petri Net model, such as (i) Hierarchical Petri Nets [MSS05], which allow several levels where a net in one level can be represented as a place or a transition in another net of a higher level. (ii) Hybrid Petri Nets [MDNM00, MTA+03], where real numbers, in addition to integers, can be used as markings for places. (iii) Stochastic Petri Nets [GP98, SPB01], in which transitions may have delays that are given by a probability distribution. (iv) Colored Petri Nets [Run04] which allow places to store tokens of different types. (v) Time Petri Nets [PZHK05], where places and transitions may be given specific time delays. Each of these variants have their own expressive powers for representing biological pathways. Amongst the variants, one of the more popular models is the Hybrid Functional Petri Net [MDNM00]. This model has discrete and continuous versions of the Petri Net com-

• ponents. Places can now represent discrete properties such as gene states, or continuous

attributes such as protein concentrations. Its transitions can be discrete - firing only after its preconditions have been satisfied for a certain time delay, or continuous - firing at a rate specified by a firing function. It also introduces two new types of arcs - the Inhibitory arc and the Test arc. The inhibitory arc performs the opposite function compared to that

• f a normal arc, preventing its connected transition from firing when its preconditions are
• satisfied. The test arc is essentially the same as a normal arc, except that no tokens or

markings are being consumed from the place it is connected from, when the transition is firing.

Figure 3.3: The components of a Hybrid Functional Petri Net

In this model, the continuous transitions can have their firing functions defined as differ- ential equations such as that given in Sections 3.2.1 and 3.2.2. Therefore, besides having a

• 3. A Review of Current Techniques

18

visual structure that is similar to the pathways biologists are used to, its dynamics are also biologically relevant from a biochemist’s point of view. There is already a software package that adopts this formalism - the Cell Illustrator, and several pathways have already been modeled, including the λ-phage genetic switch mechanism, the circadian rhythms in fruit flies and mice, and the Delta-Notch signaling pathway1. For our preliminary studies, we have also adopted this formalism in the modeling of the Akt pathway and the Wnt canonical pathway. Details on these two pathways are presented in Chapter 5.

3.4 Hybrid Automata

The hybrid automaton [Hen96] is a formal model for hybrid systems, which contains both discrete and continuous components. Cellular processes need not be seen as being totally

• continuous. Certain mechanisms such as the on-off states of genes, or the states of the

cell can be described in a discrete manner, with the cell behaving differently when certain protein concentration levels cross their threshold levels [LT04]. Like the Petri Net formalism mentioned in Section 3.3, the Hybrid Automata is a graph based model. However, each node or vertex represents a discrete state of the cell, satisfying certain invariants and performing discrete jumps to another state when the boundary con- ditions, or jump conditions are met, such as a protein concentration crossing a threshold

• level. Within each state, the continuous variables, usually representing the concentrations
• f the molecular species in the cell will have their dynamics being modeled using differential
• equations. There will also be initial, invariant and flow conditions that control the values
• f the variables within that state.

Using this formalism, both continuous (e.g. protein dynamics) and discrete (e.g. mode changes and the switching of the genes) aspects of cellular functions can be modeled. Since hybrid automata has a control graph, which essentially is a state transition system, model checking techniques [AMP+03] can be applied to discover properties such as liveness of the pathway and reachability [GT04a, GT04b]. Alur et al [ABI+01] further extends Hybrid Automata by combining it with agent based

• models. Instead of having the entire cellular system represented by a single hybrid automa-

ton, they define each cell, protein or gene type - depending on the level of abstraction, to be hierarchical agents, with each agent containing its own automaton, describing its behavior.

1 http://genome.ib.sci.yamaguchi-u.ac.jp/∼gon/

• 3. A Review of Current Techniques

19

The entire system can then be simulated as a population of concurrent agents, communicat- ing with one another to influence their behavior. For this purpose, they have developed a language CHARON [AGH+00], which incorporates ideas from concurrency theory to model such systems. Figure 3.4 shows the hybrid automata of the bacterium agent Vibrio fischeri.

fOF F

1

= η1( 1

2 c − x1)

fi

4 = η4(x2 − x4) − r4x4

fP OS

1

= η1

4 (3c + xv81

8

κv81

81 +xv81 8

− 4x1) fi

5 = η5(x2 − x5)

fNEG

1

= −η1x1 fi

6 = η5(x2 − x6)

fOF F

2

= −η2x2 fi

7 = −η7x7 + r4x4 − rmem(x7 − x9) − r37,Rx3X7

fP OS

2

= fNEG

2

= η2(

xv82

8

κv82

82 +xv82 8

− x2) fi

8 = −η8x8 + r37,Aix3x7

fi

3 = η3(x1 − x3) − r37,Aix3x7

fi

9 = −η7x9 + rmem(x7 − x9) + u

Figure 3.4: CHARON structure of the bacterium Vibrio fischeri [ABI+01]. The discrete modes are hierarchical, with two main modes OFF and POS-NEG. Within the POS-NEG mode, the agent can either be in the POS mode or the NEG mode. The dynamics of its variables xi are defined by various differential equations depending on which mode the agent is in

3.5 Process Calculi

Other than using quantitative models such as differential equations to represent chemical reactions in the cell, there have also been efforts to use qualitative models such as formal

• languages. The practise of biology has been dominated by the use of symbols, diagrams

and words to explain biological phenomena such as signaling and gene regulation. However, precise reasoning using these informal methods requires much effort and these models are not amenable to simulation to generate out all possible cell behaviors. There are researchers trying to formalize the symbols, leading to the formulation of diagrams such as the Molecular Interaction Maps [Koh99] to represent biochemical networks, defined in a way such that it

• 3. A Review of Current Techniques

20

System ::= Enzyme | Substrate Enzyme ::= (v rels relp). bindsrels, relp.EX(rels, relp) + bindprels, relp.EX(rels, relp) EX(releases, releasep) ::= releases.Enzyme + releasep.Enzyme Substrate ::= binds(erels, erelp).X(erels, erelp) Product ::= bindp(erels, erelp).X(erels, erelp) X(reles, relep) ::= reles().Substrate + relep().Product

Figure 3.5: An example of an enzyme-catalyzed reaction described using π-calculus

leaves no room for misinterpretation [CFT01, Kit03]. From the computer scientists, they themselves have also come up with a variety of languages with the aim of trying to find one that is rich enough to describe most of the processes in the cell, and at the same time, is rigorous enough to allow formal analysis [Car04, DK03, RPS+04, RSS01]. Process calculi is one such formalism. Process calculi is a family of algebraic representation that allows one to describe high- level descriptions of interaction, communication and synchronization in concurrent systems. Members include the Calculus of Communication System (CCS), Communicating Sequential Process (CSP), and more recently, the π-calculus. The premise of this family of languages is that molecules can be modeled as processes - the ‘molecules-as-computation’ abstrac- tion [RS04a]. These processes interact with one another by communicating via channels. Modifications to the molecules will then lead to state changes of the processes that they

• represent. One of the features unique to π-calculus is the mode of communication where

the channels themselves are being passed as messages such that communication capabilities

• f the processes are not specified a priori. This dynamic communication scheme is used

to model mobility. An example of a system being specified using process calculi is shown in Figure 3.5. Simulation of such models involve applying rewrite and reduction rules to the terms that make up the system until no more unique rules can be applied. The simu- lation trace can then be shown as a transition system and analyzed using model checking techniques. There are several extensions to this formalism, adding quantitative features such as time and probability [PRSS01]. Graphical representations have also been added [PC05]. How- ever, to model a system using process calculi requires an in-depth understanding of the syntax of the language and its reduction rules. Furthermore, by representing each molecule as an individual entity, the modeling of large systems will place a huge computational burden on the computer hardware. Hence it has not yet been widely adopted by biologists.

• 3. A Review of Current Techniques

21

3.6 XS-System

The XS-System [APUM02], or the S-System extended with Automata is a model that was created to automate the process of deriving a hybrid automaton (Section 3.4) from simulation traces using the S-System formalism (Section 3.2.2). Instead of analyzing a single simulation trace, the XS-System allows the analysis of multiple traces by combining them to form an automaton. The number of states are kept manageable by declaring equivalence relations and collapsing multiple states of the traces into a single node. The automaton is formed from the discretization of simulation profiles by discrete time steps, such that the collective values of its variables xi at each time t0 + hk will form a discrete state. h is the duration of a single time step and k denotes which step the simulation is in. The states are then collapsed using an equivalence relation Ri,δi such that within the collapsed state, this relation holds | ˙ xi(t0 + hk) − ˙ xi(t0 + (k − 1)h)| ≤ δi (3.5) where ˙ xi(t) is the derivative of xi at time t

Figure 3.6: Obtaining the automaton from the simulation profile. (A) The trace with each discrete time step being converted to a state in the automaton. (B) Collapsing the points on the trace into equivalent states based on an equivalence relation

After deriving the hybrid automata, the system can be queried using CTL (Computation Tree Logic) statements. For example, should the researcher wants to find out, out of the several traces, whether a particular protein p will exceed a certain value k, he would then be able to query the automata using this query eventually(p > k) by which the system should be able to reply a true, or false. Using this model, researchers will be able to reason

• 3. A Review of Current Techniques

22

about simulation traces using a precise language instead of having to visually compare them. However, the type of analysis is restricted mainly to answering queries.

3.7 Piecewise-affine differential equations

Piecewise-affine differential equations provide another means of analyzing biological sys-

• tems. Unlike the XS-System that attempts to build a transition system from simulation

traces, the piecewise-affine models generate out all possible reachable states of the system based on the relationships between parameters of their equations. Instead of using the kinetic models in Sections 3.2.1 and 3.2.2, this method uses a class

• f piecewise linear differential equations [MPO95] to model the dynamics of regulatory
• networks. For a system of n variables, the equations will have the form

˙ xi = fi(x) − gi(x)xi, xi ≥ 0, 1 ≤ i ≤ n (3.6) where x = (x1, . . . , xn) is a vector of protein concentrations, fi(x) is the function denoting rate of synthesis and gi(x)xi denotes the rate of degradation. The function fi : Rn

≥0 → R≥0

is then defined by fi(x) =

• l∈L

kilbil(x) (3.7) where kil is the synthesis rate parameter and bil is a regulation function and L is the possibly empty set of indices of regulation functions. bil describes the threshold conditions under which the protein encoded by gene i is being expressed. The conditions are essentially step functions s+, s− : R2 → {0, 1} such that s+(xj, θj) =

• 1,

xj > θj 0, xj < θj and s−(xj, θj) = 1 − s+(xj, θj) (3.8) The rate constant of degradation for protein xi is denoted by γi. The choice of using piecewise linear differential equations is three-fold. It is observed that the rate of activation

• f a gene often follows a steep sigmoidal curve [dJGH+04], much like a step function.

Also, this class of differential equations has mathematical properties that allows qualitative analysis of the system without having the need to numerically integrate the system. Lastly,

• ne does not need to know the exact values of the parameters in order to analyze the system.

What needs to be provided is just constraints of the form 0 < θ1

i < . . . < θpi i < maxi

where there are pi concentration thresholds for the gene i and maxi is a parameter that denotes the maximum concentration for the protein.

• 3. A Review of Current Techniques

23

˙ xa = κas−(xa, θ2

a)s−(xb, θb) − γaxa

˙ xb = κbs−(xa, θ1

a) − γbxb

0 < θ1

a < θ2 a < maxa

0 < θb < maxb

Figure 3.7: A simple example of a gene regulatory network and the equations that make up the

• system. Genes a and b express proteins xa and xb. There is a negative feedback such that both xa

represses the activity of genes a and b while xb suppresses gene a once they cross certain threshold

• values. The system of equations and the constraints on the threshold values are also shown above

Simulation of this model does not involve stepping through small and discrete time steps, numerically integrating the differential equations to derive the system dynamics. Instead, it uses an n-dimensional phase space box Ω = Ω1 × . . . × Ωn, which represents the entire state space. This state space is then partitioned into regular mode domains and singular mode domains [BRdJ+05] using the threshold values. This is because the protein encoded by a gene will be involved in different interactions as it crosses different threshold levels. The regions in the phase box where at least one of the variables is at one of its threshold values will form the singular mode domains. The phase space box and the domains for the network in Figure 3.7 is shown in Figure 3.8 A property of this set of differential equations is that within a regulatory domain D, all the trajectories will converge monotonically towards an equilibrium point given by x = κD/γD. Letting Φ(D) = {κD/γD}, if Φ(D)∩D = {} then all the trajectories will stay in D, otherwise they will leave D after some time. Using this property, de Jong et al [BRdJ+05] provided an algorithm to derive the transition system from a set of piecewise-affine differential equations with given constraints (Figure 3.9). The resulting transition system can then be passed into

• ther software for subsequent analysis such as reachability, stability, or oscillations.

Unlike kinetic models, this form of analysis technique is not plague by lack of informa- tion on the kinetic parameters [dJGH+04] (although some form of relationship between the parameters must be provided). Also, it does not require numerical integration of the

• 3. A Review of Current Techniques

24 Figure 3.8: The phase space box for the simple network in Figure 3.7. D denotes the domains in the phase space with D1, D3, D5, D11, D13 and D15 being the regular domains and the rest being singular mode domains Figure 3.9: The transition system derived from the differential equations. Some of the regulatory domains have to be partitioned, denoted by a number after the domain name, e.g. D2.1, D2.2. This is due to differing trajectories as they converge to a focal point within that domain

differential equations. Instead, it derives all the possible state space for the system and simulation would then be traversing through the resulting transition system - a step known as qualitative simulation. This method is much faster, and hence enable upscalability for more complex models. Currently, it is only applicable to gene regulatory networks due to its approach of approximating regulation by piecewise linear equations. Signaling pathways are more rich and complex in their interactions and hence are not suitable for such abstrac-

• tions. Nevertheless, the piecewise-affine differential equations allow one to analyze all the

possible dynamics of a gene regulatory network without having to numerically simulate all

• f them.

• 3. A Review of Current Techniques

25

3.8 Other Models

Other than the piecewise-affine differential equations (Section ??), all the models in the previous sections can be used to model metabolic pathways, signaling pathways and regu- latory networks. The following models, however, concentrates solely on the gene regulatory

• networks. The large-scale reconstruction of gene regulatory networks is another area of im-

mense research in Systems Biology. This is largely due to advent of microarray technology, which enables researchers to obtain high throughput data not only on a few sets of gene expression, but that of the entire genome. Most of the analysis techniques on such data involves statistical analysis [PSS+03]. However, researchers are increasingly interested in inferring gene network structures by reverse engineering techniques [dJ02, GF05]. There are two main strategies in such reverse engineering - the physical strategy and the influence strategy [GF05]. The models in this section use the influence strategy to reconstruct the networks. 3.8.1 Bayesian Networks A Bayesian network is a statistical formalism for representing informational or causal de- pendencies amongst variables of interests [Mur01]. In the case of gene regulatory networks, the variables would represent the genes of interest (and maybe some external conditions) and the edges will then denote regulation between the genes. A gene regulatory network modeled by Bayesian networks will form a directed acyclic graph G = (V, E). Letting each gene be represented by xi, 1 ≤ i ≤ n, the expression level of each gene xi would then fol- low a conditional probability distribution P(xi | parents(xi)) where parents(xi) denotes the genes that regulate xi. This relation is due to the Markov assumption, which states that the conditional probability distribution of a variable depends only on the values of its immediate parents. The Bayesian network is hence a graphical structure that encodes the following joint probability distribution P(x) =

n

• i=1

P(xi | parents(xi)) (3.9) An example of a Bayesian network and the probability density function that it encodes is shown in Figure 3.10. The main issue in this form of modeling is to infer the network structure and its parame- ters from a set of gene expression data. Model selection is done by evaluating the structure against cost functions, which includes the Bayesian Information Criteria and the Akaike In- formation Criteria [HTF01]. This problem is known to be NP-Hard and hence much effort has been put into designing heuristic-based algorithms for inferring an acceptable network structure [FLNP00].

• 3. A Review of Current Techniques

26

P(x) = P(x1)P(x5)P(x2 | x1)P(x6 | x5)P(x3 | x4)P(x4 | x1, x5)

Figure 3.10: A Bayesian network consisting of six genes and the joint probability distribution that it encodes

Bayesian networks provide a way to re-create the entire regulatory network consisting

• f hundreds and thousands of genes using just expression data. It is also based on sound

statistical foundation. However there are still several problems that need to be resolved. A major criticism for Bayesian networks is its lack of feedback loops. In a cell, gene regulation and signaling will have feedback loops in order to achieve homeostasis so as to maintain a stable condition. For example, in the Wnt canonical pathway, β-Catenin is known to up-regulate the expression of Axin2 [JZD+02], which inhibits the activity of β-Catenin, forming a negative feedback. Extensions to handle this issue includes the dynamic Bayesian networks [MM99, SI04] which infers the network using only time-series microarray data. Other problems include losing information due to the need to discretize the data, usually to just three values − −1, 0, 1 for downregulation, no effects and upregulation respectively. Kim et al [KIM03] alleviated this issue by providing a means to represent continuous data by a normal density function. Despite the limitations of Bayesian networks, it allows us to analyze the structure of a regulatory network with limited data and incomplete knowledge about the system. 3.8.2 Boolean Networks A Boolean network does not encode for probability. Graphically it resembles that of a Bayesian network. However, each variable is a boolean variable that denotes whether or not a gene is ‘on’ or ‘off’. The state of a gene xi at time t is hence, xi(t) = fB

i (x1(t − 1), . . . , xn(t − 1))

(3.10) where fB

i

is the boolean function for gene xi. Unlike the Bayesian networks, a Boolean network can be used to study and analyze the dynamics of gene expression. Simulation of this model occurs in discrete time-steps of arbitrary units by applying the boolean function to the inputs of all the genes. The state of the system is then the collective states of all the genes in the cell. The sequence of states obtained from simulations will form the trajectory

• 3. A Review of Current Techniques

27

• f the system. By varying the initial conditions, a set of possible state space can be obtained

by combining all the trajectories.

Figure 3.11: A Boolean network [SS96] with its functions displayed as a truth table. Within a row, the values in the unshaded regions represent the state of the regulators which will produce the state of the regulated gene, as shown in the shaded column Figure 3.12: Two possible trajectories with (A) Starting from the state xa = 1, xb = 0 and xc = 0 and (B) starting from xa = 1, xb = 1 and xc = 0. (A) results in a loop while (B) reaches a steady state where all the genes are off

Since the Boolean network is a finite, deterministic model, all trajectories will either end at a steady state or a state cycle. By analyzing this in the state space, interesting

• bservations can be made, including looking for genes that are necessary for certain cell

functions (by observing that they are ‘on’ in most trajectories). Inference of a boolean network from gene expression data can then be done from an information theoretic approach, as exemplified by the algorithm REVEAL [LFS98] or look- ing for the most parsimonious set whose expression changes are consistent with the out- put [AMK99]. Unlike the Bayesian networks, which can work quite well with incomplete data, the boolean network requires huge sets of perturbation experiments in order to infer the interaction functions [AMK00] hence limiting its practical use. 3.8.3 Association Networks Like the Bayesian networks, an Association network is also a statistical formalism that can be used for the reverse-engineering of gene regulatory networks. However, it does not

• 3. A Review of Current Techniques

28

employ bayesian statistics. Instead it uses similarity, or correlation such as the Pearson correlation, to infer gene interactions. An Association network does not form a directed

• graph. Instead, it forms an undirected graph where the edges indicate significant correlation

between its connected genes [MNB+04]. There are several algorithms for inferring association networks from a set of gene expres- sion data. How they differ is the way they prune their networks to remove redundant edges. One of the more efficient algorithms for such purpose is the ARACNE algorithm (Algorithm for the Reconstruction of Accurate Cellular Networks) [MNB+04]. Its uniqueness lies in it using the Data Processing Inequality [CT91] to remove the edges. The DPI states that if genes x1 and x3 interacts through a third gene x2 and no alternate paths exists between x1 and x3, then the least of the Mutual Information can come only as a result of indirect inter-

• action. Mutual Information is a quantity that measures the independence of two variables

x, y and is defined by I(x, y) =

• x,y

P(x, y) log2 P(x, y) f(x)g(y) (3.11) where P(x, y) is the joint probability density function of x and y and f(x), g(y) is the probability density functions of x and y respectively. With this property, the entire network is scanned for triplets and within each triplet, the pair with the least MI score will have their edge removed. Margolin et al tested the ARACNE algorithm on a synthetic network and has claimed that it performs better than Bayesian networks in reconstructing the gene regulatory net- work [MNB+04]. However a separate paper by Hartemink observes the opposite [Hus03]. Nevertheless, these three formalisms and the algorithms associated with them do allow us to handle expression data in more ways than can be achieved by clustering methods.

3.9 Standards and Tools

Together with the availability of a variety of methods to model biological systems, several tools have also emerged, each of them based on one, or more of the above models. Here in this section, we will introduce a few tools and standards that have been developed over the past few years. 3.9.1 Systems Biology Markup Language The Systems Biology Markup Language2, or SBML for short, is an ongoing international effort to formulate a universal data exchange protocol based on the popular XML (Exten- sible Markup Language) language. In short, it is a description of biological models that is

2 http://sbml.org/index.psp

• 3. A Review of Current Techniques

29

framework and platform independent. This project was started in 2000 and since then, it has evolved into its second version (Level 2). The rational for developing such a language is to allow researchers to use complementary tools for different analysis, as long as they are SBML compliant. It also allows researchers to exchange models without requiring the need for them to obtain the software. However, it must be said that SBML does not attempt to be a universal representation of biological models, reason being Systems biology, as a new field, will see several changes in the way researchers think about and view biological

• entities. Instead, it aims to encompass the diverse approaches and methods while allowing

communication about the most basic and essential parts of a model [HFB+04]. Already, there are more than 85 tools that are SBML compliant. Despite having lots of details, the SBML has a very simple hierarchical model for cellular

• systems. To encompass multiple concepts, all classes are inherited from the same base class
• SBase. Its inheritance hierarchy is shown in Figure 3.13 in a UML like formalism.

Figure 3.13: Inheritance hierarchy of the major datatypes in SBML [FH03]. They are all inherited from the same base class

A model defined using SBML will then have a basic model class which would often aggregate any of the following classes:

• a. FunctionDefinition - This defines custom functions which might be used in pathway

simulations.

• b. UnitDefinition - This defines the units that is to be used for quantitative modeling,

such as nM or s−1

• 3. A Review of Current Techniques

30

• c. Compartment - A cell can be thought of as a system containing different enclosures,
• r compartments, such as the nucleus, cytoplasm or mitochondria. However it does

not have to necessarily correspond to the internal structures of a cell. It could even be used to denote different regions within the cytoplasm

• d. Species - This will represent the different molecular entities in the system, such as

protein or genes

• e. Parameter - This class is used to declare variables that are to be used in reactions,

such as the Michaelis constants KM

• f. Rule, Reaction and Event - The dynamics and simulation rules of the system will be

defined by these three classes Using these structures and components, cellular systems (currently mostly emphasized on quantitative models) can be represented and read by several tools. Currently, SBML Level 3 is being drafted and some proposals for extensions include classes to allow modularization, diagrammatic layout, multi-component species and arrays [HFB+04]. There are other XML based languages used to represent biological systems, such as CellML, BIOpolymer Markup Language (BIOML), Chemical Markup Language (CML), just to mention a few. However, none of them has the massive attention and attempts

• f compliance enjoyed by the SBML standard. In the future, SBML may indeed be the

standard language for biological systems with its constant development and evolution. 3.9.2 Systems Biology Workbench The SBML in the previous section attempts to formulate a universal way of digitally representing a model. On the software side, there are also attempts to create a universal platform with the standard functionality from which other programs can be built upon - the Systems Biology Workbench (SBW)3. If the SBML is viewed as an approach to allow researchers to share models, then the SBW would then be a way to allow them to share tools. SBW is a free, open source framework. It has a simple and lightweight message passing architecture (Broker-based architecture - Figure 3.15) and is language neutral. Program- mers can choose to build their tools in whichever language they are comfortable with (Java, C++, Python) and the different tools can work with one another by communicating via the SBW. For instance, one can build a model creation tool, and for model execution, he

3 http://sbw.kgi.edu/research/sbwIntro.htm

• 3. A Review of Current Techniques

31

can choose to use a simulation engine created by others. This saves a lot of time as certain tools such as simulation engines and plotting tools need not be recreated.

Figure 3.14: The structure of the SBW framework [HFS+01]. It consists of a basic Biological Modeling Framework, from which other modules can be added on. The SBW contains some basic modules collection. The other modules can be developed and added separately

The SBW consists of two layers. At the core is a skeletal software known as the Biological Modeling Framework (BMF). It provides the scaffold for which the rest of the components can be added as plug-ins (Figure 3.14). The entire SBW package is hence the BMF and a basic collection of such plug-ins. Developers can then add in their own application specific layers. The advantages of using such a framework for creating tools for Systems Biology include

• a. Modularity - The modules, or plug-ins are hidden from one another via a common set
• f Application Programmer Interface. Hence changes and modifications of the tools

will not have much impact on the entire SBW framework and its associated tools

• b. Reusability - Development time can be reduced when solutions to common problems

can be re-used, such as simulation engines and integrators

• c. Extensibility - Similar to the SBML, the designers of SBW recognize that this field will

see immense changes and improvements in the coming years and hence extensibility, allowing new modules to fit in, becomes a major concern when designing it 3.9.3 Cell Illustrator Originally called the Genomic Object Net [DNF+03, NDMM03], the Cell Illustrator is a commercial stand-alone software by Gene Network International4 to model and simulate

4 http://www.gene-networks.com/

• 3. A Review of Current Techniques

32 Figure 3.15: Messages from the various tools are routed to one another using the SBW Bro- ker [HFS+02]. Developers will not be aware of the brokers, being able to pass messages through a set of standard APIs

biological pathways. Its modeling formalism is based on the Petri Net methodology (Sec- tion 3.3). Although not entirely SBML compliant, this software uses its own format - CSML,

• f which SBML is a subset of, hence enabling it to read models encoded in SBML [NNK+04].

The Cell Illustrator has an easy-to-use and intuitive user interface (Figure 3.16). The components of the Petri Net - places and transitions, are accessible from the toolbar and the user can just drag and drop the components onto the drawing canvas. The dynamics for various chemical reactions can then be edited on the canvas in the form of mathematical formulae. The Cell Illustrator has all the tools from modeling to analysis. After building the Petri Net pathway, the user can simulate the system by specifying the initial conditions and the total amount of time to simulate, as well as the length of each time step for numerical integration. As the model is being simulated, graphs denoting states of the places or concentrations of the various molecules will be plotted in real time (Figure 3.17). In addition, the user can choose to perturb the system by changing the values of the places in the middle of a simulation run by pausing it. The results can also be saved as comma- delimited text files for visualization and analysis by other software. Accompanying this software is an animation program called the Cell Animator. Using this software, the user is able to create custom animations and graphs using the data generated from simulating the model. However, its animation features are still limited and it requires the user to manually key in all the image coordinates in an XML-like format and hence its usefulness is only restricted to simple presentations.

• 3. A Review of Current Techniques

33 Figure 3.16: Screen capture of the Cell Illustrator showing the modeling of biological systems. The above model is a simple model depicting two enzyme-catalyzed reactions. The enzyme for deactivation is not shown but the dynamics is modeled using Michaelis Menton kinetics as shown in the equations next to the transitions Figure 3.17: Concentration profiles of the proteins or molecules can be viewed in real-time as the model is being simulated. The figure shows the profiles of the substrate S and product P

Using the Cell Illustrator, large systems can be created easily. As mentioned in the pre- vious sections, Petri Net is an attractive formalism for it is a graph-based model and it resembles the signaling diagrams that biologists are used to. Hence the Cell Illustrator can be used by biologists without the need for rigorous knowledge on programming and

• mathematics. There is already a growing database of pathway models created by the de-

signers of Cell Illustrator5. The pathways include the λ-phage, the lac-operon, apoptosis and Delta-Notch signaling pathway.

5 http://genome.ib.sci.yamaguchi-u.ac.jp/∼gon/

• 3. A Review of Current Techniques

34

3.10 Summary

In this chapter, we have introduced some methods and tools in the area of Systems Biol-

• gy. However, the list is not exhaustive. There are several more modeling techniques and

tools that have the goal of creating a cell in silico. The kinetic models presented here are largely deterministic, with Michaelis-Menton and S-System forming the underlying dynam- ics driving the various models. There are also stochastic dynamics used to represent and simulate biological systems such as the Gillespie algorithm [Gil76, RA03] and the Chemical Master Equations [Gil92]. As more and more researchers from various fields recognize Systems Biology as a necessary step to further the understanding of cellular functions [IGH01], we will, in the future, see more creative methods and powerful software to aid them in uncovering the secrets of life itself.

• 4. PROBLEM FORMULATION

4.1 Framework for Cell Modeling

A cell can be viewed as a complex reactive system. As mentioned in Chapter 2, the heart

• f the cell circuitry are the biological pathways. In our research, we will only concentrate on

signaling pathways and gene regulatory networks. The major contribution of this research will be the specification of a framework for the modeling and simulation of the biopathways. The framework should be hybrid, being able to encompass both continuous and discrete components of the cell. It may also contain a heterogenous mix of existing models, each representing a different type of pathway.

4.2 Data Collection and Interpretation

To build a mathematical model, one requires sufficient data to develop hypotheses and infer the underlying models. The main concern with experimental data is whether or not the data is appropriate, accurate and sufficient [Sar03]. Biological data is inherently noisy. Although there are several high-throughput technologies such as microarray and flow cytom- etry, these data are indirect observations [WSW+05] of cellular activities and are subjected to several possible sources of noise. For example, in quantitative western blotting, there are several steps that need to be done before any readings can be taken, including

• a. Preparation of the cell lysate
• b. Extraction of the desired proteins by antibody immunoprecipitation
• c. Separation of the precipitate from the lysate by electrophoresis
• d. Extraction of the precipitate of the precipitate from the gel by transferring to a filter
• e. Blotting
• f. Imaging and densitometric analysis

Errors and inaccuracies can accumulate and propagate as the researcher follows through all the steps. What is observed at the end of the quantitation process is then, a quali- tative comparison of the concentration levels between various protein types. Similarly for

• 4. Problem Formulation

36

microarray data [HQA+00], the information that we obtain is the level of luminescence produced by the binding of expressed mRNA molecules to their complements on the DNA

• chip. The information that can be deduced from these measurements is whether or not

a gene has been up-regulated or down-regulated and its relative intensity of change. It does not give information such as the actual concentration of the mRNA molecules that are being synthesized. One point to note for microarray data is that the mRNA levels have to be amplified to obtain observable results. In short, an understanding of the experimental setup and initial conditions is mandatory so that systematic errors can be identified and taken into account when interpreting the data.

4.3 Validation

Model validation is “the substantiation that a computerized model within its domain of applicability possesses a satisfactory range of accuracy consistent with the intended appli- cation of the model” [Sch79]. It is an important issue which has largely been glossed over. Some researchers recognize model validation as an important factor [FHL+04, SHW04] but most do not specifically offer practical, systematic and reliable means of achieving it [CSK+03, HKW04]. Most validation in current research is based on subjective, visual comparison with experimental data such as Western blot images [HCC+05, PSH03]. To be able to create a model with a relatively high level of confidence in its correctness, there is a need to quantify the comparison of its structure and parameters with experimental

• data. Also, the usefulness of a model lies in its ability to predict the results of experiments

that have yet to be performed. On this, we have identified two topics for further exploration. 4.3.1 Metrics The purpose of metrics is to allow researchers to quantitatively evaluate the model struc- ture and parameters. Its importance lies not only in validation but also for the model building process such as parameter estimation, where the estimated parameters need to be constantly assessed for correctness [Kas05, KS04]. For reverse engineering of gene regula- tory networks, there are already some metrics that are being used. Mutual Information and Bayesian Information Criteria [HTF01] are some examples of metrics being used for assessing the correctness of the network structure based on gene expression data. For sig- naling and metabolic pathways, the more common metrics is the weighted root mean square distance [MMB03] between simulated and experimental data, which is often in the form of protein or metabolite concentrations. As mentioned previously, quantitative data obtained from experiments may contain a lot

• f noise and inaccuracy. Hence the weighted root mean square distance may not be suitable.

• 4. Problem Formulation

37

Other forms of metrics should be formulated to address this problem. It must also be noted that there cannot be one universal metric for assessing all the components and modules that make up the cellular system. This is due to the diverse nature and the way we choose to model and represent the various pathways. 4.3.2 Methods and Paradigms for Validation Validation of a model involves several steps. Sargent [Sar03] has formalized a paradigm (Figure 4.1) for model validation and presented a series of validation techniques. In model building, one or more of those validation techniques will be used. For the purpose of this research, we have adopted what he has defined as Historical Data Validation and Predictive Validation to build and assert the correctness of our models.

Figure 4.1: The validation paradigm as formulated by Sargent [Sar03]. Most of the validation efforts in this research will be on operational validation, where the simulated output is being compared with experimental ones

Under this paradigm, the historical data validation will use a part of the experimental data to build the model and estimate its parameters. The remaining data is then used to test whether the derived model will behave as expected. There are several methods that are based on this concept, such as the holdout technique [Koh95] and bootstrapping [HMM+02]. In Chapter 5, we will show how we use this method to validate the first iteration of the Akt pathway model.

• 4. Problem Formulation

38

4.4 Model Representation and Integration

Most signaling pathways are constructed using the Petri Net formalism or ordinary dif- ferential equations. On the other hand, gene regulatory networks are often represented by boolean or bayesian networks. To fully understand how the cell behaves, these pathways must be integrated together. This will allow biologists to study how the regulatory network respond to external signals, with the signaling pathways directing the signals and turning the genes ‘on’ or ‘off’. This integrated network can be viewed as having two layers which roughly mimics the actual compartments in the cell - the cytoplasm and the nucleus. In the cytoplasm layer, signaling pathways will be functionally separated into their respective modules, with their target proteins - usually a transcription factor - turning ‘on’ the genes which are in the nucleus layer. Interaction, however, will not be one way. There will be feedback loops from the regulatory networks back into the signaling pathways. Instead of formulating a new model that can encompass all the features of the three types of networks, it would be more practical to use existing hybrid models such as Hybrid Functional Petri Nets, whose continuous components can model the signaling pathways, while its discrete components can be used to represent the genes in the gene regulatory

• network. The models can then be extended if necessary. However, to successfully ‘hook

up’ these disparate pathways, communication issues such as time semantics and mapping

• f quantitative measures (e.g. the level of concentration an activated transcription factor

must reach before its target genes are turn ‘on’) have to be resolved.

Figure 4.2: Issues in biological modeling that are to be addressed in this research. Metabolic pathway is greyed out as it is not the main focus of our work

• 4. Problem Formulation

39

4.5 Issues to be addressed by proposed research

Given the issues mentioned above, the research proposes to study and contribute to the following areas

• a. Provide a means to represent data, especially western blot data, quantitatively or

semi-quantitatively, with appropriate assumptions

• b. Formulate metrics such that pathway and parameter correctness can be assessed
• c. Formalize the validation techniques for various types of pathways after model building
• d. Provide a hybrid framework in which different types of pathways can interact with
• ne another

• 5. PRELIMINARY STUDIES

This chapter reports our preliminary work on two widely studied pathways - the Wnt signaling pathways and the Akt pathway. Both pathways have important roles in the proper functioning of a cell, with the Wnt signaling pathway being involved in stem cell differentiation and the Akt pathway being one of the regulators of programmed cell death

• apoptosis.

In studying these two pathways, issues mentioned in Chapter 4 regarding model building process, data interpretation and validation have surfaced. In the sections that follow, we present our current techniques to deal with them.

5.1 The Wnt Pathways

The Wnt proteins are a large family of cysteine-rich1 secreted ligands that are involved in the development of cells. In addition, this class of proteins is also implicated in the genesis

• f human cancer.

The Wnt genes are defined by a sequence homology to the original members Wnt-1, first discovered in the mouse (then called the Int-1) and the wingless (Wg) in Drosophilia, hence the name Wnt. There are several members in this family of proteins (more than 16 mammalian family members) but not all their functions are known. These Wnt proteins are involved in more than one signaling pathways. The most well studied pathway is known as the Wnt canonical pathway, or the Wnt/β-catenin pathway. The central player in this pathway is β-catenin, which is a transcriptional co-factor with the T-cell factor/lymphoid enhancement factor (TCF/LEF) in the cell nucleus. Recently, other pathways that involve the Wnt proteins are being discovered, including those that control certain aspects of gastrulation move- ments, cochlear hair cell morphology, heart induction, dorsoventral patterning and tissue separation [VAM03]. These pathways are collectively known as the non-canonical pathways 5.1.1 The Wnt Canonical Pathway Wnt proteins are secreted proteins and they trigger off cell activities by first binding to a seven-transmembrane-spanning superfamily of cell-surface receptors, the Frizzled (Frz)

1 Cysteine is a naturally occurring hydrophilic(“water loving”) amino acid which is found in most proteins

• 5. Preliminary Studies

41

receptors. These receptors have an extracellular Wnt-binding domain and an intracel- lular C-terminal tail. Recent evidence also suggests that the co-receptors Low Density Lipoprotein receptor related protein (LRP5/6) are also involved in the Wnt canonical path- way [PBM+00]. Immediately downstream of this pathway is the Dishevelled (Dsh) protein. When the Wnt ligand binds to the Frz receptor in the extracellular domain, it recruits intracellular Dsh to the plasma membrane where it is activated via phosphorylation on the serine and threonine residues [WN98]. β-catenin is a cytosolic protein whose level is kept low due to phosphorylation-dependent degradation [He03]. It is phosphorylated by a degradation complex consisting of Glycogen Synthase Kinase (GSK-3β), Axin and Adenomatous Polyposis Coli (APC). In the presence

• f Wnt signaling, the activated Dsh prevents the degradation complex from phosphorylating

β-Catenin, hence enabling it to translocate into the nucleus and becoming the co-activator

• f several genes that are needed in cell growth and differentiation.

5.1.2 The Model For our preliminary study, we have constructed a HFPN model depicting the Wnt signal- ing pathway. This model is adapted and modified from [LSK+03]. A HFPN schematic of the model is shown in Figure 5.1. The kinetic equations, parameters and initial concentration levels of the various proteins are shown in Table 5.1 and Table 5.2. Using our model, simulations such as the effects of Wnt signaling on the levels of β-catenin in the nucleus can be performed (Figure 5.2). Using this model, other simulations such as the knocking down of the degradation complex components can also be done in silico. 5.1.3 Downstream Gene Regulatory Network Most studies on signaling pathways concentrate only on the effects of certain perturba- tions on downstream protein concentration levels [AK03, HKN+03, LSK+03]. For the Wnt canonical pathway, most simulation studies stop at the concentration level of free β-catenin in the cytoplasm. However, it is known that β-catenin is a transcription factor and can regulate gene expression, some of which are components of the Wnt pathway itself, forming a feedback control. Some of these genes have been elucidated and are shown in [Nus05]. 5.1.4 Stem Cells and Differentiation Stem cells, especially embryonic stem cells, have the potential to differentiate into other types of cells, such as mesodermal cells and endodermal cells. One of the challenges of stem cell research is to discover the control mechanisms that will influence stem cells toward certain cell types such that they can be used for transplantation to treat disorders such as

• 5. Preliminary Studies

42 Figure 5.1: The Hybrid Functional Petri Net model of the Wnt canonical pathway. (A) shows the binding of Wnt ligands to the Frz receptors. (B) depicts the formation of the degradation complex consisting of APC, Axin and GSK-3β. Phosphorylated Dsh will activate the transition to dissociate GSK-3β from the degradation complex. (C) is the degradation cycle where β-catenin gets phospho- rylated by the degradation complex, hence marking it for degradation by the proteosome (abstracted away by a degradation transition) and (D) shows the cytosolic β-catenin being translocated to the nucleus and binding to its co-factor TCF to start the transcription of its downstream targets Figure 5.2: Profiles of the concentration of β-catenin in the cytoplasm and nucleus in the absence and presence of Wnt signaling

Parkinson’s disease and stroke [LK05]. The Wnt canonical pathway is a crucial component in the differentiation of stem cells. It is required in the formation of the primitive streak as well as paraxial mesoderm development in the mouse [YTY+99], hence explaining the attention that it has been getting from the scientific community.

• 5. Preliminary Studies

43

No Rate Equation (nM.min−1) 1 0.1 × [Wnt] × [Frz] 2 0.1 × [Wnt.Frz] 3 0.18 × [Dsh] 4 0.018 × [Dshp] 5 0.01 × [Axin] × [APC] 6 0.5 × [Axin.APC] 7 8.2 × 10−5 8 0.17 × [Axin] ×

[APC] 50×[APC]

9 0.91 × [Axin.APC.GSK − 3β] 10 0.091 × [Axin.APC] × [GSK − 3β] 11 3 × [Axin.APC.GSK − 3β] × [Dshp] 12 0.27 × [Axin.APC.GSK − 3β] 13 0.13 × [Axinp.APCp.GSK − 3β] 14 1 × [βCatenin] × [Axinp.APCp.GSK − 3β] 15 120 × [Axinp.APCp.GSK − 3β.βCatenin] 16 0.42 17 2.6 × 10−4 × [βCatenin] 18 0.42 × [βCateninp] 19 210 × [Axinp.APCp.GSK − 3β.βCatenin] 20 210 × [Axinp.APCp.GSK − 3β.βCateninp] 21 0.1 × [TCF] × [βCatenin] 22 3 × [TCF.βCatenin] 23 0.1 × [βCatenin] × [APC] 24 120 × [βCatenin.APC]

Table 5.1: Rate reactions of the Wnt canonical pathway

Reactant Concentration (nM) Wnt 10 Frz 10 Axin 0.02 APC 100 βCatenin 34 βCateninp 1 Dsh 100 GSK−3β 50 TCF 15

Table 5.2: Initial concentration of the various reactants in the Wnt canonical pathway

To understand the importance of the Wnt signaling pathways in cell development, it is imperative to elucidate the mechanisms of how they affect the states of the underlying genes, or in some cases, how the gene regulatory network is being ‘reprogrammed’ [ALR03]. 5.1.5 Issues and Non canonical pathways Besides the canonical Wnt pathway, pathways involving other members of the Wnt family have been discovered. Examples include the Ca2+ and the planar polarity pathway [PP00, WM03], each playing its own role in the development of the organism. The schema of the Wnt/Ca2+ pathway is shown in Figure 5.4. Our future work will include modeling of the

• 5. Preliminary Studies

44 Figure 5.3: A simple cell lineage showing the possible ways stem cells can be differentiated. The Wnt pathway plays a role in guiding the stem cell towards specific lineages during different phases

• f development. Hence finding out how β-catenin affects the epigenetic state of the gene regulatory

network is an important and challenging issue

non-canonical pathways and integrating them together.

Figure 5.4: The scheme of the Wnt/Ca2+ pathway. Wnt5a, after binding to Rfz2, will cause the cytosolic G-Protein to break off into its α and β/γ components. The Gα subunit will activate phosphodiesterase (PDE), which will inhibit cyclic guanosine monophosphate (cGMP) activity in the cell. On the other hand, the Gβ/γ subunits will activate phospholipase C (PLC), which in turns hydrolyzes phosphatidylinositol 4,5-bisphosphate (PIP2) into inositol 1,4,5-trisphosphate(IP3) and diacylglycerol (DAG). DAG will activate protein kinase C (PKC) while IP3 increases the Ca2+ level in the cell, activating calcineurin and Ca2+-calmodulin - dependent protein kinase II (CamKII). Finally, calcineurin activates nuclear factor of activated T cells (NF-AT), influencing gene expression [WM03]

• 5. Preliminary Studies

45

5.2 AKT and ERK Pathways

The Akt pathway has been studied extensively due to its role in apoptosis and cancer. In the prostate cancer cell line LNCaP, genetic defects cause such cells to be lacking in the active lipid phosphatase - PTEN, which is an inhibitor of the Akt pathway [HLC+80]. Hence this cell line has unusually high amounts of activated Akt, promoting cell survival. In experiments to knockdown 3-Phosphoinositide-dependent Protein Kinase 1 (PDK1), one

• f the components of the Akt pathway, it was noticed that activity of other pathways are

also affected. One of them is the Mitogen-activated Protein Kinase (MAPK) pathway, also known as the Extracellular Signal-Regulated Kinase (ERK) pathway, which plays a part in cell growth and proliferation. We have created a model which has these two pathways as the main modules - the Akt and the ERK pathways. These two pathways are stimulated by external growth factors, which is modeled as a discrete element - serum. Several models representing these two pathways have already been created separately [HKN+03, PSA+04, PSH03]. In our model, we are interested not only in the individual pathways but also in their cross interaction and effects on common downstream targets. The mechanism behind the activation of the Akt pathway is the membrane recruitment and activation of the signaling proteins. Binding of growth factors to the cell surface recep- tors will cause them (the receptors) to act as scaffolds for specific binding interactions be- tween cytosolic proteins. This leads to the activation of Phosphoinositide 3-Kinase (PI3K), which will catalyze the phosphorylation of Phosphotidylinositol 4,5-biphosphate (PIP2) at the inositol ring, forming PIP3. This reaction can be reversed by the lipid phosphatase PTEN. Activated PIP3 will then recruit proteins containing the Pleckstrin homology (PH) do- mains, which includes the protein Akt. At the cell membrane, Akt will then be able to interact with PDK1. Through this interaction, Akt is activated by a sequential phosphory- lation of Thr308 and Ser473 [FHC05, NA00], the former being catalyzed by PDK1 and the latter by a yet-to-be-identified kinase, and is presently dubbed PDK2 [KNK+04]. Activated Akt will then be released into the cytoplasm where it will further activate downstream tar- gets such as the Bcl-2 Antagonist of Cell Death (Bad) and the Forkhead transcription factor (FKHR). Activated Akt is regulated by Protein Phosphatase 2A (PP2A), which deactivates it by removing the phosphate groups. Similar to the Akt pathway, the ERK pathway is stimulated by the binding of external growth factors to receptors on the cell membrane. This causes the proto-oncogene Ras to be

• 5. Preliminary Studies

46 Figure 5.5: The Hybrid Functional Petri Net model of the Akt and ERK pathways. (A) is the module showing the Akt pathway, starting from the activation of PI3K by Serum. (B) shows the ERK cascade from Ras down to ERK. (C)-(E) shows the other modules that are related to the two main pathways, with (C) representing the ROS module, (D) depicting the RSK module and finally (E), the Bcl-2 family module

converted to its active conformation, having Guanosine 5-triphosphate (GTP) bounded to it instead of Guanosine diphosphate (GDP). In this conformation, its binding affinity to the Raf kinase is increased, causing Raf to be recruited to the membrane and become activated. This, in turn, causes Raf to activate the dual specific protein MAPK/ERK Kinase (MEK) by phosphorylating it at its serine/threonine residues. The activated MEK then carries

• n to phosphorylate ERK [Kol00].

Finally, activated ERK will induce cell growth and proliferation by activating a variety of downstream transcription factors. The Akt pathway directly influences the ERK pathway in two ways. Phosphorylated Akt negatively regulates the ERK pathway by phosphorylating Raf at Ser259, thereby inhibiting its activity along the ERK pathway [RZS+01]. The second influence is hypothesized by

• bservations through wet-lab experiments. Knocking down of PDK1 expression levels lead

to decreased MEK and ERK activity, indicating that there is a positive regulation of PDK1

• 5. Preliminary Studies

47

• n the ERK pathway in a MEK dependent manner.

Other than these two cross interactions, there are also other modes of interaction. Ex- periments have shown that prostate cancer cells produce substantial amounts of Reactive Oxygen Species (ROS) such as Hydrogen Peroxide and Superoxide and one of the sources could be NAD(P)H oxidase (NOX5) [BCK+03]. This generation of ROS may be inhib- ited by the flavoprotein-dependent NAD(P)H oxidase inhibitor - diphenylene iodonium (DPI) [BCK+03]. ROS is known to mediate MEK and ERK activity, possibly via acti- vation of the p21-activated kinase 1 (PAK1) [DWVL98]. Activated PAK1 can then influ- ence the ERK pathway by regulating the phosphorylation of Raf [CKM+00]. Related to the PAK1 mediated ERK activation is the role of PI3K, adding on another possible cross interaction between the Akt and the ERK pathways. Experiments have shown that PI3K can regulate the activity of PAK1 by phosphorylation [CKM+00] and we have already seen previously that PAK1 regulates the ERK pathway, hence showing yet another indirect in- teraction between the Akt pathway and the ERK pathway. PAK1 is also shown to be able to phosphorylate Bad at Ser112 and Ser136 both in vitro and in vivo [SMS+00]. Bad is a protein that is involved in regulating apoptosis by the competitive heterodimer- ization of Bcl-2 and Bcl-XL with itself and Bax. Bcl-2 and Bcl-XL are anti-apoptotic proteins which inhibit the release of cytochrome c from the mitochondria [FBV00]. By releasing mitochondrial cytochrome c into the cytosol, procaspase 9, a protein directly in- volved in cell death, will be activated, leading to apoptosis. Phosphorylation of Bad at Ser136 by Akt and Ser112 by ERk via Ribosomal S6 Protein Kinase (P90RSK) will prevent it from forming dimers with Bcl-2 or Bcl-XL, enabling them to sequester with Bax, thereby preventing the release of cytochrome c, and hence encouraging cell-survival. The entire HFPN model, including the above-mentioned pathways and their cross inter- actions is shown in Figure 5.5. However, the model is still not completed as we do not know the rates of the various reactions that drive the model. This brings us to the issue of Parameter Estimation. 5.2.1 Parameter Estimation In any modeling endeavour, the bottleneck is the estimation of the parameters for the various rate reactions. Technical difficulties and huge resource requirements make the ex- perimental determination of all the parameters impossible and hence they have to be esti-

• mated. This amounts to a nonlinear global optimization problem and an in-depth review
• f the various approaches has been reported in [MMB03]. Several parameter estimation

algorithms have been considered and one of the best performing algorithm is a variant of

• 5. Preliminary Studies

48

the Evolutionary Strategies approach [RY00]. Every algorithm to estimate the parameters require searching through a large solution space to find suitable parameter values. Our proposed method is to exploit the structure of the network whenever possible, such that instead of searching the entire set of parameters, we concentrate only on a few parameters at any given time. This form of optimization is also known as Alternating Optimization [BH02]. The motivation for such a method is not only to reduce the size of the search space, but also the nature of the Akt pathway allows such grouping. Although not analyzed in depth, there have been instances where Alternating Optimization out-performs other techniques [BH02], having faster convergence

• time. However, there must be a basis for deciding which parameters are to be grouped

together and their order of estimation. To do so, there is a need to look at the pathway topology. 5.2.2 Pathway Topology Looking at the structure of the network in Figure 5.5, we can see that most of the reactions are catalytic reactions. Since we are adopting the Michaelis-Menten kinetics for enzymatic reactions, such reactions do not consume the enzyme, i.e. the concentration of the enzymes do not change with respect to the reactions that they catalyze. The concentration flow

• f the molecules in Figure 5.5 can be viewed as being contained within separate modules,

e.g. PI3K activity is contained within its own loop consisting of the places PI3K and PI3Ka and molecules do not flow into other places. Hence the places PI3K, PI3Ka and their associated transitions (excluding those that PI3K and PI3Ka are the enzymes of) will constitute a module. Dependency relation between the modules are determined by the catalytic influence of the places in a module on the transitions of another module. We let X be the set of variables (or places) representing protein concentration and V be the set of rate reactions (or transitions). Each variable is denoted by xi and its rate of change will then follow Equation 3.2 where the rate reaction is denoted by vj and cij is the stoichiometric coefficient of xi in reaction vj. Definition 1: (Module) A module is a smallest subset M′ = (X′, V ′) such that it satisfies the following properties. For each variable xi where xi ∈ X′ and its rate of change is specified by dxi/dt = |V |

j=1 cijvj, it must be the case that if cij = 0 then vj ∈ V ′. Also for

any two modules M′ = (X′, V ′), M′′ = (X′′, V ′′), X′ ∩ X′′ = Ø and V ′ ∩ V ′′ = Ø Hence we can simplify the entire network into a directed graph G = (V, E) where the nodes Vi ∈ V are the individual modules and the edges Ei ∈ E are the catalytic activa-

• 5. Preliminary Studies

49

tion/inhibition between the modules. An example of how the scheme in Figure 5.5 can be simplified into a graph is shown in Figure 5.6. In biopathways, some molecules regulate their activity by having feedback loops, which may be direct or indirect (by regulating the activity of upstream molecules). For the graph

• f such pathways, their cyclic subgraphs have to be fused together to form a larger module

such that the resultant graph forms a Directed Acyclic Graph. Having a DAG will allow us to specify the order for the estimation of the modules. However, this limits the way the pathway structure can be exploited and in the worse scenario, the modules are so large that there are no advantages in partitioning it. For our model of the Akt pathway, there are no such loops and hence the modules can be kept relatively small.

Figure 5.6: Directed acyclic graph of the individual modules that make up the pathway. For brevity, enzymes that do not change in concentration such as the PP2A and PDK1 are not shown in the diagram

5.2.3 Topological Ordering Estimating the parameters of the pathway involves fitting the model with available exper- imental data, usually steady state protein concentration levels, or less frequently, time-series

• profiles. Most parameter estimation algorithms regard the system in its entirety. In our

model, we shall focus only on optimization via Evolutionary Strategies. Such algorithms require several model evaluations to assess the fitness of the parameters or solutions. In Evolutionary Strategies, current solutions are being ‘recombined’ and ‘mutated’ in every iteration. Hence for each iteration, or generation, the algorithm has to work with a

• 5. Preliminary Studies

50

solution space with exponential dimensionality. Whilst this is inevitable, we can reduce it by fitting the separate modules individually. Taking the pathway in Figure 5.6 as an example, one can see that simulation of the PI3K module does not depend on the Pak1 module. However the same cannot be said for the reverse. Hence the PI3K module must be fitted prior to Pak1, leading to the notion of precedence constraints. The nodes in the DAG can therefore be topologically

• rdered and their parameters will then be estimated in that order.

There are already several efficient algorithms to perform topological sort on a DAG and hence it will not be further elaborated [Ski97]. 5.2.4 Rank Despite being arranged in topological order, it is sometimes not possible to estimate the parameters due to unavailable data for fitting them with. As such, there is a need to group multiple modules together and assign them a rank to denote which modules are to be fitted together. Usually for fitting, there will not be one, but a series of experiments, each providing information on certain variables. For a variable xi, we denote (if available) its profile for experiment j as x(e)

ij .

Definition 2: (Rank) The rank is a grouping of the modules in their linear extension such that for each group {M′} of the same rank, there exists exactly one M′

k that has at least a

variable xi with an experimental profile x(e)

ij

Additional requirements are that other than the last rank, the module with the experimental profile must be the last module within its group when arranged in topological order. When representing the network as a graph G = (V, E), where each vertex Vk represents the module Mk, vertices with no experimental profiles take on the rank of the closest vertex Vl where there is a path Vk...Vl. This will impose a slight restriction on the partial order of the modules that do not have experimental profiles to fit them to. The requirement for each rank to have only one module with a profile is to minimize the number of modules in a rank, so that the solution space to search from is also kept to a minimal. For the graph in Figure 5.6, we have the experimental profiles for the following protein types - AKT, Raf, MEK, ERK, P90RSK and Bad. Hence the modules and their grouping will be in this order - {PI3K, AKT}, {ROS, Pak1, Ras, Raf }, {MEK}, {ERK}, {P90RSK}, {Bad}, with the highest rank of 1 being the group {PI3K, AKT} and the lowest rank of 6 for the group {Bad}.

• 5. Preliminary Studies

51

5.2.5 Evolutionary Strategy The choice of estimation algorithm being used is independent on the ranking, i.e. one can choose to perform deterministic optimization algorithms such as the Levenberg - Marquardt method, or stochastic methods such as Genetic Algorithms, or even a mix of them for different ranks, depending on the characteristics of the modules within that rank. We have chosen Evolutionary Strategies for our work, specifically the (µ+λ)−ES algorithm [BS02]. The following steps show the basic implementation of the algorithm.

• 1. Generate µ parent vectors, each being a set of parameters Pi = (pi1, . . . , pin).
• 2. Create λ new offsprings, with each child being formed from the recombination of

randomly selected parent parameters.

• 3. Mutate the offsprings.
• 4. For each child parameter set, evaluate the model to assess its fitness score.
• 5. Select the µ most fit parameters from the (µ + λ) parameters to form the next gener-

ation.

• 6. Repeat steps 2 to 5 for a predetermined number of generations, or when no better

parameters could be obtained. This process is repeated in rank order. In addition, we need to maintain a set of fixed parameters obtained after estimating parameters of the previous ranks or those that have already been determined beforehand. We denote this set as Pf and for each generated parent parameters Pi in step (1), Pf ⊆ Pi. Recombination and mutation are then performed on the parameters pi where pi / ∈ Pf. To assess the fitness of the parameters, the cost function must also take the current rank k into consideration. J =

• i,j,t

wij

• (xij(t) − x(e)

ij (t))2

(5.1) where xij(t) is the predicted value of xi in experiment j at time t. x(e)

ij (t) is its corresponding

experimental value. wij is the weight used to normalize the contributions of each term to the cost function and it is usually taken to be wij = 1/max(xij). Also, for all xij used in the cost function, Rank(xij) ≤ k. After estimating the parameters for a particular rank k, we append the best parameters with rank k to the set of fixed parameters Pf. Based on this algorithm, we estimate the parameters for the rate reactions that gov- ern the network in Figure 5.5. The possible values for the parameters are bounded within

• 5. Preliminary Studies

52

physiological range, with association constants ki ∈ [10−6, 10−2] nM−1.s−1, dissociation con- stants k−i ∈ [10−6, 10−1] s−1 and catalytic constants kd ∈ [10−6, 1000] s−1. For Michaelis Menton parameters, KM ∈ (0, 109] nM and Vmax ∈ (0, 105] nM.s−1. Production and degra- dation rate constants are kept between 0 and 1. The estimated parameters and the total concentration of the various reactants are shown in Table 5.3 and Table 5.4. 5.2.6 Experimental Data For our pathway, we estimated the parameters based on three sets of experimental data, two of them being steady state values and the third, time-series data. Experiment Set 1 - Knock-down Experiments The first set of data is obtained from knockdown experiments. LNCaP cells are being transfected with small interfering RNA (siRNA) to knockdown the genes via RNA inter-

• ference. In this set of experiments, cells are treated with 20 nM, 50 nM and 100 nM of

siRNA for the proteins Akt and PDK1 to reduce their expression levels. Higher levels of Akt and PDK1 siRNA will lead to lower expression levels of Akt and PDK1 respectively (for experiments where only one level of siRNA is being used, such as Figure 5.7C, the presence

• f that treatment is indicated by a ‘+’ while ‘-’ means no treatment). The siRNAs are

transfected using the calcium phosphate method and the cells are then incubated for either 48 or 72 hours. The expression levels of selected proteins are then assayed and visualized via Western blotting. (A), (B) and (C) of Figure 5.7 show the expression levels of Akt and its activation levels under various knock-down conditions with darker bands indicating higher concentration levels. The three experiments are carried out independently. The effects of both treatments on the expression levels of Akt is already known, and such effects are seen from the experiments. Needless to say, Akt siRNA is suppose to reduce the total amount of Akt concentration, as can be seen in the figure. 20 nM of Akt siRNA is enough to reduce the total concentration of Akt to an almost negligible level. PDK1 is an activator of Akt, phosphorylating it at its Thr308 residue. Reducing PDK1 expression levels via siRNA will down-regulate Akt activation levels, as can be seen in Figure 5.7. The next set of data shows the expression levels of the various proteins involved in the ERK pathway. From Figure 5.8, it is observed that variations of PDK1 levels has not much effects on the activation of Raf. However, it has been reported in [ZM99] that Akt can regulate the ERK pathway by phosphorylating Raf at Ser259. Although this is not readily observed in (A) of Figure 5.8, its downstream effects on MEK and ERK can be seen, with increasing

• 5. Preliminary Studies

53

No Rate Equation Parameter Ref 1 V1[PI3K]/(K1 + [PI3K]) V1 = 130, K1 = 57800 Estimate 2 k2[PI3K] k2 = 1.5 × 10−4 Estimate 3 V3[PI3K∗]/K3 + [PI3K∗]) V3 = 7, K3 = 8000 Estimate 4 k4[PI3K∗][PIP2]/K4 + [PIP2]) k4 = 2.52, K4 = 85200 Estimate 5 k5[PTEN][PIP3]/K5 + [PIP3]) k5 = 15, K5 = 306 Estimate 6 k6[PIP3][AKTcyto] − k−6[AKT] k6 = 1.25 × 10−5, k−6 = 4.22 × 10−3 Estimate 7 k7[PDK1][AKT]/(K7 + [AKT]) k7 = 20, K7 = 80000 [BCC+00] 8 k8[PP2A][AKT∗]/(K8 + [AKT∗]) k8 = 35, K8 = 34000 Estimate 9 k9[PDK2][AKT∗]/(K9 + [AKT∗]) k9 = 20, K9 = 80000 [BCC+00] 10 k10[PP2A][AKT∗∗]/(K10 + [AKT∗∗]) k10 = 30, K10 = 64000 Estimate 11 k11[PP2A][AKT∗∗]/(K11 + [AKT∗∗]) k11 = 8.5, K11 = 45100 Estimate 12 V12[Ras]/(K12 + [Ras]) V12 = 0.825, K12 = 74500 Estimate 13 V13[Ras∗]/(K13 + [Ras∗]) V13 = 0.198, K13 = 5.52 Estimate 14 k14[Pak1∗][Raf]/(K14 + [Raf]) k14 = 0.18, K14 = 185 Estimate 15 k15[Ras∗][Raf]/(K15 + [Raf]) k15 = 0.47, K15 = 150 Estimate 16 k16[AKT∗∗][Raf∗]/(K16 + [Raf∗]) k16 = 3.3, K16 = 77.5 Estimate 17 V17[Raf∗]/(K17 + [Raf∗]) V17 = 66, K17 = 16.7 [SEJGM02] 18 k18[PDK1][MEK]/(K18 + [MEK]) k18 = 5.53 × 10−2, K18 = 14700 Estimate 19 k19[Raf∗][MEK]/(K19 + [MEK]) k19 = 3.5, K19 = 317 [SEJGM02] 20 k20[Raf∗][MEK∗]/(K20 + [MEK∗]) k20 = 2.9, K20 = 263 [SEJGM02] 21 k21[PDK1][MEK∗]/(K21 + [MEK∗]) k21 = 0.055, K21 = 14470 Estimate 22 k22[PP2A][MEK∗]/(K22 + [MEK∗]) k22 = 0.058, K22 = 2232 [SEJGM02] 23 k23[PP2A][MEK∗∗]/(K23 + [MEK∗∗]) k23 = 0.058, K23 = 60 [SEJGM02] 24 k24[MEK∗∗][ERK]/(K24 + [ERK]) k24 = 16, K24 = 1.46 × 105 [SEJGM02] 25 k25[MEK∗∗][ERK∗]/(K25 + [ERK∗]) k25 = 5.7, K25 = 5.21 × 104 [SEJGM02] 26 k26[MKP3][ERK∗]/(K26 + [ERK∗]) k26 = 0.3, K26 = 160 [SEJGM02] 27 k27[MKP3][ERK∗∗]/(K27 + [ERK∗∗]) k27 = 0.27, K27 = 60 [SEJGM02] 28 k28[ERK∗∗][p90RSK]/(K28 + [p90RSK]) k28 = 2.5 × 10−5, K28 = 97.6 Estimate 29 k29[ROS][p90RSK]/(K29 + [p90RSK]) k29 = 1.6 × 10−5, K29 = 81.4 Estimate 30 V30[p90RSK∗]/(K30 + [p90RSK∗]) V30 = 0.058, K30 = 5.8 Estimate 31 k31[p90RSK∗][Bad]/(K31 + [Bad]) k31 = 0.002, K31 = 346 Estimate 32 k32[Pak1∗][Bad]/(K32 + [Bad]) k32 = 0.46, K32 = 710 Estimate 33 k33[Pak1∗][Bad]/(K33 + [Bad]) k33 = 0.125, K33 = 1310 Estimate 34 k34[AKT∗∗][Bad]/(K34 + [Bad]) k34 = 3.46, K34 = 307 Estimate 35 V35[Bads112]/(K35 + [Bads112]) V35 = 5.8, K35 = 3450 Estimate 36 V36[Bads136]/(K36 + [Bads136]) V36 = 180, K36 = 797 Estimate 37 k37[Bad][Bcl − 2] − k−37[Bcl − 2.Bad] k37 = 4.78 × 10−2, k−37 = 1.35 × 10−2 Estimate 38 k38[Bax][Bcl − 2] − k−38[Bcl − 2.Bax] k38 = 2 × 10−3, k−38 = 0.02 [HCC+05] 39 k39[PI3K∗][Bax]/(K39 + [Bax]) k39 = 640, K39 = 38900 Estimate 40 V40[Baxcyto]/(K40 + [Baxcyto]) V40 = 79800, K40 = 29200 Estimate 41 k41[NOX5] k41 = 1.26 × 10−3 Estimate 42 k42[ROS] k42 = 1.7 × 10−3 Estimate 43 k43[ROS][PI3K∗][Pak1]/(K43 + [Pak1]) k43 = 0.33, K43 = 22600 Estimate 44 V44[Pak1∗]/(K44 + [Pak1∗]) V44 = 270, K44 = 895 Estimate

Table 5.3: Rate reactions and the parameters being used. The first order and second order rate constants are given in s−1 and nM −1.s−1 respectively while Michaelis Constants V and K are given in nM.s−1 and nM respectively.

levels of MEK and ERK activation when the expression level of Akt decreases (B to D). On the other hand, decreasing the levels of PDK1 does have an effect of attenuating MEK

• 5. Preliminary Studies

54

Reactant Concentration (nM) Ref PI3K 30 [CAZ04] PTEN 0.1 Estimate PIP2 7000 [CAZ04] AKT 200 [BDJF+05] PDK1 1000 Estimate PDK2 1000 Estimate PP2A 150 [KFK+03] NOX5 2000 Estimate PAK1 500 [LNLM04] ROS 2000 Estimate Ras 18900 [SEJGM02] Raf 66.4 [SEJGM02] MEK 36500 [SEJGM02] ERK 34900 [SEJGM02] MKP3 44000 Estimate P90RSK 5 [BJ00] Bad 100 Estimate Bax 100 [JXD+98] Bcl2 100 Estimate

Table 5.4: Initial concentration of the cellular components Figure 5.7: Akt and PDK1 knockdown on Akt expression and activation levels

and ERK activity. P90RSK plays an essential role in cell growth by activating several transcription factors, including the Na+/H+ exchangers. This protein is also being implicated by the activation

• 5. Preliminary Studies

55 Figure 5.8: Effects of the knockdown experiments on the various parts of the ERK pathway. (A) shows Raf expression and activation levels upon knocking down Akt and PDK1 while (B) shows that of MEK. (C) and (D) are two independent experiments for the activation levels of ERK.

• f the MAPK pathway, and hence is one of the proteins that are being assayed in the
• experiments. One of the downstream targets of P90RSK is Bad. P90RSK blocks Bad-

mediated cell death by phosphorylating Bad at Ser112 [TRDC99]. At the same time, Bad is also one of the downstream targets of Akt, being phosphorylated at Ser136 in a PI3K- dependent manner [TRDC99, ZM99].

Figure 5.9: Akt and PDK1 knockdown effects on P90RSK (A) and Bad (B, C, D) activation

Experiment Set 2 - PP2A Knockdown and Post Transfection Treatment The second set of experiments involves not only RNA interference (for this set, PP2A siRNAs are being used instead), but also the treatment of cells with reagents LY294002 and DPI. The cells are exposed to the various conditions for 1 hour before being assayed for protein expression.

• 5. Preliminary Studies

56

PP2A is a family of phosphatase that reverses the activation of the various proteins by removing the phosphate groups. Some of these proteins include Akt and MEK. It is also reported that PP2A deactivates ERK. However, this regulation is dominated by the highly selective activity of MAP Kinase Phosphatase 3 (MKP3) [Key00]. Nevertheless, PP2A is still crucial in the regulation of the ERK pathway via MEK dephosphorylation.

Figure 5.10: Effects of the various treatment on the activity of Akt

The LNCaP cell line lacks the active lipid phosphatase PTEN, a negative regulator of PIP3 and hence will display high levels of Akt activity [NLM+01]. With that, we can see that even in the absence of serum, there will be some activation of the Akt pathway, with its level probably being kept moderate by PP2A. The importance of PI3K in the Akt pathway can be seen in the treatment of the cells by LY294002, a selective PI3K inhibitor. In Figure 5.10, adding of LY294002 appears to abrogate Akt phosphorylation, hence inhibiting this pathway. DPI on the other hand, does not affect Akt activity much, showing that ROS does not have any profound impact on the activity of the Akt pathway. The role of PP2A in deactivating the Akt pathway is shown in Figure 5.10. For almost every experiment, its corresponding counterpart with PP2A siRNA appears to have a higher Akt phosphorylation, which would show that PP2A indeed plays a role in negating the effects of Akt activation. Next are the effects of PP2A knockdown and the various treatments on the components

• f the ERK pathway (Figure 5.11). The effects of serum on MEK activation is very obvious,

with its activity increasing as much as 3-fold in the presence of serum. The knocking down

• f PP2A also has the effect of increasing MEK phosphorylation. The addition of DPI will

affect the activation for both MEK and ERK, inhibiting their activity. This re-affirms the role of ROS as second messengers in signaling pathways, although we only show them in the context of the ERK pathway.

• 5. Preliminary Studies

57 Figure 5.11: Effects of the various treatments on MEK and ERK activity. Greyed out regions means that the protein types were not assayed for that particular experiment

Experiment Set 3 - Time Series Experiment The third and final set of experiments is aimed primarily at finding out how protein activ- ity changes with time under various conditions (Different combinations of serum, LY294002 and DPI). For the cells, they are initially kept in a serum starved state. They are then exposed to the different treatments and the readings of the various protein activation levels are taken at times 0 min, 10 min, 20 min, 30 min, 45 min, 60 min and 120 min after the start of the treatment. The proteins of interest in this set of experiments are MEK, ERK, P90RSK and Bad phosphorylation at Ser112. The relative activation levels of the proteins are shown in Figure 5.13 as triangles and squares. The relative endpoints are already ex- pected from the model, but having time profiles will enable us to have a tighter fit of the parameters. Using all these three sets of experiments, we then estimate the parameters of the model such that it can reproduce most, if not all, of the observations made. 5.2.7 Simulation The model is constructed using Cell Illustrator on a Pentium IV personal computer. To estimate the parameters, we broke the pathway down into modules separately and assigned them ranks according to the scheme presented in the previous sections. We then run the estimation algorithm, written in C++, on a PC cluster system. The algorithm was run a few times, and the most representative results (i.e. whose simulation profile matches the experimental ones the best) are used.

• 5. Preliminary Studies

58

Effects of PDK1 siRNA on MEK and ERK Activation The observation that prompted this investigation is the treatment of LNCaP cells with PDK1 siRNA. Usually not thought to interfere directly with the ERK pathway, decreased concentration of PDK1 shows a significant drop in both MEK and ERK activity, suggesting that PDK1 might have been affecting the ERK pathway via an MEK dependent manner. This effect has been captured in our model. Figure 5.12 shows the western blot of MEK and ERK activation with various levels of PDK1 siRNA treatment, as well as their corresponding simulation profiles. We simulate the effect of PDK1 siRNA treatment by reducing its total concentration from 1000 nM (Control) down to 0 nM (100 nM siRNA treatment). In the control experiment, the level of ERK phosphorylation is moderate but once siRNA treatment is being administered, it shows an acute inhibition, dropping to nearly 0 nM, as can be seen in our simulations.

Figure 5.12: Profiles of MEK and ERK activation levels compared to experimental data

Effects of LY294002 and DPI The experimental time series profiles are obtained mainly from cells that have been treated with LY294002 and DPI, both in the presence and absence of serum. Experiments show that both DPI and LY294002 affects the ERK pathway in similar ways, suggesting that other than the PDK1 and Akt interaction, there might still be some other forms of interaction. In our model, we realise this effect via an indirect action of activated PI3K. It has been reported in [CKM+00] that PI3K is one of the regulators of Pak1 activation, which in turns regulates Raf activation. The profiles of the simulation, compared to experimental data are shown in Figure 5.13. 5.2.8 Model Validation To validate the model and parameters, we simulate the system again, this time com- paring it with another separate set of experimental data, or the validation data set. This

• 5. Preliminary Studies

59 Figure 5.13: Profiles of the various activity levels (ratio of activated to total concentration for that protein type) under different treatment. Each graph shows the activity of the particular protein type as time progresses. Solid lines indicate the simulations where serum is present while dashed lines represent those where serum is absent. In addition, some of the profiles include data points from their corresponding experiments for ease of comparison. The triangles indicate experimental data points for the cells that are grown in serum while the squares indicate data points for those that are serum starved.

• 5. Preliminary Studies

60

experiment is a re-creation of the cell treatment with LY294002. However, due to experi- mental variations, the initial conditions of the various protein types are different and this is reflected in the simulation. This form of validation is a simplified form of model validation - the holdout technique [Koh95]. Although it is not the best validation method, we use it due to the little and sparse nature of data sets that are obtained from laboratory experiments. The profiles of the following activated proteins Raf, MEK, ERK, AKT and P90RSK are shown in Figure 5.14.

Figure 5.14: Profiles of the model simulation compared to the validation set

From the validation profiles, it is hard to conclude that the parameters are an exact fit. The AKT profile shows an almost exact fit while there are deviations in the activation levels

• f MEK and ERK. However, similarities in qualitative effects can still be seen, lending some

level of confidence to the structure of the model. However, more work still needs to be done to explain some of the dynamic profiles that can be seen from the experiments.

• 5. Preliminary Studies

61

5.3 Summary

We have looked at two different pathways that play important roles in cell differentiation as well as cancer and apoptosis - the Wnt pathways and the Akt pathway. The models are still being iteratively updated as more information and data becomes available to us. We have identified some of the important issues associated with the modeling of such pathways, namely data interpretation, model construction, parameter estimation and lastly, model validation.

• 6. FUTURE PLANS

To understand the issues associated with the modeling biological pathways, we have modeled two major pathways, the Wnt and the Akt pathways using the Hybrid Functional Petri Net methodology. We have looked into issues regarding data collection and interpre- tation, model building and more importantly, parameter estimation. Thereafter, we have formulated a roadmap (Figure 4.2) that will guide our research The first goal of our future work is on the development of proper validation techniques. We will mainly explore Historical Data Validation and Predictive Validation, formalizing the data representation and metrics in the process. This will allow us to quantify the correctness of the model and its parameters instead of using subjective experience such as visual comparison to determine whether or not a model reflects the biological phenomena. The next part will be the development of a hybrid modeling framework that will allow signaling pathways to interact with regulatory networks. This requires the incorporation of some graphical models such as the Hybrid Functional Petri Nets, exploiting not only their hybrid nature, but also the introduction of more modularity features so that the pathways can be modeled separately and pieced together to form a larger network. One additional issue which can be further explored is the epigenetic state of the gene regulatory networks. Most gene regulatory networks representation are static in the sense that the structure or the states of the genes do not change, and that if the available transcription factors are present, the gene will express its transcripts. However, this is just part of the picture as it still depends on the state of the gene, whether or not it is transcriptionally active before it allows transcription. Finally, we will look into ways the simulation results can be presented such that biologists can understand and reason about them, comparing the models with their hypotheses and experimental set-up. Upon completion of the above tasks, they will be then be used to iteratively build up the models of the Wnt and the Akt pathways. The models will then be used to predict cell activity so as to assess the feasibility as well as the accuracy of such computational models.

• 6. Future Plans

63

This will be done with the help of our collaborators from the Genome Institute of Singapore and the Department of Biochemistry, NUS.

BIBLIOGRAPHY

[ABI+01] Rajeev Alur, Calin Belta, Franjo Ivanˇ ci´ c, Vijay Kumar, Max Mintz, George J. Pappas, Harvey Rubin, and Jonathan Schug. Hybrid Modeling and Simulation

• f Biomolecular Networks. In Proc of the 4th International Workshop on Hybrid

Systems: Computation and Control, volume 2034 of Lecture Notes in Computer Science, pages 19–32. Springer-Verlag, 2001. [AGH+00] Rajeev Alur, Radu Grosu, Yerang Hur, Vijay Kumar, and Insup Lee. Modular Specification of Hybrid Systems in CHARON. In Proc of the 3rd Interna- tional Workshop on Hybrid Systems: Computation and Control, volume 1790

• f Lecture Notes in Computer Science, pages 6–19, 2000.

[AK03] Isil Aksan and M. Levent Kurnaz. A Computer-Based Model for the Regulation

• f Mitogen-Activated Protein Kinase (MAPK) Activation. J. Recept. Signal

Transduction, 23(2-3):197–209, 2003. [AL04] Liliana Attisano and Etienne Labb´

• e. TGFβ and Wnt pathway cross-talk. Can-

cer and Metastasis Reviews, 23:53–61, 2004. [ALR03] K Mark Ansel, Dong U Lee, and Anjana Rao. An epigenetic view of helper T cell differentiation. Nature Immunology, 4(7):616–623, 2003. [AMK99] Tatsuya Akutsu, Satoru Miyano, and Satoru Kuhara. Identification Of Genetic Networks From A Small Number Of Gene Expression Patterns Under The Boolean Network Model. Pac Symp of Biocomputing, 4:17–28, 1999. [AMK00] Tatsuya Akutsu, Satoru Miyano, and Satoru Kuhara. Inferring qualitative rela- tions in genetic networks and metabolic pathways. Bioinformatics, 16(8):727– 734, 2000. [AMP+03]

• M. Antoniotti, B. Mishra, C. Piazza, A. Policriti, and M. Simeoni. Modeling

Cellular Behavior with Hybrid Automata: Bisimulation and Collapsing. In Computational Methods in Systems Biology, volume 2602 of Lecture Notes in Computer Science. Springer-Verlag, 2003. [APUM02] Marco Antoniotti, Alberto Policriti, Nadia Ugel, and Bud Mishra. XS-Systems: eXtended S-Systems and Algebraic Differential Automata for Modeling Cellular

Bibliography 65

• Behavior. In Proceedings of the International Conference on High Performance

Computing, HiPC 2002, volume 2552 of Lecture Notes in Computer Science, pages 431–442. Springer-Verlag, Dec 2002. [APUM03] Marco Antoniotti, Alberto Policriti, Nadia Ugel, and Bud Mishra. Model Building and Model Checking for Biochemical Processes. Cell Biochemistry and Biophysics, 38:271–286, 2003. [BCC+00] Ricardo M. Biondi, Peter C.F. Cheung, Antonio Casamayor, Maria Deak, Richard A. Currie, and Dario R. Alessi. Identification of a pocket in the PDK1 kinase domain that interacts with PIF and the C-terminal residues of PKA. The EMBO Journal, 19(5):979–988, 2000. [BCK+03] Sukhdev S. Brar, Zachary Corbin, Thomas P. Kennedy, Richelle Hemendinger, Lisa Thornton, Bettina Bommarius, Rebecca S. Arnold, A. Richard Whorton, Anne B. Sturrock, Thomas P. Huecksteadt, Mark T. Quinn, Kevin Krenetsky, Kristia G. Ardie, J. David Lambeth, and John R. Hoidal. NOX5 NAD(P)H

• xidase regulates growth and apoptosis in DU 145 prostate cancer cells. Am J

Physiol Cell Physiol, 285(2):C353–C369, 2003. [BD02] Hamid Bolouri and Eric H. Davidson. Modeling DNA Sequence-Based cis- Regulatory Gene Networks. Developmental Biology, 246:2–13, 2002. [BDJF+05] Stanley F. Barnett, Deborah Defeo-Jones, Sheng Fu, Paula J. Hancock, Kath- leen M. Haskell, Raymond E. Jones, Jason A. Kahana, Astrid M. Kral, Karen Leander, Ling L. Lee, John Malinowski, Elizabeth M. McAvoy, Debbie D Na- has, Ronald G. Robinson, and Hans E. Huber. Identification and Character- ization of Pleckstrin Homology Domain Dependent and Isozyme Specific Akt

• Inhibitors. Biochem J, 385:399–408, 2005.

[BF00] Frances A. Brightman and David A. Fell. Differential feedback regulation of the MAPK cascade underlies the quantitative differences in EGF and NGF signaling in PC12 cells. FEBS Letters, 482(3):169–174, 2000. [BH02] James C. Bezdek and Richard J. Hathaway. Some Notes on Alternating Opti-

• mization. In Lecture Notes in Computer Science, volume 2275, pages 289–300.

Springer-Verlag, 2002. [BJ00] Ramesh R. Bhatt and James E. Ferrell Jr. Cloning and Characterization of Xenopus Rsk2, the Predominant p90 Rsk Isozyme in Oocytes and Eggs. J Biol Chem, 275(42):32983–32990, 2000.

Bibliography 66

[BRdJ+05] Gr´ egory Batt, Delphine Ropers, Hidde de Jong, Johannes Geiselmann, Michel Page, and Dominique Schneider. Qualitative Analysis and Verification of Hy- brid Models of Genetic Regulatory Networks: Nutritional Stress Response in Escherichia coli. In Proc of the 8th International Workshop on Hybrid Sys- tems: Computation and Control, volume 3414 of Lecture Notes in Computer Science, pages 134–150. Springer-Verlag, 2005. [BS02] Hans-Georg Beyer and Hans-Paul Schewefel. Evolutionary strategies - A com- prehensive introduction. In Natural Computing, pages 3–52. Springer-Verlag, 2002. [Car04] Luca Cardelli. Brane Calculi - Interactions of Biological Membranes. In Com- putational Methods in Systems Biology, volume 3082 of Lecture Notes in Com- puter Science, pages 257–280, 2004. [CAZ04] David Chodniewicz, Abdullatif M. Alteraifi, and Doncho V. Zhelev. Ex- perimental Evidence for the Limiting Role of Enzymatic Reactions in Chemoattractant-induced Pseudopod Extension in Human Neutrophils. J Biol Chem, 279(23):24460–24466, 2004. [CFT01] Daniel L. Cook, Joel F. Farley, and Stephen J. Tapscott. A basis for a visual language for describing, archiving and analyzing functional models of complex biological systems. Genome Biology, 2(4):research0012.1–0012.10, 2001. [CKM+00]

• A. Chaudhary, W.G. King, M.D. Mattaliano, J.A. Frost, B. Diaz, D.K. Morri-

son, M.H. Cobb, M.S. Marshall, and J.S. Brugge. Phosphatidylinositol 3-kinase regulates Raf1 through Pak phosphorylation of serine 338. Current Biology, 10(9):551–554, 2000. [CSK+03] Kwang-Hyun Cho, Sung-Young Shin, Hyun Woo Kim, Olaf Wolkenhauer, Brian McFerran, and Walter Koch. Mathematical Modeling of the Influence

• f RKIP on the ERK Signaling Pathway. In Proceedings of the 1st Interna-

tional Workshop on Computational Methods in Systems Biology, volume 2602

• f Lecture Notes In Computer Science, pages 127–141. Springer-Verlag, 2003.

[CT91] Thomas M. Cover and Joy A. Thomas. Elements of Information Theory. John- Wiley, 1991. [dJ02] Hidde de Jong. Modeling and Simulation of Genetic Regulatory Systems: A Literature Review. Journal of Computational Biology, 9(1):67–103, 2002. [dJGH+04] Hidde de Jong, Jean-Luc Gouz´ e, C´ eline Hernandez, Michel Page, Tewfik Sari, and Johannes Geiselmann. Qualitative Simulation of Genetic Regulatory

Bibliography 67

Networks Using Piecewise-Linear Models. Bulletin of Mathematical Biology, 66(2):301–340, 2004. [DK03] Vincent Danos and Jean Krivine. Formal Molecular Biology done in CCS-R. In BIOCONCUR ’03, Electronic Notes in Theoretical Computer Science, 2003. [DNF+03] Atsushi Doi, Masao Nagasaki, Sachie Fujita, Hiroshi Matsuno, and Saturo

• Miyano. Genomic object net: Ii. modelling biopathways by hybrid functional

petri net with extension. Applied Bioinformatics, 2(3):185–188, 2003. [DRO+02] Eric H. Davidson, Jonathan P. Rast, Paola Oliveri, Andrew Ransick, Cristina Calestani, Chiou-Hwa Yuh, Takuya Minokawa, Gabriele Amore, Veronica Hin- man, C´ e Arenas-Mena, Ochan Otim, C. Titus Brown, Carolina B. Livi, Pei Yun Lee, Roger Revilla, Maria J. Schilstra, Peter J.C. Clarke, Alistair G. Rust, Zhengjun Pan, Maria I. Arnone, Lee Rowen, R. Andrew Cameron, David R. McClay, Leroy Hood, and Hamid Bolouri. A Provisional Regulatory Gene Network for Specification of Endomesoderm in the Sea Urchin Embryo. Devel-

• pmental Biology, 246:162–190, 2002.

[DWVL98] Ami A. Deora, Terrance Win, Bart Vanhaesebroeck, and Harry M. Lander. A Redox-triggered Ras-Effector Interaction - Recruitment of Phosphatidylinositol 3-Kinase to Ras by Redox Stress. J Biol Chem, 273(45):29923–29928, 1998. [ESM04] E.J.Crampin, S. Schnell, and P.E. McSharry. Mathematical and computational techniques to deduce complex biochemical reaction mechanisms. Progress in Biophysics and Molecular Biology, 86:77–112, 2004. [FBV00] Martin Fussenegger, James E. Bailey, and Jeffrey Varner. A mathematical model of caspase function in apoptosis. Nature Biotechnology, 18:768–774, July 2000. [FDG+03] Daniel B. Forger, Dennis A. Dean II, Katherine Gurdziel, Jean-Christophe Leloup, Choogon Lee, Charlotte Von Gall, Jean-Pierre Etchegaray, Richard E. Kronauer, Albert Goldbeter, Charles S. Peskin, Megan E. Jewett, and David R.

• Weaver. Development and Validation of Computational Models for Mammalian

Circadian Oscillators. OMICS A Journal of Integrative Biology, 7(4):387–400, 2003. [FH03] Andrew Finney and Michael Hucka. Systems Biology Markup Lan- guage (SBML) Level 2: Structures and Facilities for Model Definitions. http://sbml.org/index.psp, June 2003.

Bibliography 68

[FHC05] Christine Farrar, Carolyn R. Houser, and Steven Clarke. Activation of the PI3K/Akt signal transduction pathway and increased levels of insulin receptor in protein repair deficient mice. Aging Cell, 4:1–12, Feb 2005. [FHL+04] Anthony Finkelstein, James Hetherington, Linzhong Li, Ofer Margoninski, Pe- ter Saffrey, Rob Seymour, and Anne Warner. Computational challenges of systems biology. IEEE Computer, 37(5):26–33, 2004. [FLNP00] Nir Friedman, Michal Linial, Iftach Nachman, and Dana Pe’er. Using Bayesian networks to analyze expression data. In 4th Annual International Conference

• n Computational Molecular Biology, 2000.

[Fri04] Nir Friedman. Inferring Cellular Networks Using Probabilistic Graphical Mod-

• els. Science, 303:799–805, Feb 2004.

[Gen00] Gene Ontology Consortium. Gene Ontology: tool for the unification of biology. Nature Genetics, 25:25–29, May 2000. [GF05] Timothy S. Gardner and Jeremiah J. Faith. Reverse-engineering transcription control networks. Physics of Life Reviews, 2:65–88, 2005. [Gil76] Daniel T. Gillespie. A general method for numerically simulating the stochas- tic time evolution of coupled chemical reactions. Journal of Computational Physics, 22(4):403–434, 1976. [Gil92] Daniel T. Gillespie. A rigorous derivation of the chemical master equation. Physica A, 188:404–425, 1992. [GKV01] Hartmann Genrich, Robert K¨ uffner, and Klaus Voss. Executable Petri Net models for the analysis of metabolic pathways. International Journal on Soft- ware Tools for Technology Transfer, 3(4):394–404, 2001. [GP98] Peter J.E. Goss and Jean Peccoud. Quantitative modeling of stochastic systems in molecular biology using stochastic Petri nets. Proc Natl Acad Sci USA, 95:6750–6755, June 1998. [GT04a]

• R. Ghosh and C. Tomlin. Symbolic reachable set computation of piecewise

affine hybrid automata and its application to biological modelling: Delta-Notch protein signalling. Syst Biol, 1(1):170–183, 2004. [GT04b] Ronojoy Ghosh and Claire Tomlin. An Algorithm for Reachability Computa- tions on Hybrid Automata Models of Protein Signaling Networks. In Proceed- ings of the 43rd IEEE Conference on Decision and Control, 2004.

Bibliography 69

[HCC+05] Fei Hua, Melanie G. Cornejo, Michael H. Cardone, Cynthia L. Stokes, and Dou- glas A. Lauffenburger. Effects of bcl-2 levels on fas signaling-induced caspase-3 activation: Molecular genetic tests of computational model predictions. Jour- nal of Immunology, 175:985–995, 2005. [He03] Xi He. A Wnt-Wnt Situation. Developmental Cell, 4:791–797, 2003. [Hen96] Thomas A. Henzinger. The Theory of Hybrid Automata. In Symp on Logic in Computer Science, pages 278–292. IEEE Computer Society Press, 1996. [HFB+04]

• M. Hucka, A. Finney, B.J. Bornstein, S.M. Keating, B.E. Shapiro, J. Matthews,

B.L. Kovitz, M.J. Schilstra amd A. Funahashi, J.C. Doyle, and H. Kitano. Evolving a lingua franca and associated software infrastructure for compu- tational systems biology: the Systems Biology Markup Language (SBML)

• project. Systems Biology, 1(1):41–53, June 2004.

[HFS+01] Michael Hucka, Andrew Finney, Herbert Sauro, Hamid Bolouri, John Doyle, and Hiroaki Kitano. The ERATO Systems Biology Workbench: An Integrated Environment for Multiscale and Multitheoretic Simulations in Systems Biology. Foundations in Systems Biology, pages 125–139, 2001. [HFS+02]

• M. Hucka, A. Finney, H.M. Sauro, H. Bolouri, J. Doyle, and H. Kitano. The

ERATO Systems Biology Workbench: enabling interaction and exchange be- tween software tools for computational biology. In Proc Pacific Symp Biocom- put, pages 450–461, 2002. [HJ96] Chi Ying F. Huang and James E Ferrell Jr. Ultrasensitivity in the mitogen- activated protein kinase cascade. Proc Natl Acad Sci, 93(19):10078–10083, Sept 1996. [HKN+03] Mariko Hatakeyama, Shuhei Kimura, Takashi Naka, Takuji Kawasaki, Noriko Yumoto, Mio Ichikawa, Jae-Hoon Kim, Kazuki Saito, Mihoro Saeki, Mikako Shirouzu, Shigeyuki Yokoyama, and Akihiko Konagaya. A computational model on the modulation of mitogen-activated protein kinase (MAPK) and Akt pathways in heregulin-induced ErbB signaling. Biochem J, 373(2):451– 463, July 2003. [HKW04] Monika Heiner, Ina Koch, and J¨ urgen Will. Model validation of biological pathways using petri nets - demonstrated for apoptosis. BioSystems, 75:15–28, 2004. [HLC+80] J.S. Horoszewicz, S.S. Leong, T.M. Chu, Z.L. Wajsman, M. Friedman, L. Pap- sidero, U. Kim, L.S. Chai, S. Kakati, S.K. Arya, and A.A. Sandberg. The

Bibliography 70

LNCaP cell line: a new model for studies on human prostatic carcinoma. Prog Clin Biol Res, 37:115–132, 1980. [HMM+02] Time Hesterberg, David S. Moore, Shaun Monaghan, Ashley Clipson, and Rachel Epstein. Bootstrap Methods and Permutation Tests, chapter 14, pages 14.1–14.70. W.H. Freeman and Company, 5 edition, 2002. [HQA+00] Priti Hegde, Rong Qi, Kristie Abernathy, Cheryl Gay, Sonia Dharap, Renee Gaspard, Julie Earle-Hughes, Erik Snesrud, Norman Lee, and John Quack- enbush. A Concise Guide to cDNA Microarray Analysis. BioTechniques, 29(3):548–562, 2000. [HTF01] Trevor Hastie, Robert Tibshirani, and Jerome Friedman. The Elements of Statistical Learning: Data Mining, Inference and Prediction. Springer-Verlag, 2001. [Hus03]

• D. Husmeier. Reverse Engineering of genetic networks with Bayesian networks.

Biochemical Society Transactions, 31(6):1516–1518, 2003. [IGH01] Trey Ideker, Timothy Galitski, and Leroy Hood. A new approach to decoding life: systems biology. Annu Rev Genomics Hum Genet, 2:343–372, 2001. [IWR+05] Takao Inoue, Minqin Wang, Ted O. Ririe, Jolene S. Fernandes, and Paul W.

• Sternberg. Transcriptional network underlying caenorhabditis elegans vulval
• development. PNAS, 102(14):4972–4977, Apr 2005.

[JXD+98] Juliane M J¨ urgensmeier, Zhihua Xie, Quinn Deveraux, Lisa Ellerby, Dale Bre- desen, and John C. Reed. Bax directly induces release of cytochrome c from isolated mitochondria. Cell Biology, 95(9):4997–5002, 1998. [JZD+02] E.H. Jho, T. Zhang, C. Domon, C.K. Joo, J.N. Freund, and F. Costantini. Wnt/beta-catenin/Tcf signaling induces the transcription of Axin2, a negative regulator of the signaling pathway. Mol Cell Biol, 22(4):1172–1183, 2002. [Kas05] Samuel Kaski. From Learning Metrics Towards Dependency Exploration. In Proc WSOM 05, 2005. To appear. [KEG] KEGG: Kyoto Encyclopedia

• f

Genes and Genomes. http://www.genome.ad.jp/kegg/. [Key00] Stephen M. Keyse. Protein phosphatases and the regulation of mitogen- activated protein kinase signalling. Current Opinion in Cell Biology, 12(2):186– 192, 2000.

Bibliography 71

[KFK+03] Shinichi Kikuchi, Kenji Fujimoto, Noriyuki Kitagawa, Taro Fuchikawa, Michiko Abe, Kotaro Oka, Kohtaro Takei, and Masaru Tomita. Kinetic simulation of signal transduction system in hippocampal longterm potentiation with dynamic modeling of protein phosphatase 2A. Neural Networks, 16(9):1389–1398, 2003. [Kho00] Boris N. Kholodenko. Negative feedback and ultrasensitivity can bring about

• scillations in the mitogen-activated protein kinase cascades. Eur J Biochem,

267:1583–1588, 2000. [KIM03] Sun Yong Kim, Seiya Imoto, and Satoru Miyano. Inferring gene networks from time series microarray data using dynamic Bayesian networks. Briefings in Bioinformatics, 4(3):228–235, Sept 2003. [Kit02] Hiroaki Kitano. Systems Biology: A Brief Overview. Science, 295:1662–1664, Mar 2002. [Kit03] Hiroaki Kitano. A graphical notation for biochemical networks. Biosilico, 1(5):169–176, Nov 2003. [KNK+04] Yuko Kawakami, Hajime Nishimoto, Jiro Kitaura, Mari Maede-Yamamoto, Roberta M. Kato, Dan R. Littman, David J. Rawlings, and Toshiaki Kawakami. Protein Kinase Cβ II Regulates Akt Phosphorylation on Ser-473 in an Cell Type and Stimulus specific Fashion. J Biol Chem, 279(46):47720–47725, 2004. [Koh95] Ron Kohavi. A Study of Cross-Validation and Bootstrap for Accuracy Esti- mation and Model Selection. In Proc 14th Int Joint Conference on Artificial Intelligence, pages 1137–1143, 1995. [Koh99] Kurt W. Kohn. Molecular Interaction Map of the Mammalian Cell Cycle Control and DNA Repair Systems. Molecular Biology of the Cell, 10:2703– 2734, Aug 1999. [Kol00] Walter Kolch. Meaningful relationships: the regulation

• f

the Ras/Raf/MEK/ERK pathway by protein interactions. Biochem J, 351:289– 305, 2000. [KS04] Samuel Kaski and Janne Sinkkonen. Principle of Learning Metrics for Exploratory Data Analysis. Journal of VLSI Signal Processing Systems, 37(2):177–188, 2004. [LDB05] William J.R. Longabaugh, Eric H. Davidson, and Hamid Bolouri. Computa- tional representation of developmental genetic regulatory networks. Develop- mental Biology, 283:1–16, 2005.

Bibliography 72

[LFS98] Shoudan Liang, Stefanie Fuhrman, and Roland Somogyi. REVEAL, A General Reverse Engineering Algorithm for Inference of Genetic Network Architectures. Pac Symp on Biocomputing, 3:18–29, 1998. [LK05] Olle Lindvall and Zaal Kokaia. Stem cell therapy for human brain disorder. Kidney International, 68:1937–1939, 2005. [LNLM04] Tsui-Han Loo, Yuen-Wai Ng, Louis Lim, and Ed Manser. GIT1 Activates p21- Activated Kinase through a Mechanism Independent of p21 Binding. Molecular and Cell Biology, 24(9):3849–3859, May 2004. [LSK+03] Ethan Lee, Adrian Salic, Roland Kr¨ uger, Reinhart Heinrich, and Marc W.

• Kirschner. The Roles of APC and Axin Derived from Experimental and The-
• retical Analysis of the Wnt Pathway. PLoS Biology, 1(1):116–132, Oct 2003.

[LT04] Patrick Lincoln and Ashish Tiwari. Symbolic Systems Biology: Hybrid Mod- eling and Analysis of Biological Networks. In R. Alur and G. Pappas, editors, Hybrid Systems: Computation and Control HSCC, volume 2993 of LNCS, pages 660–672. Springer, Mar 2004. [MDMK00] G. Marnellos, G.A. Deblandre, E. Mjolsness, and C. Kintner. Delta-Notch lateral inhibitory patterning in the emergence of ciliated cells in Xenopus: ex- perimental observations and a gene network model. In Proc Pacific Symp Biocomput, pages 329–340, 2000. [MDNM00] Hiroshi Matsuno, Atsushi Doi, Masao Nagasaki, and Satoru Miyano. Hybrid Petri Net Representation of Gene Regulatory Network. Pac Symp Biocomput, pages 341–352, 2000. [Mic] G. Michal. Roche Applied Science: Biochemical pathways chart. http://www.expasy.ch/cgi-bin/show thumbnails.pl. [MM99] Kevin Murphy and Saira Mian. Modelling Gene Expression Data using Dy- namic Bayesian Networks. Technical report, Computer Science Division, Uni- versity of California, Berkeley, CA, 1999. [MMB03] Carmen G. Moles, Pedro Mendes, and Julio R. Banga. Parameter Estimation in Biochemical Pathways: A Comparison of Global Optimization Methods. Genome Research, 13(11):2467–2474, Nov 2003. [MNB+04] Adam A. Margolin, Ilya Nemenman, Katia Basso, Chris Wiggins, Gustavo Stolovitzky, Riccardo Dalla Favera, and Andrea Califano. ARACNE: An Al- gorithm for the Reconstruction of Gene Regulatory Networks in a Mammalian Cellular Context. eprint arXiv:q-bio/0410037, Oct 2004.

Bibliography 73

[MP03]

• B. Mishra and A. Policriti. Systems biology and automata. In Proceedings
• f the 3rd Workshop on Computation of Biochemical Pathways and Genetic

Networks, Oct 2003. [MPO95] Thomas Mestl, Erik Plahte, and Stig W. Omholt. A mathematical framework for describing and analysing gene regulatory networks. Journal of Theoretical Biology, 176(2):291–300, Sept 1995. [MSS05] Wolfgang Marwan, Arumugam Sujatha, and Christine Starostzik. Recon- structing the regulatory network controlling commitment and sporulation in Physarum polycephalum based on hierarchical Petri Net modelling and simu-

• lation. J Theor Biol, 236:349–365, 2005.

[MTA+03] Hiroshi Matsuno, Yukiko Tanaka, Hitoshi Aoshima, Atsushi Doi, Mika Matsui, and Satoru Miyano. Biopathways representation and simulation on hybrid functional Petri net. In Silico Biology, 3(3):389–404, 2003. [Mur01] Kevin P. Murphy. An introduction to graphical models. Technical Report, Intel Research Technical Report, 2001. [NA00] Karleen M. Nicholson and Neil G. Anderson. The protein kinase B/Akt sig- naling pathway in human malignancy. Cellular Signaling, 14(5):381–395, May 2000. [NDMM03] Masao Nagasaki, Atsushi Doi, Hiroshi Matsuno, and Satoru Miyano. Genomic Object Net: I. A platform for modelling and simulating biopathways. Applied Bioinformatics, 2(3):181–184, 2003. [Nic] Donald Nicholson. IUBMB-Nicholson Metabolic Maps, Minimaps and Ani-

• maps. http://www.tcd.ie/Biochemistry/IUBMB-Nicholson/.

[NLM+01] Alexandre Nesterov, Xiaojun Lu, Gary J Miller, Yuri Ivashchenko, and An- drew S. Kraft. Elevated AKT Activity Protects the Prostate Cancer Cell Line LNCaP from TRAIL-induced Apoptosis. J Biol Chem, 276(14):10767–10774, 2001. [NNK+04]

• M. Nakano, R. Noda, H. Kitakaze, H. Matsuno, and H. Miyano. XML pathway

file conversion between Genomic Object Net and SBML. In Proc of The 3rd International Conference on Information, pages 585–588, 2004. [Nob02] Denis Noble. Modelling the heart: insights, failures and progress. BioEssays, 24(12):1155–1163, 2002. [NP00] Stephane Noselli and Norbert Perrimon. Signal Transduction: Are there Close Encounters Between Signaling Pathways? Science, 290(5489):68–69, Oct 2000.

Bibliography 74

[Nus05] Roel Nusse. The Wnt Homepage. http://www.stanford.edu/∼rnusse/wntwindow.html, 2005. [PBM+00] Kathleen I. Pinson, Jane Brennan, Susan Monkley, Brian J. Avery, and William C. Skarnes. An LDL-receptor related protein mediates Wnt signalling in mice. Nature, 407:535–538, Sep 2000. [PC05] Andrew Phillips and Luca Cardelli. A Graphical Representation for the Stochastic Pi-Calculus. In Proceedings of Concurrent Models in Molecular Bi-

• logy (Bioconcur’05), 2005.

[PP00] Mark Peifer and Paul Polakis. Wnt Signaling in Oncogenesis and Embyrogen- esis - a Look Outside the Nucleus. Science, 287:1606–1609, 2000. [PP05] Daniela Palacios and Pier Lorenzo Puri. The epigenetic network regulating muscle development and regeneration. Journal of Cellular Physiology, Epub ahead of print, Sept 2005. [PRSS01] Corrado Priami, Aviv Regev, Ehud Shapiro, and William Silverman. Applica- tion of a stochastic name-passing calculus to representation and simulation of molecular processes. Information Processing Letters, 80(1):25–31, 2001. [PSA+04] Yves Pommier, Olivier Sordet, Smitha Antony, Richard L. Hayward, and Kurt W. Kohn. Apoptosis defects and chemotherapy resistance: molecular interaction maps and networks. Oncogene, 23(16):2934–2949, 2004. [PSH03] Chang Shin Park, Ian C. Schneider, and Jason M. Haugh. Kinetic Analysis

• f Platelet-derived Growth Factor Receptor/Phosphoinositide 3-Kinase/Akt

Signaling in Fibroblasts. J Biol Chem, 278(39):37064–37073, 2003. [PSS+03] Tomi Pasanen, Janna Saarela, Ilana Saarikko, Teemu Toivanen, Martti Tolva- nen, Mauno Vihinen, and Garry Wong. DNA Microarray Data Analysis. CSC

• Scientific Computing Ltd, 1st edition, 2003.

[PWM03] J.W. Pinney, D.R. Westhead, and G.A. McConkey. Petri Net representations in systems biology. Biochemical Society Transactions, 31(6):1513–1515, 2003. [PZHK05] Louchka Popova-Zeugmann, Monika Heiner, and Ina Koch. Time Petri Nets for Modelling and Analysis of Biochemical Networks. Fundamenta Informaticae, 67:149–162, 2005. [RA03] Christopher V. Rao and Adam P. Arkin. Stochastic chemical kinetics and the quasi-steady state assumption: Application to the Gillespie algorithm. Journal

• f Chemical Physics, 118(11):4999–5010, 2003.

Bibliography 75

[Rae02]

• L. Raeymaekers. Dynamics of Boolean Networks Controlled by Biologically

Meaningful Functions. J Theor Biol, 218:331–341, 2002. [RN05] Tim C. Roloff and Ulrike A. Nuber. Chromatin, epigenetics and stem cells. European Journal of Cell Biology, 84:123–135, 2005. [RPS+04] Aviv Regev, Ekaterina M. Panina, William Silverman, Luca Cardelli, and Ehud

• Shapiro. BioAmbients: An abstraction for biological compartments. Theoret-

ical Computer Science, Special Issue on Computational Methods in Systems Biology, 325(1):141–167, 2004. [RS04a] Aviv Regev and Ehud Shapiro. The π-calculus as an abstraction for biomolec- ular systems. Modelling in Molecular Biology, pages 219–266, 2004. [RS04b] Johann Rohwer and Jacky Snoep. Concepts in Systems Biology: Kinetics, Control and Regulation of Cellular Systems, Apr 2004. [RSD03]

• W. Reik, F. Santos, and W. Dean. Mammalian epigenomics: reprogramming

the genome for development and therapy. Theriogenology, 59:21–32, 2003. [RSS01] Aviv Regev, William Silverman, and Ehud Shapiro. Representation and sim- ulation of biochemical processes using the π-calculus process algebra. In Proc Pacific Symp Biocomput, volume 6, pages 459–470, 2001. [Run04] Thomas Runge. Application of Coloured Petri Nets in Systems Biology. In Proceedings of the 5th Workshop and Tutorial on Practical Use of Coloured Petri Nets and the CPN Tools, pages 77–96, Oct 2004. [RY00] T.P. Runarsson and X. Yao. Stochastic ranking for constrained evolutionary

• ptimization. IEEE Trans Evol Comput, 4:284–294, 2000.

[RZS+01]

• H. Peter Reusch, Sven Zimmermann, Michael Schaefer, Martin Paul, and Karin
• Moelling. Regulation of Raf by Akt Controls Growth and Differentiation in

Vascular Smooth Muscle Cells. J Biol Chem, 276(36):33630–33637, 2001. [Sar03] Robert G. Sargent. Verification and Validation of Simulation Models. In Proc

• f the Winter Simulation Conference, 2003.

[Sav69a] M.A. Savageau. Biochemical Systems Analysis. I: Some mathematical prop- erties of the rate law for the component enzymatic reactions. J Theor Biol, 25:365–369, 1969. [Sav69b] M.A. Savageau. Biochemical Systems Analysis. II: The steady-state solutions for an n-pool system using a power-law approximation. J Theor Biol, 25:270– 279, 1969.

Bibliography 76

[Sav70] M.A. Savageau. Biochemical Systems Analysis 3: Dynamic Solutions using a power-law approximation. J Theor Biol, 26:215–226, 1970. [Sch79]

• S. Schlesinger. Terminology for Model Credibility. Simulation, 32(3):103–104,

1979. [SDZ02] Ilya Shmulevich, Edward D. Dougherty, and Wei Zhang. From Boolean to Probabilistic Boolean Networks as Models of Genetic Regulatory Networks. Proceedings of the IEEE, 90(11):1778–1792, Nov 2002. [SEJGM02] Birgit Schoeberl, Claudia Eichler-Jonsson, Ernst Dieter Gilles, and Gertraud M¨

• uller. Computational modeling of the dynamics of the MAP kinase cascade

activated by surface and internalized EGF receptors. Nature Biotechnology, 20(4):370–375, 2002. [SHNI02] S.V. Sabau, S. Hashimoto, Y. Nemoto, and S. Ihara. Cell Simulation for Circa- dian Rhythm Based on Michaelis-Menten Model. Journal of Biological Physics, 28:465–469, 2002. [SHW04] R.H. Smallwood, W.M.L. Holcombe, and D.C. Walker. Development and val- idation of computational models of cellular interaction. Journal of Molecular Histology, 35(7):659–665, 2004. [SI04] Naoya Sugimoto and Hitoshi Iba. Inference of Gene Regulatory Networks by Means of Dynamic Differential Bayesian Networks and Nonparametric Regres-

• sion. Genome Informatics, 15(2):121–130, 2004.

[Ski97] Steve V. Skiena. The Algorithm Design Manual. Springer, 1st edition, 1997. [SMS+00]

• A. Sch¨

urmann, A.F. Mooney, L.C. Sanders, M.A. Sells, H.G. Wang, J.C. Reed, and G.M. Bokoch. p21-Activated Kinase 1 Phosphorylates the Death Agonist Bad and Protects Cells from Apoptosis. Mol Cell Biol, 20(2):453–461, 2000. [SPB01]

• R. Srivastava, M.S. Peterson, and W.E. Bentley. Stochastic Kinetic Analysis
• f the Escherichia coli Stress Circuit Using σ32-Targeted Antisense. Biotechnol

Bioeng, 75(1):120–129, Oct 2001. [SS92]

• F. Shiraishi and M.A. Savageau. The tricarboxylic acid cycle in Dictyostelium

discoideum I. Formulation of alternative kinetic representations. J Biol Chem, 267(32):22912–22918, 1992. [SS96]

• R. Somogyi and C.A. Sniegoski. Modeling the complexity of genetic networks:

Understanding multigenic and pleitropic regulation. Complexity, 1(6):45–63, 1996.

Bibliography 77

[THT+99] Masaru Tomita, Kenta Hashimoto, Kouichi Takahashi, Thomas Simon Shimizu, Yuri Matsuzaki, Fumihiko Miyoshi, Kanako Saito, Sakura Tanida, Katsuyuki Yugi, J. Craig Venter, and Clyde A. Hutchinson III. E-CELL: soft- ware environment for whole-cell simulation. Bioinformatics, 15(1):72–84, 1999. [TK05] Hiroaki Takagi and Kunihiko Kaneko. Dynamical Systems basis of metamor- phosis: diversity and plasticity of cellular states in reaction diffusion network. J Theor Biol, 234:173–186, 2005. [TRDC99] Yi Tan, Hong Ruan, Matthew R. Demeter, and Michael J. Comb. p90RSK Blocks Bad-mediated Cell Death via a Protein Kinase C-dependent Pathway. J Bio Chem, 274(49):34859–34867, 1999. [TSD04] C.M. Tan, S. Somani, and P. Dhar. Modeling and simulation of biological systems with stochasticity. In Silico Biology, 4(3), 2004. [TYH+02] Kouichi Takahashi, Katsuyuki Yugi, Kenta Hashimoto, Yohei Yamada, Christopher J.F. Pickett, and Masaru Tomita. Computational Challenges in Cell Simulation: A Software Engineering Approach. IEEE Intelligent Systems, 17(5):64–71, 2002. [VAM03] Michael T. Veeman, Jeffrey D. Axelrod, and Randall T. Moon. A Second Canon: Functions and Mechanisms of β-Catenin- Independent Wnt Signaling. Developmental Cell, 5:367–377, Sep 2003. [Voi00] Eberhard O. Voit. Computational Analysis of Biochemical Systems : A Prac- tical Guide for Biochemists and Molecular Biologists. Cambridge University Press, 2000. [WBI99] Gezhi Weng, Upinder S. Bhalla, and Ravi Iyengar. Complexity in Biological Signaling Systems. Science, 284(5411):92–96, Apr 1999. [WM03] Hsien Yu Wang and Craig C. Malbon. Wnt Signaling, Ca2+, and Cyclic GMP: Visualizing Frizzled Functions. Science, 300:1529–1530, 2003. [WN98] Andreas Wodarz and Roel Nusse. Mechanisms of Wnt signaling in development. Annu Rev Cell Dev Biol, 14:59–88, 1998. [WSW+05] O. Wolkenhauer, S.N. Sreenath, P. Wellstead, M. Ullah, and K.H. Cho. A systems- and signal- oriented approach to intracellular dynamics. Biochemical Society Transactions, 33(3):507–515, 2005. [YMM+05] Chen-Hsiang Yeang, H. Craig Mak, Scott McCuine, Christopher Workman, Tommi Jaakkola, and Trey Ideker. Validation and refinement of gene-regulatory

Bibliography 78

pathways on a network of physical interactions. Genome Biology, 6:R62.1 – R62.10, 2005. [YTY+99] Terry P. Yamaguchi, Shinji Takada, Yoshiaki Yoshikawa, Nongying Wu, and Andrew P. McMahon. T(Brachyury) is a direct target of Wnt3a during paraxial mesoderm specification. Genes and Development, 13(24):3185–3190, 1999. [ZM99] Sven Zimmermann and Karin Moelling. Phosphorylation and Regulation of Raf by Akt (Protein Kinase B). Science, 286(5445):1741–1744, 1999. [ZOS03] Ionela Zevedei-Oancea and Stefan Schuster. Topological Analysis of Metabolic Networks Based on Petri Net Theory. In Silico Biology, 3:323–345, 2003.