WO2005111905A2 - Prediction of the dynamic behavior of a biological system - Google Patents

Abstract

A method for predicting the dynamic behavior of a biological system and/or the quantitative and qualitative prediction of the dependencies of the biological objects included in said biological system is provided. According to the method of the present invention, a computational model of the biological system is generated. Said model comprises one or more models, objects and associated model reactions, said model reactions being provided with kinetic parameters. The computational model may be provided as a hierarchically object oriented model. Further, a computer program for performing the method is provided.

Description

Prediction of the dynamic behavior of a biological system

The present invention relates to a method for the prediction of the dynamic behavior of a biological system and the quantitative and qualitative prediction of the dependencies of the biological objects included in said biological system.

Current functional genomics and health applications lack computational strategies that are capable of explaining and simulating biological processes that are the basis of most of disease and developmental processes. On the one hand, the field of kinetic modelling is well developed in the area of small metabolic and signalling models (for example kinetic models for glycolysis, Wnt pathway). These models rely on the exact identification of the kinetic parameters and the subsequent experimental measurement of these parameters. The exact calculation and modelling of metabolic processes via systems of differential and algebraic equations require the implementation of kinetic rate equations known in the literature for the relevant steps in the network and the literature values of all parameters entering these equations. This has been done for several models, for example glycolysis. Such exact calculations are not feasible for large models simply because of two main reasons: 1.) the number of model parameters is too high, 2.) the measurement of a realistic number of parameters that determine model behaviour is too expensive. On the other hand, the fields of genomics, proteomics and bioinformatics yield a large amount of primarily unstructured data and a set of computational and statistical tools to analyse these data that give information on many biological objects in parallel. In contrast to the kinetic modelling approach that delivers exact models for the quantitative prediction of small biological systems, the bioinformatics approach delivers information for the comparison of the qualitative prediction of systems of any size.

In other words, a challenging goal in biology is the integration of the mass of knowledge about genomic sequence data, functional genomics, proteomics and metabolics into large models of cellular systems. Figure 1 illustrates several fundamental cellular processes, which form the complex interaction network of a cellular system. A computational approach offers the possibility to use simulations for the prediction of the dynamical behavior of biological systems according to the defined models, and to test the validity of the underlying assumptions (Kitano 2002). Therefore the construction of computer executable models is needed, that are consistent with experimental observations. The development of such a model is an iterative process of (1.) model design based on existing knowledge, (2.) simulation and model-analysis, which results in (3.) the generation of new hypothesis that can be proven by experiments in the wet lab and used anew for model-refinement. This hypothesis driven approach based on in-silico experiments will support the experimental design or help to investigate questions that are not accessible to experimental inquiry.

The computation of time courses of a biochemical reaction system based on a given pathway structure and its kinetic scheme, that is required for simulations, has already been reviewed by Garfmkel et al. (1970). It arises from fundamental research on biochemical reaction kinetics (e.g. Michaelis and Menten 1913). The first simulation of a biochemical system for the peroxidase reaction was done by Chance (1943), who used a mechanical differential analyzer to solve his mathematical equations.

Cellular systems involve various components and complex interaction networks. To successful develop models for these systems computational tools for the model design, simulation and analysis are required. A brief overview of current software platforms and projects that face up to this as well as an overview about computational requirements for this purpose is given in Takahashi et al. (2002). Common systems — among others — for this purpose are Gepasi (Mendes 1993, 1997, Mendes and Kell 1998), E-Cell (Tomita et al. 1999, Takahashi et al. 2003), ProMoT/Diva (Ginkel et al. 2003), Virtual Cell (Schaff et al. 1997, Loew and Schaff 2001, Slepchenko et al. 2003), or the Systems Biology Workbench (SBW) and its add-ons (Hucka et al. 2002).

Gepasi comes up with a user-friendly interface for the simulation and analysis of biochemical systems. It supports the definition of compartments. Common kinetic types as well as user defined kinetic types are available. It provides time course and steady state simulation and the ability to explore the behavior of the model over a wide range of parameter values using a parameter scan, that runs one simulation for each parameter combination. Gepasi characterizes steady states using Metabolic Control Analysis (MCA) and linear stability analysis and is capable of doing parameter estimation with experimental data.

E-Cell is based on the modeling theory of the object-oriented Substance-Reactor Model. Models are constructed with three object classes, Substance, Reactor, and System. Substances represent state-variables, Reactors describe operations on state variables and Systems represent logical or physical Compartments. Version 2.25 contains 18 different classes of standard Reactors (e.g. Michaelis-Menten formula). Time course calculation is done by the use of a simulation engine. Numerical integration is supported by first-order Euler or fourth-order Runge-Kutta method.

ProMoT/Diva consists of the modeling tool ProMoT and the simulation environment Diva. The workbench deals with modular models and can handle Differential Algebraic Equation (DAE) systems. Modeling is supported by a graphical user interface and a modeling language. The modeling tool provides the possibility to use ready made modeling entities out of knowledge-bases. The authors claim, that their system is well suited for simulation and parameter estimation problems with up to 10000 differential equations.

The Virtual Cell is a web-based client-server architecture with a central database of user models. It provides a formal framework for modeling biochemical, electrophysiological, and transport phenomena while considering the subcellular localization of the molecules that take part in them (Slepchenko et al. 2003).

The SBW provides a server that acts as a broker between different modeling and simulation tools (clients) via a common interface. These clients (add-ons) cover graphical tools for model population, deterministic and stochastic simulators and analysis tools like the integration of MetaTool (Pfeiffer et al. 1999). Closely related to the SBW is the development of the Systems Biology Markup Language (SBML, Hucka et al. 2003) — a standard to store and exchange computer based biological models on the basis of XML (Exensible Markup Language). A similar approach is defined by CellML ( Lloyd et al., 2004).

SBML comes up with an open model repository, which comprises models for signal transduction, metabolics, and others. Another model repository is the JWS Online Cellular Systems Modeling software, that supports online simulation via a Java based interface and is maintained by B. Olivier and J. Snoep, Stellenbosch University and the Free University of Amsterdam (Snoep and Olivier, 2002, http : / / j j j . biochem . sun . ac . za/).

Several of these tools mentioned above use equation systems based on rate laws for the simulation. This widely accepted approach is often used for the simulation of metabolic networks but also for signal transduction (Kholodenko 2000, Levchenko et al. 2000), gene expression (Bhartiya et al. 2003) or other biological processes (e.g. Goldbeter 2002). None of the above mentioned tools supports direct access to metabolic or genomic databases or information resources about signal transduction or functional genomics.

The object of the present invention is to provide an improved method for predicting a dynamic behavior and, in particular, the time courses of a biological system. It is a further object of the present invention to provide a method for the quantitative and qualitative prediction of the dependencies of the biological objects included in said biological system. The method should provide a computer-executable model that is consistent with experimental observations. The method should further be adapted to analyze large models.

These objects are achieved with the features of the claims.

According to the present invention, the method for predicting the dynamic behavior, i.e., the spatial and temporal alterations, of a biological system and/or the quantitative and qualitative prediction of the dependencies of the biological objects included in said biological system, i.e., the relationship between objects in a biological system, is provided. According to the method of the present invention, a computational model of the biological system is generated. Said model comprises one or more model objects and associated model reactions, said model reactions being provided with kinetic parameters. The computational model may be provided as a hierarchically object oriented model. The method further comprises populating said model based on publicly available information using automated imported/export-functions, for example via a web-interface, via python scripts, or via SBML level 1. The model is further populated based on experimental data with respect to initial conditions and kinetic parameters of the model and its model topology. The term "model topology" generally describes the structure of a computational model of a biological system comprising objects, the network interconnecting the objects and the stoichiometry. A differential equation system, which may be an ordinary differential equation system or a differential algebraic equation system, describing the time dependency of the concentrations of at least one of the model objects is generated. The differential equation system is solved and thereby, said prediction of the dynamic behavior and particularly the time course of biological system and/or the quantitative and qualitative prediction of the dependencies of the biological objects are obtained.

The present invention further provides a computer program for performing the method of the present invention.

PyBioS, a computer program adapted to carry out the method according to the present invention has been developed, PyBioS being an object oriented environment for the design and simulation of biological models. This environment serves as a hierarchical object oriented database to store models of cellular systems. Each model is represented by a SimulationEnvironment, that comprises the model objects in a hierarchical object oriented manner. Model objects are entities of the abstract BioObject-class that represents biological objects. Derived of this class are concrete classes for biological entities that are subdivided into container-like objects (Environments), that can contain other BioObjects, and non container-like objects. This hierarchical structure is illustrated in Figure 4 A. Container-like object classes are Cell, Compartment, Complex and Clrromosome. Non container-like object classes are Gene, pre-mRNA, mRNA, Polypeptide, Protein and Enzyme (of which also Polymerase, Spliceosome, RNase, Ribosome and Protease are derived from). Additional information such as annotation, sequence-data, parameters and initial concentrations are stored as properties to an object. Actions, which describe reactions between different objects, are attached to BioObjects, e.g. a metabolic reaction is bound to its catalyzing enzyme.

Certainly, small subsystems can be modeled and analyzed to some extent in isolation by assuming steady and simplified boundary conditions. But as soon as these boundary conditions become variable — as given for a complex cell — it is clear that also this subsystem might behave differently in the context of a more comprehensive model. For instance ATP, one of the most important energy sources in the cell, is involved in diverse cellular processes and for example a massive consumption of ATP by a single process might have an important impact on other processes; this impact will not be discovered as long as the ATP concentration is handled as a constant parameter or variable of the isolated subsystem.

Thus, PyBioS is particularly developed for the analysis of large models. Here, automated import/export functions to populate models and design ODE systems are essential, h the following, an overview of the PyBioS system is given. The features are described using a published well-characterized metabolic process — glycolysis, published by Hynne et al. (2001). Population, simulation and analysis of the model is shown. It is a reasonably large model mvolving 36 reactions (forward- and backward-reactions are handled here as separate reactions) and 22 variables (metabolites). It can be accessed from the PyBioS model repository.

Automated generation of ODE systems, automated interface to the KEGG database (Goto et al. 1998, Kanehisa and Goto 2000) and scaling of the simulation system are demonstrated by exemplary networks that are populated via the interface to the KEGG database.

In the following, the invention is described in more detail with reference to the attached Figures, wherein

Fig.l illustrates cellular processes form a complex net of interacting objects,

Fig. 2 shows a graphical representation of hexokinase action,

Fig. 3 shows a simplified class diagram of the main PyBioS object-classes that represent biological objects,

Fig. 4A shows the hierarchical structure of the glycolysis model, Fig 4B a listing of all model reactions, and Fig 4C shows the view of the hexokinase action.

Fig. 4D: (1): Model repository. (2): Hierarchical structure of the selected model. (3): Automatically generated graphical representation of a model. (4): Listing of model reactions. (5): Graphical representation of simulation results. (6): User friendly representation of kinetic laws. Fig.s 4E and 4F: General database interface to biological databases E: Search interface. F: From the list of found reactions, some of them or several can be selected for population.

Fig. 4G: Parsing of the individual database information and automatic generation of a respective model.

Fig. 5 shows the network diagram of the glycolysis model from Hynne et al. (2001) generated by PyBioS,

Fig. 6 shows an artificial model; the artificial model shown in A illustrates the parameter- scanning functionalities of PyBioS, in B parameter k_t is scanned along the interval [0,10] and the steady-state concentrations are plotted versus their respective A,. -value, in C parameter k₂ is scanned along the interval [0.5,10] and the steady-state fluxes are plotted versus their respective k₂ -value, and

Fig. 7 shows the scaling behavior for systems in different sizes.

PyBioS provides functionalities for model design and population, simulation and analysis via a web interface. It has the following features, which will be discussed in this section in detail:

1. Model population: creation of model objects and reactions • python classes for biological entities are used for the formulation of hierarchical object-oriented cellular models • models are stored in an object-oriented database (ZOPE Object Database) • model population is done via the web-interface, python scripts or SBML • objects, that are involved in a reaction (enzyme, substrate, product and modifier lists), its kinetic law and further data and annotations stored as object properties are described • general database interface in order to exploit database information and to structure this information in a way that it can be fed into kinetic models of any size, for example a database interface to KEGG or Reactome • interface(s) to experimental functional genomics data

2. Visualization of the model-network

3. Automatic generation of the deterministic mathematical model (ODE system) 4. Simulation using numerical integrators

5. Methods for model analysis, computation of conservation relations, detection of steady-states, stability analysis, parameter scanning 6. Repository of models and kinetics

Model structure

Objects PyBioS employs an object-oriented strategy, that was firstly introduced by Stoffers et al. (1992). The autliors used classes for metabolic entities and biochemical reactions for the modeling and simulation of metabolic systems.

Models in PyBioS have hierarchical object-oriented structures. Each model is stored in a separate SimulationEnvironment object that contains the objects which represents the biological entities. Biological entities, like gene, mRNA, protein, compound, enzyme, complex or compartment, are derived from the same class BioObject. BioObjects might have different properties. For instance, if a Michaelis-Menten kinetic is used, parameters like K_M or V_τrax are properties of an Enzyme instance. Properties have a value, e.g. a floating point number. Properties might also be annotations, sequence-data, etc. Furthermore, one or several actions can be bound to a BioObject. An action describes a biochemical reaction or physical process or a group of similar reactions or processes. For instance a translation- action which is bound to a ribosome can be defined that operates on all mRNAs and produces the corresponding polypeptides. Actions are described in more detail below. Figure 3 gives an overview of the defined classes of biological objects and the information that is expected to be stored by these objects. Figure 3 shows a simplified UML class diagram of the main PyBioS object-classes that represent biological objects. In UML notation arrows point on classes from which other classes are derived. Central classes in the current version of PyBioS are the abstract classes BioObject and Environment. All classes which represent biological objects are derived from BioObject. BioObject has properties (in object-oriented programming denoted attributes) and methods (the diagram shows only some attribute and method examples). Properties are for instance the name of a BioObject (id) or its initial concentration (concentration). Methods are functions which belong to a certain class and operate on its attributes or other objects, that are passed along by the method call; e.g. getld() or getActions() return the object name or a list of the object's actions, respectively. Concrete classes, like Gene, Enzyme or Cell inherit from BioObject (and possibly other classes as well), which means that they have the same attributes (but likely other values) and offer the same methods as their parent classe(s). In the same way concrete classes can posses further methods. These methods might implement specific biological functions and thus they can be used during model population; e.g. a Ribosome has a translate method that generates for each mRNA object in the model (in this context the object-oriented model is given by a SimulationEnvironment instance) an appropriate Polypeptide instance, thus this methods reflects the real behavior of ribosome in a cell, and similarly other biological functions are implemented.

Some of these BioObjects, like Compartment, Complex or Chromosome, are container-like objects (derived from the abstract class Environment) which can hold other BioObjects.

Actions

One or several actions can be assigned to a BioObject. Each action holds a directed reaction and a kinetic law. Table 1 shows an exemplary action data structure. E, S, P, and Mare lists of size k, I, m, n, respectively. These lists reference BioObjects, that are involved in the reaction and their stoichiometric coefficients. Table 1:

Action attributes Representation of . . . Example Name (id): Reaction identifier e.g. Phosphorylation Enzyme (E): E . . . E_k e.g. Hexokinase Substrates (S): i_\ Sι . . . iι Si ATP, Glucose Products (P): jl Pi ■ ■ ■ jm Pm ADP, Glucose 6-phosphate Modifiers (M): i . . . M_n Reaction: k Si -ϊ ii Si ATP + Glucose → h Pi + ^{■ ■} " T" Jm -tm — r ADP + Glucose 6-phosphate

Kinetics: v = f(S) ., _ V_m [ATP][GIcl ^υ - Km)_cK_m.+K_Gk[AΥP]+K_AΥP[Glc]+[Glc\[K P}

The directed reaction describes the mass flow from the substrate(s) to the product(s), as well as the molecularity (stoichiometry) with which they take part in the reaction.. Substrates and

Products (denoted S and P respectively) as well as the catalyzing enzyme (denoted E) are stored in lists to the action. The elements of these lists are references to BioObjects together with their respective stoichiometric coefficients. Further lists for modifiers etc. can be added, if necessary. The term "modifier" defines a class of objects corresponding to compounds capable of interfering with the modeled reaction. Modifiers include activiators and inhibitors of an enzyme. Reversible reactions are either constructed by defining the backward reaction as a separate action or by using a rate law that already considers this behavior and thus might become negative. In the latter case reversibility is indicated by a flag, which is an attribute of the used rate law. Kinetic laws are formulated in an abstract fashion by a list that consists of parameter and variable references and fundamental mathematical operations (+, -, *, /, log, exp, ...) and parentheses. The final kinetic law-term is constructed from the lists of substrates, products, the enzyme and other modifiers. This makes it possible to come up with a database of predefined kinetic laws — that is available in PyBioS — of which kinetics can be chosen during the model-population as desired. There is no rule that lay down to which BioObject an action should be assigned to but for instance for enzymatic reactions it makes sense to attach it to the enzyme that catalyses the given reaction. Autocatalytic reactions can be assigned to the substrate itself and chemical reactions that take place in the absence of any catalyst can be assigned to the compartment-object to which the reactants belong to or to a pseudo-object whose only task then is to represent this action. An action that describes a transport process can be bound to a membrane or transporter complex-instance. Locally defined actions have in contrast to global action list the advantage, that similar reactions in different compartments don't have to be denoted differently. This is very useful for large models.

Model population

The first step of modeling is the collection of objects and reactions and appropriate kinetics that are relevant to the model. Using these data a draft can be populated. PyBioS provides three methods for model population: (1.) via the web-interface, (2.) via python-scripts, or (3.) via SBML Level 1.

In the following the creation of a new model via the web-interface is described. As an example the hexokinase reaction of glycolysis is used (Figure 2). First, an appropriate web- page is entered, e.g. the PyBioS homepage, using a web-browser and a respective model- folder on said homepage is selected. A list of models that are available in this folder is then displayed. To create a new model, a model identifier and optionally a short title for the model is requested. A new SimulationEnvironment is thus created. This empty simulation environment is the place for a single hierarchically ordered cellular model (cf. the object tree of the glycolysis model in Fig. 4A). New BioObjects can be added, for instance one can select "Cell" from a menu and press the create button to create a new Cell-instance. Then, an identifier for the new cell may be entered and the "Cytoplasm" from the compartment hierarchy can be selected. This will create a new cell-object and the Compartment-object cytoplasm. The cytoplasm should be entered by following its link and then "Enzyme" from the object menu should be selected. PyBioS will search all possible reactions for a certain EC-number that has been entered after selecting "Enzyme", and present the results. When, e.g., the reaction ATP + alpha-D-Glucose -> ADP + alpha-D-Glucose 6- phosphate is selected, the enzyme hexokinase and the substrate and product compound objects will be created. The link of the new enzyme hexokinase directs to the View-tab of this object. The Actions-tab shows a list of all actions, that are bound to this object, and each action identifier directs to a detailed action interface (Fig. 4C). This interface has four frames: The upper left one shows the Enzyme-, Substrate-, Product- and optionally any modifier-lists. Each list holds references to the involved BioObjects and their corresponding stoichiometric coefficients. The upper right-hand frame enables the user to add BioObject-references from the objects available in model tree. The bottom left frame offers the possibility to choose a kinetic from a list of predefined laws, hi case of requiring a kinetic-law that is not available in the kinetics database, the user can simply define an individual one (cf. description of Action).

Figure 4D shows another view of the PyBioS user interface. A particular model can be inspected via different views that are accessible by several tabs. These views provide a representation of the hierarchical model structure, a listing of model reactions, a graphical representation of the model network, user interfaces for simulation and analysis and other functionalities. The population-tab provides the access to the generic database interface and supports an easy model design. The user can search for reactions a specific gene or metabolite is involved in, or simply for a whole pathway (Figure 4E). From the list of results the user can select several or all reactions for model population (Figure 4F). The list of reactions that are dedicated for population can be extended step by step with reactions of further database searches. All objects and reactions populated by the generic database interface still provide references to their originating database entry. This feature enables an automated annotation of all model components and an easily available model extension by further database requests. Additionally, these database references make it easier to merge models with other models and thus it support the integration and re-usability of well defined models.

As seen by these examples, interfaces to the metabolic data of the KEGG-database and the pathway information available from the Reactome database (Josbi-Tope et al. 2005, ) enable an easy population of the model with biochemical reactions. Using this, PyBioS is able to take care of common compound names as they are used in biochemistry.

An interface to the Systems Biology Markup Language (SBML level 1; Hucka et al. 2003) enables model-exchange with other systems. The PyBioS database contains several models that are available on the SBML repository page (see, e.g., http : / /www. sbw-sb l . org/ModelsWebPages/ModelRepository . htm). This repository comprises models for the cell-cycle, metabolic and circadian oscillators, MAPK cascade signal transduction, glycolysis and others.

Models of cellular systems presently scale up to about 30 reactions and as much reactants (e.g. Hynne et al. 2001), but even systems with up to 127 genes and their corresponding reactions have been implemented (Tomita et al. 1999). Since the population of large models is labor intensive the availability of an appropriate application programmer interface for this task is helpful. Therefore Pybios makes it feasible to populate models by python-scripts using PyBioS' programming interface.

In a preferred embodiment, step (b) of the method of the invention is effected using a general database interface which permits exploitation of database information and structurmg of this information such that it can be fed into kinetic models of any size.

In a preferred embodiment, step (c) of the method of the invention is effected using interface(s) to experimental functional genomics data.

The purpose of both features is to gain information on biological processes and diseases via complex mathematical model that use comprehensive information on all levels of cellular information (genome, proteome, interactome, metabolome etc.). Models are hybrid in the sense that they are structured in modules of defined input and output. The individual model underlying a module can be of different mathematical nature (ODE, PDE, stochastic, qualitative etc.) and can have a flexible input/output structure.

Functional genomics data encompass DNA sequence data, gene expression data based on quantitative (DNA arrays, PCR, ChJJ?) and qualitative data (SAGE, WISH, CGH), proteomics data (mass spectrometry, ICAT, protein interaction data, 2D gels), kinetic measurements (Biacore) and measurements and information based on particular phenotypes, environmental, clinical and physiological data and provide initial conditions for the simulation run and/or may be employed for validation purposes, e.g. comparison of simulated time-dependent bahaviour or steady-state concentrations to experimental data. General database interface in order to exploit database information

A new database interface to biological databases has been developed. Many of current knowledge remains unexploited for kinetic modelling since it cannot be adequately translated from the database to the modelling system. PyBioS supports this exploitation by 1. the definition of a unique model ontology 2. the development of parser for the access of relevant database information 3. the development of a user interface that connects the modelling system with the databases the development of an applied programming interface that enables the user to import the database information into the computational models.

Figures 4E and 4F show examples of the database interface. The user can specify the database and query that database with search terms. The example shown here relates to the Reactome database at EBI.

hi the next step the user can parse the individual database information (here reactions) and select what should be included in the mathematical PyBioS model (see Figure 4G). The information is then translated into a model directly.

Interface(s) to experimental functional genomics data

Large networks for several organisms have been generated. These networks are metabolic networks, signalling networks, gene regulatory networks and physiological networks. In order to compare two different states of the network (for example a normal and a disease state) we incorporate the experimental data obtained from biological material of control and disease groups and direct these data into the models. A pre-processing performs data normalisation, statistical data analysis in order to eliminate technical bias. By comparing two large-scale hybrid models we are able to draw qualitative predictions on different biological information levels.

The different steps of this strategy are: 1. Experimental design The definition of the experimental questions and techniques that will be applied combined with a strategy for experimental design and methods for the pre-processing of data. For example, RNA can be extracted form a disease state and a normal state of the organism. DNA arrays are hybridised with labelled RNA and microarray data is analysed by specific statistical methods. Results are lists of signal values for ~700 enzymes active in the >110 processes in the human metabolic network 2. Definition of a simulation set up The second step consists of a strategy for the simulation. This involves defining the parameters of the model and methods for the estimation of these parameters. Here, we incorporate knowledge from known reaction kinetics, dosage effects, etc. Alternatively or additionally, in particular in those cases where kinetic data are not available, said kinetic data may be repeatedly generated randomly. Subsequently, several simulations runs, e.g. about 10 to about 100, using the different randomly generated kinetic data sets may be performed. 3. Prediction The third step yields the model prediction. For example, the model predicts the concentration curves of the metabolites of the metabolic network. By comparing both conditions (normal and disease state) we derive candidates for the subsequent metabolomics experimental work. This strategy reduces the cost factor tremendously since it allows a directed approach to metabolic targets and possible drugs.

Visualization

PyBioS model networks are given by the BioObjects and their actions. Since already a model of even 10 or 20 BioObjects easily becomes complex by diverse substrates, products and modifiers, a visualization of the underlying model structure is of substantial benefit. Therefore, the biological network can be made more easily accessible by a graphical representation. PyBioS uses a representation in which nodes are connected via arrows. Figure 5 shows a network diagram of the glycolysis model from Hynne et al. (2001) generated by PyBioS using the GraphViz library (see http : / /www , research , att . co / sw/ tools /crraphviz). The cut-out at the right side shows the forward and reverse reaction of the isomerization of glucose-phosphate to fructose-phosphate. Line styles differ between mass- and information-flow. Two kinds of nodes, which stands for either BioObjects (visualized by rectangles) or actions (visualized by ellipses), are alternating connected via directed arrows. The arrow color or style indicates either mass- (black) or information- flow (any other color or line-style). The direction of the mass-flow-arrow indicates the mass-flow. Information-flow-arrows always point from BioObject- to Action-nodes, since they represent BioObjects which catalyze or modify the particular action.

Quantitative simulation

Stoichiometric equations and rate laws that are given by the actions can be used for the automatic generation of ordinary differential equation (ODE) systems. The time change in the concentration of species is given by the following balance equation:

dS_i /dt = ∑n_iJv .(S) ι = l,. ,m (1) =ι

Where r is the number of reactions, m the nmnber of species, and n_f] the stoichiometric coefficient of S. in the reaction j, which is positive for the production of S_t , negative for its degradation and otherwise defaults to zero. v_j denotes the velocity of reaction y, that is given by its rate law.

PyBioS supports deterministic simulation by numerical integration of first order ODE- systems. It offers the use of the solvers LIMEX and LSODA (denoted SciPy in the modeling interface) to get the numerical solution of the initial value problem. LSODA (Hindmarsh 1983, Petzold 1983) is a solver for ordinary differential equations written in Fortran and a variant of the LSODE package. The algorithm used in this solver switches between stiff and nonstiff methods, automatically.

PyBioS uses the interface to LSODA which is available from SciPy

(http : / /www . scipy . org) in a similar way as described by Olivier et al. (2002). The solver LIMEX (Deuflhard and Nowak 1987) is an extrapolation integrator for the solution of linearly-implicit differential-algebraic systems (DAEs) written in Fortran (ftp : / /elib . zib . de/pub/elib/codelib/LIMEX4 2A1 , Deuflhard and

Nowak, 1987). It combines an implicit one step method with stepsize extrapolation to permit an adaptive control of stepsize and order.

For LIMEX the ODE system is converted into an appropriate Fortran source code that is compiled. Pybios then starts the executable program and parses the caculated solution from the output. The simulation results are reported in equidistant dense points, (in the current version 800 dense points are used). Time course concentration and time course fluxes are available as text and as graphic. The user has the possibility to choose a representation of the results in an HTML page or to get the results by download. Graphical diagrams of the results prepared by interfaces to Gnuplot or Dislin (http : / /www . gnuplot . info or http : / /www. linmpi . mpg . de/dislin) or simple tables in text form are available.

Steady state search

When the system has reached a steady state the concentration of the metabolites does not change any longer by time. The steady state of a system of reactions is characterized by ,_.(.is

— = O..S defined in equation 1). In a nonlinear system this equation can have several dt solutions. As PyBioS creates a system of autonomous differential equations, stability analysis is performed by linarizing the mostly nonlinear system at the steady state. The analysis method computes the eigenvalues of the Jacobian matrix of the ODE system at the steady state. From the eigenvalues information about stability and the structure of the phase plot can be derived. A stable steady state is given, if the real-parts of all eigenvalues are below zero (Walter 1998). Two different methods for steady state search are available: First an approach that depends on a root finding method to get the steady state by numerical algorithm and second a progressive simulation, that we denote the direct search. Starting with simulation results of a user-defined time interval, the root finding approach computes the roots of the ODE system by using the MINPACK subroutine HYBRID 1

(http : / /www. netlib . org/minpack), which is a modification of the method described by Powell (1970). The direct search performs a progressive time course calculation: Starting with the time interval [t₀ , t_n ] given by the user a series S_t , • • • , S_t is calculated, where

is a vector with the concentrations of all species at time t_i . The direct search algorithm checks whether a steady state is reached by regarding

S. +J+l - S. < ε j = ,-,k-l (3)

for a user defined thresholds . Here, ||- • denotes tl e Euclidean norm. If this equation is satisfied for twenty consecutive time points in the interval [t„,t„_+it], it is assumed, that the steady state is reached. Otherwise another evaluation step starts for the interval [t, +k ' ^ n+k

This is repeated up to / times (in the current version / = 10) until the steady state is found. In case of an unsuccessful search, the algorithm aborts and reports this.

The use of the direct search is recommended, if the Jacobian Matrix of the ODE system is singular and the root-finding might fail. By the threshold value and the size of the simulation time interval the accuracy and performance of the 'direct search'- algorithm can be controlled. For instance, if the algorithm fails for a certain parameter value, a decrease of the threshold and/or an increase of the time interval enhances the accuracy but with a decline in performance.

The Jacobian matrix is computed by forward difference formulas. PyBioS calculates the eigenvalues of this matrix by using the LinearAlgebra module of NumericalPython (http : / /www .pfdubois . co /numpy/) which contains an interface to the according LAPACK functions (Anderson et al. 1999). Eigenvalues of the Jacobian matrix, steady state concentrations and steady state fluxes are displayed by an HTML table. Analysis modules

The object-oriented model of the biological system as well as its derived mathematical model can be used for further analyses and consistency checks. PyBioS offers the possibility to compute conservation relations and perform parameter scans.

Stoichiometric Matrix and Conservation Relations

Frequently, the quantity of material of several species involved in cellular reaction networks is conserved, e.g. n(AT?) + n(ADP) + n(AMP) - const. (4)

Where n() denotes the amount of substance.

Such conservation relations can be computed by the network topology, that is given by the reactions and their stoichiometry. This topology of the reaction network describes the mass flow and is embodied in the stoichiometric matrix:

A column of the matrix corresponds to a distinct action in the network, a row corresponds to a single biological species (BioObject). r is the number of actions and m the number of species. An element n_tj indicates whether a certain BioObject takes part in a particular reaction or not. The conservation matrix Tcan be obtained by computing the nullspace (kernel) of the stoichiometric matrix Ν using the relation TΝ = 0.

Since it is more conventional to compute the right nullspace, T^τ (where T denotes the transposed matrix) is calculated from N^ΓΓ^Γ= O (6)

using the block diagonalization algorithm described by Schuster and Schuster (1991).

Parameter scan

Parameter scan can be performed to analyze the behavior of the model. The possibility to consider the effect of one parameter on the concentrations of the metabolites and on the fluxes of the reactions is given by regarding the system in the steady state. In the steady state the system is independent of time and an implicit dependence of the concentrations and fluxes on a parameter can be viewed. One parameter is varied in a given interval and the according steady states are computed by the direct search or root finding method. PyBioS offers two possibilities: Representation of either the fluxes or the concentrations in steady state for the selected parameter in the given interval. The results of parameter scan are available in the same various manner like the simulation results and the stability analysis is available in text form.

For the analysis of characteristics of biological interest stability analysis (as described above) over the given parameter range can be done. Here, in some cases the user has also the chance to draw conclusions about the possible dynamic behavior of the biological system.

The parameter scanning is illustrated by the artificial model shown in Figure 6A. A scanning for parameter k_t in the interval [0,10] indicates that the concentration of S s independent and S₀ and S₂ are dependent of this parameter (Figure 6B). This is confirmed by its analytical solution. Similarly, flux changes can be analyzed. In Figure 6C, the parameter k₂ is scanned along the interval [0.5, 10] and the steady-state fluxes are plotted versus their respective k₂ -value. Both scans were performed at 46 points in the given interval by using the root finding method. Non-varied parameters and initial-values are set to 1 by default. Below the diagrams the respective rate-laws, ODEs and analytical solutions for the concentrations and fluxes at steady-state are shown. Test of systems performance

Since molecular biological data becomes massively available through the Internet and by enormously evolving high-throuput techniques, strategies and methods for the integration of these data into biological models are required. Small systems of 20 objects or less can directly be translated into mathematical models by hand. But the population of models with several dozens, hundreds or even thousands of objects become no more feasible without an automation of this process. Therefore the huge amount of experimental data as well as textbook data — which becomes more and more available in a computational amenable manner — are excellent sources for this purpose. PyBioS supports functionalities for the integration of external data sources. An interface to the metabolic data of the KEGG database enables the automated population of a model with a single or several pathways. Since the model population is one central step in the process of model design, its scaling behavior is of interest. Therefore, metabolic models of different sizes in the number of reactions and objects were populated, exemplarily, using the interface to the KEGG database. This ranges from models with 20 objects and 11 reactions upto 1668 objects and 2365 reactions. For simplification all reactions take place in the same compartment and are modeled by mass action rate laws. Kinetic parameters and initial concentrations are initialized with 1. The CPU time required for this population process was measured. Figure 7A shows that the duration of the population process scales linearly with the model size in this example with metabolic systems derived from the KEGG database. In parallel, the duration required for the simulation of the time-interval [0,10] and [0,1000] (arbitrary units) using the SciPy-solver was measured for each model. Therefore, a quadratic relation of time versus model size (given by the number of reactions) was found (Figure 7B). The inserted graphic in Figure 7B shows the scaling behavior of some models from the PyBioS models repository: (A) CellCycle-1991Tys-2, (B) CellCycle-1991Gol, (C) CellCycle-1991Tys, (D) lVLAPKcasc-2000Kho, (E) CircClock-2002Vil, (F) Metabolism-2000Teu, (G) CircClock- 1999Lel, (H) CellCycle-1997Nov, (I) Hynne; (A)-(H) are imported via SBML from the SBML-model repository; (I) is described in Hynne et al. (2001).It should be noted, that the simulation time depends strongly on the complexity of the used kinetic laws. The scaling behavior of some published models is also illustrated in Figure 7B. Summary

PyBioS is an object-oriented framework for the modeling, simulation and analysis of cellular systems. It comes up with repositories for cellular models and predefined kinetic laws.

The integration of interfaces to external databases offers functionalities for the automated population of large models as it has been shown by using the KEGG database. The definition of certain methods for special object-classes further extends the model population strategies. These methods enable the implementation of general models for transcription, translation or enzymatic reactions and they can acquire further object or reaction data from the available databases. The system can be extended by interfaces for the integration of high-throughput- data.

The automatic generation of mathematical models from the object-oriented models and the subsequent functionalities for simulation and graphical representation of the results are the basis for the hypothesis driven approach. Performing simulations with PyBioS permits the formulation of hypotheses and provides guidance and decision support for the design of wet- lab experiments, e.g. for knock-out or RNAi experiments. Furthermore, PyBioS is not only able to support ODE-systems, but also DAE-systems.

Computation of conservation relations and the parameter scanning methods offer some fundamental functionalities for in-silico experimentation. These can be extended by further analysis-modules:

• Bifurcation analysis

• Metabolic control analysis (MCA)

• Simplification of models, e.g. by using knowledge about conserved moieties

• Optimization and parameter fitting

• Computation of elementary modes • Generation and simulation of stochastic models

• Implementation of support for system parameters other than species concentration to be handled as dependent variables (e.g. volume of a compartment, tungor-pressure, pH, or electrophysiological quantities). References

Anderson E., Bai, Z., Bischof, C.,Blackford, S., Demmel, J., Dongarra, J.,Du Croz, J., Greenbaum, A., Hammarling, S., McKenney, A. and Sorensen, D. 1999. LAPACK users' guide. Society for Industrial and Applied Mathematics, Philadelphia, ISBN 0-89871-447-8.

Bhartiya S., Rawool S., and Venkatesh K.V. 2003. Dynamic model of Escherichia coli tryptophan operon shows an optimal structural design. Eur. J. Biochem. 270: 2644-2651. Chance B. 1943. The kinetics of the enzyme-substrate compound of peroxidase. J. Biol. Chem. 151: 553-577.

Deuflhard P. and Nowak U. 1987. Extrapolation integrators for quasilinear implicit ODEs. In Prog. Sci. Comp. 7, eds. P. Deuflhard and B. Engquist.

Garfinkel D., Garfinkel L., Pring M., Green S., and Chance B. 1970. Computer Applications to Biochemical Kinetics. Annual Review of Biochemistry 39: 473-498.

Ginkel M., Kremling A., Nutsch T., Rehner R., and Gilles E. 2003. Modular modeling of cellular systems with ProMot/Diva. Bioinformatics 19: 1169-1176.

Goldbeter A. 2002. Computational approaches to cellular rhythms. Nature 420: 238- 245.

Goto S., Nishioka T., and Kanehisa M. 1998. LIGAND: chemical database for enzyme reactions. Bioinformatics 14: 591—9.

Hindmarsh A. 1983. Odepack, a systematized collection of ode solvers. In Scientific computing, eds. R. Stepleman and et al. North-Holland, Amsterdam, pp. 55-64. Hucka M., Finney A., Sauro H.M., Bolouri H., Doyle J., and Kitano H. 2002. The ERATO Systems Biology Workbench: enabling interaction and exchange between software tools for computational biology. Pac Symp Biocomput : 450-61.

Hucka M., Fimiey A., Sauro H.M., Bolouri H., Doyle J.C., Kitano H., Arkin A.P., Bomstein B.J., Bray D., Comish-Bowden A., Cuellar A.A., Dronov S., Gilles E.D., Ginkel M., Gor V., Goryanin I., Hedley W.J., Hodgman T.C., Hofmeyr J.H., Hunter P.J., Juty N.S., Kasberger J.L., Kremling A., Kurnmer U., Le Novere N., Loew L.M., Lucio D., Mendes P., Minch E., Mjolsness E.D., Nakayama Y., Nelson M.R., Nielsen P.F., Sakurada T., Schaff J.C., Shapiro B.E., Shimizu T.S., Spence H.D., Stelling J., Takahashi K., Tomita M.,Wagner J., andWang J. 2003. The systems biology markup language (SBML): a medium for representation and exchange of biochemical network models. Bioinformatics 19: 524-31.

Hynne F., Dano S., and Sorensen P.G. 2001. Full-scale model of glycolysis in Saccharomyces cerevisiae. Biophys Chem 94: 121-63.

Joshi-Tope G, Gillespie M, Vastrik I, D'Eustachio P, Schmidt E, de Bono B, Jassal B, Gopinath GR, Wu GR, Matthews L, Lewis S, Birney E, Stein L. 2005. Reactome: a knowledgebase of biological pathways. Nucleic Acids Res. 33 Database Issue:D428-32. Kanehisa M. and Goto S. 2000. KEGG: kyoto encyclopedia of genes and genomes. Nucleic Acids Res 28: 27-30. Kholodenko B.N. 2000. Negative feedback and ultrasensitivity can bring about oscillations in the mitogen-activated protein kinase cascades. Eur. J. Biochem. 267: 1583-1588. 21

Kitano H. 2002. Computational systems biology. Nature 420: 206-10. Levchenko A., Bruck J., and Stemberg P.W. 2000. Scaffold proteins may biphasically affect the levels of mitogen-activated protein kinase signaling and reduce its threshold properties. Proc. Natl. Acad. Sci. 97: 5818-5823.

Loew L.M. and Schaff J.C. 2001. The Virtual Cell: a software environment for computational cell biology. Trends Biotechnol 19: 401-6.

Lloyd CM., Halstead M.D., Nielsen P.F. 2004. CellML: its future, present and past. Prog. Biophys. Mol. Biol. 85: 433-50. Mendes P. 1993. GEPASI: a software package for modelling the dynamics, steady states and control of biochemical and other systems. Comput Appl Biosci 9: 563- 71.

Mendes P. 1997. Biochemistry by numbers: simulation of biochemical pathways with Gepasi 3. Trends Biochem Sci 22: 361-3.

Mendes P. and Kell D. 1998. Non-linear optimization of biochemical pathways: applications to metabolic engineering and parameter estimation. Bioinformatics 14: 869-83.

Michaelis L. and Menten M. 1913. Die Kinetic der Invertinwirkung. Biochemische Zeitschrift 49: 334-369.

Olivier B.G., Rohwer J.M., and Hofmeyr J.H. 2002. Modelling cellular processes with Python and Scipy. Mol Biol Rep 29: 249-54. Petzold L. 1983. Automatic selection of methods for solving stiff and nonstiff systems of ordinary differential equations, siamj. sci. stat. comput. 4: 136-148.

Pfeiffer T., Sanchez- Valdenebro I., Nuno J., Montero F., and Schuster S. 1999. METATOOL: For Studying Metabolic Networks. Bioinformatics 15: 251-257.

Powell M. 1970. A Hybrid Method for Nonlinear Equations. In Numerical Methods for Nonlinear Algebraic Equations, ed. P. Rabinowitz. Gorden and Breach, pp. 87-114.

Schaff J., Fink C.C., Slepchenko B., Carson J.H., and Loew L.M. 1997. A general computational framework for modeling cellular structure and function. Biophys J 73: 1135- 46.

Schuster S. and Schuster R. 1991. Detecting Strictly Detailed Balanced Subnetworks in Open Chemical Reaction Networks. J. Math. Chem. 6: 17-40. Slepchenko B.M., Schaff J.C., Macara I., and Loew L.M. 2003. Quantitative cell biology with the Virtual Cell. Trends Cell Biol 13: 570-6.

Snoep J.L. and Olivier B.G. 2002. Java Web Simulation (JWS); a web based database of kinetic models. Mol. Biol. Rep. 29: 259-63.

Staffers H.J., Sonnhammer E ., Blommestijn G.J., Raat N.J., and Westerhoff H.V. 1992. METASEVI: object-oriented modelling of cell regulation. Comput Appl Biosci 8: 443-9.

Takahashi K., Ishikawa N., Sadamoto Y., Sasamoto H., Ohta S., Shiozawa A., Miyoshi F., Naito Y., Nakayama Y., and Tomita M. 2003. E-Cell 2: Multi-platform E-Cell simulation system. Bioinformatics 19: 1727-1729. 22 Takahashi K., Yugi K., Hashimoto K., Yamada Y., Picket! C, and Tomita M. 2002. Computational challenges in cell simulation. IEEE Intelligent Systems 17: 64-71.

Tomita M., Hashimoto K., Takahashi K., Shimizu T.S., Matsuzaki Y., Miyoshi F., Saito K., Tanida S., Yugi K., Venter J.C., and Hutchison C. A. r. 1999. E-CELL: software environment for whole-cell simulation. Bioinformatics 15: 72-84.

Walter W. 1998. Ordinary Differential Equations, chap. Stability and Asymptotic Behavior. Springer Verlag, New York. Current international consortia including PyBioS

EU:

EMI-CD STREP, 01.01.2004-31.12.2006

EMBRACE Network of Excellence, 01.02.2005-31.01.2010

Germany: NGFN, SMP Protein, 01.11.2004-31.10.2007

BMBF, BioProfile Nutrigenomik, 01.03.2005-31.07.2007

Book Publication including PyBioS:

Klipp, E., Herwig, R., Kowald, A., Wierling, C, and Lehrach, H. (2005) Systems biology in practice. Wiley VCH Weinheim.

Claims

Claims:

1. A method for the prediction of the dynamic behavior of a biological system and/or the quantitative and qualitative prediction of interdependencies of the biological objects included in said biological system, said model comprising the steps: (a) generating a computational model of said biological system, said model comprising one or more model objects and associated model reactions, said model reactions being provided with kinetic parameters; (b) populating said model based on publicly available information using automated import/export functions; (c) populating said model based on experimental data with respect to initial conditions and kinetic parameters of said model and its model topology; (d) generating a differential equation system describing the time dependency of the concentrations of at least one of the model objects generated using the parameters provided in step (a); (e) solving the differential equation system, thereby obtaining said prediction of the dynamic behavior of the biological system and/or said quantitative and qualitative prediction of interdependencies of the biological objects contained in the biological system.

2. The method according to claim 1, further comprising the step: (f) analyzing the model system with respect to stoichiometry, conservation relations, and steady states, as well as performing parameter scans.

3. The method according to claim 1 or 2, further comprising the step: (g) refining the model system by comparing and optimizing the predictions computed in step (e) with experimental data.

4. The method according to claim 1, 2 or 3, wherein said computational model is a hierarchically object oriented model, comprising said model objects in a hierarchically object oriented manner.

5. The method according to claim 4, wherein said model objects are entities of an abstract class representing biological objects.

6. The method according to claim 5, wherein from said class representing biological objects, concrete classes for biological entities are derived.

7. The method according to any one of the preceding claims, wherein said model population is done via a web-interface, python scripts and/or SBML.

8. The method according to anyone of the preceding claims, wherein the step of generating a differential equation system comprises generating an ordinary differential equation system.

9. The method according to anyone of the preceding claims, wherein the step of generating a differential equation system comprises generating a differential algebraic equation system.

10. The method of anyone of the preceding claims, wherein the model is adapted to predict the reaction to a stimulus, a compound and/or a mutation.

11. The method according to anyone of the preceding claims, further comprising the step of visualizing the model structure.

12. Computer program comprising program code means for performing the method of anyone of the preceding claims when the program is run on a computer.