In Silico Biology 3, 0029 (2003); ©2003, Bioinformation Systems e.V.  
Special Issue: Petri Nets for Metabolic Networks

Quantitative Petri net model of gene regulated metabolic networks in the cell

Ming Chen* and Ralf Hofestaedt

Bioinformatics / Medical Informatics, Technische Fakultaet, Universitaet Bielefeld
Postfach 10 01 31, D-33501 Bielefeld

*  corresponding author

Edited by E. Wingender; received April 28, 2003; revised and accepted April 29, 2003; published June 07, 2003


A method to exploit hybrid Petri nets (HPN) for quantitatively modeling and simulating gene regulated metabolic networks is demonstrated. A global kinetic modeling strategy and Petri net modeling algorithm are applied to perform the bioprocess functioning and model analysis. With the model, the interrelations between pathway analysis and metabolic control mechanism are outlined. Diagrammatical results of the dynamics of metabolites are simulated and observed by implementing a HPN tool, Visual Object Net ++. An explanation of the observed behavior of the urea cycle is proposed to indicate possibilities for metabolic engineering and medical care. Finally, the perspective of Petri nets on modeling and simulation of metabolic networks is discussed.

Key words: metabolic network, gene regulation, Petri nets, quantitative model, urea cycle, modeling and simulation


With the success of the human genomic project, we have experienced increasing floods of data, both in terms of volumes and in terms of new databases and new types of data. More and more experimental data both on the genetic and cellular level are systematically collected and stored in specific databases that are also available to public via the Internet [Baxevanis, 2003]. Some well-known databases are: gene sequence (e. g. GenBank, EMBL, DDBJ), protein (e. g. SWISS-Prot, PIR, BRENDA), biochemical reactions (e. g. KEGG, WIT/MPW), transcription factors (TRANSFAC) and signal induction reactions (e. g. CSNDB, TRANSPATH, GeneNet). In the post-genomic era, the focus is now shifting to the so called "from the sequence to the function", i. e., in addition to completing genome sequences, we are learning about gene expression patterns and protein interactions on the genomic scale. Undoubtedly, analysis of metabolic networks is becoming a promising field, which requires providing new algorithms and tools to fulfill this task.

The study of gene regulated metabolic networks plays an important role in the detection of genetic/metabolic defects as well as drug research. Genetic/metabolic defects often lead to metabolic blockades, resulting in metabolic diseases. Many inborn errors of metabolism result from a single gene encoded enzyme deficiency. Regarding drug research, it is necessary to first understand the reaction pathways that are affected by the drug, directly and indirectly, and to know the effect of the modification of specific reaction steps on reaction networks. For instance, it is known that hyperammonemia is a hereditary disease concerning the urea cycle, resulting from an enzyme deficiency. Ornithine transcarbamylase (OTC), a deficiency of an X-linked enzyme disorder of urea synthesis, leads to a disease whose clinical manifestations include lethargy, coma, and cerebral edema. The identification of sites in a metabolic pathway that result in such diseases, caused by inborn errors in metabolism, would be useful so that the relevant enzyme or metabolite could be substituted or suitably modified. In addition, simulating the related complex metabolic networks will additionally help to understand the impact of various factors (e. g. enzyme insufficiency, metabolic blockade, drugs effects, etc.) on metabolic systems. This is particularly useful in the pharmaceutical industry for designing site-directed drugs to target mutant enzymes.

In order to understand the molecular logic of the cell, methods of modeling and simulation are of importance. Different models are available in the literature, commonly falling into two different categories: the descriptive models (discrete approaches) and the analytical models (differential equations). Traditionally, metabolic pathways are regarded as coherent sets of enzymatic reactions and can be interpreted as relational graphs. Each node represents a metabolite and each edge represents a biochemical reaction that is catalyzed by a specific enzyme. Kohn and Letzkus [1983] expanded the graph theory by a specific function that allowed the modeling of dynamic processes. Then, the first application of Petri net on modeling metabolic pathways was introduced by Reddy et al. [1993]. In contrast to naive graphs, Petri net is a graphically oriented language of design, specification, simulation and verification of systems. It offers a formal way to represent the structure of a discrete event system, simulate its behavior, and draw certain types of general conclusions on the properties of the system. Ordinary Petri net models do not have such functions as quantitative aspects, so there are some extension of Petri nets that can support dynamic change, task migration, superimposition of various levels of activities and the notion of mode of operations. Various extensions of PNs, such as (Stochastic) Timed PNs [Wang, 1998; Wang, 1999], Colored PNs [Kurt, 1997], Predicate/Transition Nets [Genrich, 1987] and Hybrid PN [David and Alla, 1992], allow for qualitative and/or quantitative analyses of resource utilization, effect of failures, and throughput rate. Hofestaedt and Thelen [1998] also presented an extension formalization, a self-modified Petri net, which allows the quantitative modeling of regulatory biochemical networks. During the last years some more papers appeared [Genrich et al., 2001; Goss and Peccoud, 1999; Hofestaedt et al., 2000a; Matsuno et al., 2000; Matsuno et al., 2001], indicating that the Petri net methodology seems to be useful in modeling and simulation of metabolism. We are motivated to exploit the methodology of Petri net to model gene regulated metabolic networks in the cell, explain the importance of sustaining core research, and identify promising opportunities for future research.

Hybrid Petri nets

We suppose that readers have some background knowledge of Petri nets, otherwise, they are strongly referred to read W. Reisig's Petri Nets: An Introduction [Reisig, 1985] or some basic reference books at Petri nets world at where a large amount of investigations on Petri nets have been compiled in the literature, and various applications have chosen Peti nets as their control models due to the intuitively understandable graphical notation of Petri nets. Herewith a brief description of hybrid Petri nets is presented as the following context.

Definition 1: A hybrid Petri net is a six tuple Q = (P, T, Pre, Post, h, M) such that:
P = {P1, P2,...,Pn} is a finite, not empty, set of places;

T = {T1, T2,...,Tm} is a finite, not empty, set of transitions;

P  T = , i. e. the sets P and T are disjointed;

h : P T {D, C} , called "hybrid function", indicates for every node whether it is a discrete node (sets PD and TD) or a continuous node (sets PC and TC);

Pre : P T + or , is the input incidence mapping (+ denotes the set of positive real numbers, including zero, and N denotes the set of natural numbers);

Post : P T + or is the output incidence mapping;

M : P + or is the initial marking.

We denote by M(t) = (m1t, m2t,...,mnt) the vector which associates with each place of P its marking at the instant t. M0 = M(t0) = (m10, m20,...,mn0) is the initial marking. At any time the present marking M is the sum of two markings Mr and Mn, where Mr is the reserved marking and Mn is the non-reserved marking. If h(Pi) = D or C then, mi(t) = mir(t) + min(t).

When a variable dTj (called the delay time of Tj) is assigned to each discrete transition Tj(h(Tj) = D), Tj is fired at time t + dTj

Pi ºTjTj denotes the set of input places of transition Tj),
mi(t) > Pre(Pi, Tj),
mi(t + dTj)  = mi(t)-Pre(Pi,Tj)

Pi Tjº (Tjº denotes the set of output places of transition Tj),
mi(t + dTj) = mi(t) + Post(Pi,Tj)

When a variable vTj (called the speed of Tj) is assigned to each continuous transition Tj(h(Tj) = C), Tj is fired at time t during a delay dt

Pi ºTj, min(t) > Pre(Pi, Tj),
mi(t + dt) = mi(t) - vj(t) × Pre(Pi,Tj) × dt

Pi Tjº,
mi(t + dt) = mi(t) + vj(t) × Post(Pi,Tj) × dt

where vj(t) is the instantaneous firing flow of Tj at time t.

When the concept of that an inhibitor arc of weight r from a place Pi to a transition Tj allows the firing of Tj only if the marking of Pi is less than r is used in a hybrid Petri net, we can extend the above-defined hybrid Petri net. If the inhibitor arc has its origin at a discrete place and has a weight r = 1, the corresponding transition can be fired only if mi > 1, actually, only if mi = 0 since mi is an integer. If the origin place is continuous, then a conventional value 0+ is introduced to represent a weight infinitely small but not nil. The new definition of an extended hybrid Petri net is similar to the definition of a hybrid Petri net (Def. 1), except that:

One can have, in addition, inhibitor arcs;

The weight of an arc (inhibitor or ordinary) whose origin is a continuous place takes its value in + {O+} instead of +;

The marking of a continuous place takes its value in + {O+} instead of +.

So far, the defined hybrid Petri net turns to be a flexible modeling process that makes sense to model biological processes, by allowing places using actual concentrations and transitions using functions.

In this paper, we exploit a hybrid Petri net tool named Visual Object Net++ (VON++) to model and simulate gene regulated networks. VON++ is a small, quick, uncomplicated and intuitive Petri net tool that supports both discrete event and timed event/conditin PN. Documentations to this tool can be downloaded via its web site at Figure 1 shows the basis elements of VON++, discrete place, continuous transition, continuous place and discrete transition connected with test arc, normal arc and inhibitor arc, respectively.

Figure 1: Elements of VON++.

The discrete transition is the active element in discrete event Petri nets. Transition can fire if all places connected with input arcs contain equal or more tokens than the input arcs specify. It can be assigned with a delay time. The continuous transition differs from the traditional the discrete transition; its activity is not comparable with the abrupt firing of discrete transition. The firing speed assigned to a continuous transition describes its firing behavior. It can be a constant or a function, i. e. transport of tokens according to v(t), in Figure 1, v(t) = 1.

The rates of bioprocesses are not defined within a Petri net, they should be specified separately. In automated control systems represented by Petri nets, execution of transitions usually depends on the presence of a specific number of tokens in all staring places. However, in most chemical and biological systems the rate of a process (transition) is defined by the mass action law. The change of tokens (or concentration) is proportional to the number of tokens (or concentration) in all starting places as expressed in the Figure 2. V is the rate of firing of the transition; k is a constant (called a rate coefficient in chemical kinetics); m3, m4 are the concentration of place S1 and place S2. Coefficient k varies with temperature, pressure, solvent, and other factors. As a result, v will become a function of several variables. Figure 3 indicates how token number (or concentration) of reaction intermediate P1 changes.

Figure 2: Presentation of transition rate in continuous systems.

Figure 3: Presentation of intermediate reaction rate.

Traditionally, kinetics has been taught in biochemistry courses in terms of the steady-state kinetics. This corresponds to a detailed study of the local properties of individual enzymes. However, one can go further and create kinetic models of the whole pathway. Such models are composed of coupled ordinary differential (for time courses) or algebraic (for steady states) equations. These equations are non-linear and most often without analytical solutions. This means that they can only be studied through numerical algorithm, such as the Newton method for solving non-linear equations and numerical integrators. After many years of development, now Petri nets have a mature mathematical algorithm and can solve NAEs and ODEs and stoichiometric matrices. But biochemical systems are also rich in time scales and thus require sophisticated methods for the numerical solution of the differential equations that describe them.

Mathematical methods and algorithms

Kinetic models of metabolic networks are becoming imperative not only the knowledge of more and more metabolic pathways has been acquired, but also the complexity of metabolic pathways requires an analytical and quantitative solution. Kinetic models can provide us spatiotemporal scale approaches and serve to check the consistency of metabolic theories with observed behaviors.

Related simulation environments

Many attempts have been made to simulate molecular processes in both cellular and viral systems. Several software packages for quantitative simulation of biochemical / metabolic pathways, based on the numerical integration of rate equations, have been developed. Table 1 shows a comparison of the most well-known metabolic simulation systems.

Table 1: A comparison of metabolic simulators with Petri nets approach
Tools Gepasia Jarnacb DBsolvec E-Celld
stoichiometry matrix presentation + + + +
Core algorithm and method MCA MCA MCA SRM, MCA
Pathway DB retrievable - - WIT/MPW, EMP KEGG, EcoCyc
Pathways graphic editor - + + -
Kinetic types +++ ++ ++ +
Virtual cell model - - - +
Simulation graphic display +++ ++ + +
Mathematical model accessible and modifiable + + ++ +
Data XML export * SBML* * *
User interface ++ + ++ +
Programming language C++ Delphi 5 C++ C++
a Gepasi (
b Jarnac (
c Dbsolve (
d E-Cell (
* SBML (Systems Biology Markup Language) ( is a description language for simulations in systems biology. It is oriented towards representing biochemical networks that are common in research on a number of topics, including cell signaling pathways, metabolic pathways, biochemical reactions, gene regulation, and many others. SBML is the product of close collaboration between the teams developing BioSpice (, Gepasi, DBSolve, E-Cell, Jarnac, StochSim ( and Virtual Cell (

Each tool possesses some prominent features while others little or no present. After a decade's development, Gepasi is widely used both for research and education purposes to simulate biochemical systems due to its powerful simulation engine and user-friendly interface. Jarnac, as a replacement of SCAMP, has a nice pathway graphic editor, called Jdesigner, enabling users to interactively draw a biochemical network. It has an SBW interface (System Biology Workbench), providing simulation capabilities for alternative clients. DBsolve is good at model analysis and optimization. By using numerical procedures for the integration of ODEs and/or NAEs to describe the dynamics of these models, DBsolve offers explicit solver, implicit solver and bifurcation analyzer. The primary focus of E-Cell is to develop a framework for constructing simulatable cell models based on the gene sets that are derived from completed genomes. Contrast to other computer models that are being developed to reproduce individual cellular processes in detail, E-CELL is designed to "paint a broad-brush picture of the cell as a whole". There is another program, named DynaFit (, which is also useful in the analysis of complex reaction mechanism.

In predicting cell behavior, the simulation of a single or a few interconnected pathways can be useful when the pathway being studied is relatively isolated from other biochemical processes. However, in reality, even the simplest and best studied pathway, such as glycolysis, can exhibit a complex behavior due to the high connectivity of metabolites. In fact, the more interconnections exist among different parts of a system, the harder it gets to predict how the system will react. Moreover, simulations of metabolic pathways alone cannot account for the longer time-scale effects of processes such as gene regulation, cell division cycle and signal transduction. When the system reaches a certain size, it will become unmanageable and non-understandable unless with decomposition of modules (hierarchical models) and/or presentation of graphs. In this sense, tools mentioned above appear faint. In comparison, Petri nets capture the basic aspects of concurrent systems of metabolism both conceptually and mathematically. The major advantages of Petri nets comprise a graphical modeling representation with sound mathematical background making it possible to analyze and validate the qualitative characteristics and quantitative behavior of a concurrent system, and to describe the system on different levels of abstraction (hierarchical models). In addition, the development of computer technology enables Petri net tools to have more friendly interfaces and possibility of standard data import/export supporting. We are motivated to exploit the Petri net methodology to model and simulate gene regulated metabolic networks.

The normal discrete system is easy to understand, so we emphasize here on the continuous one that proves useful to dynamic systems. We will describe some mathematical formulations that occur frequently in biology models. A general differential equation for a single state variable is dx/dt  =  flowin - flowout, while the expressions for the flowin and flowout can be quite complex as every bioprocess gives rise to its own system of differential equations involving many dependent variables (species concentrations) and many free parameters (reaction rate constants). The mass action law assumes that particles move incessantly. However, cellular metabolites are not like gas molecules. A metabolic reaction is very complex; interaction delay or saturation effect often exists in a metabolic system. In these cases, the mass action law becomes violated and should be replaced by equations that better describe the metabolic interaction while the rest of the algorithm remains the same.

PN model of metabolic reactions

In biochemistry, the most commonly used expression that relates the enzyme catalyzed formation rate of the product to the substrate concentration is the Michaelis-Menten equation, which is given as v = vmax · S / Km + S. An example of its Petri net model and simulation result is shown as below.

Figure 4: Petri net model of a simple enzyme catalyzed biochemical reaction (Michaelis-Menten reaction).

It is clear that such an enzymatic reaction is characterized by these two parameters: Vmax and Km, and biochemists are interested in determining these parameters from experiments. Fortunately, there are several biochemical reaction databases available for public such as BRENDA that provides enzymatic reaction kinetics. However, only for a subset of the well known pathways those kinetic parameters are complete, and moreover an enzymatic reaction can be affected by the presence of other compounds, i. e., the simplest form of the Michaelis-Menten equation does not account for the higher than first order substrate concentration dependence found in many allosteric enzymes. In the first case, we can introduce a general function v = Kapp · S to meet the lack of unknown parameters, where Kapp is the apparent rate constant. As we know the Michaelis-Menten equation is only valid when the concentrations of substrate and enzyme meet the precondition that [E] is not less than 0.001[S]. Considering the effect of enzyme concentration on the reaction rate in case the enzyme is sensitively regulated, i. e., the enzyme concentration is a variable of the model, the Michaelis-Menten equation can be written as , where kcat is known as turnover number. When there is more than one substrate involved in an enzymatic reaction, and its kinetic type is unknown, one gets processes more complicated than we discussed in the previous section. As the Michaelis-Menten equation is obviously invalid at this time, we simply apply the following function: .
For instance, given a two-substrate biochemical reaction, . Fortunately, if a two or more substrate biochemical reaction is already determined as one of the kinetic types such as competitive kinetics, Ping-Pong kinetics, etc. and available in the literature, the corresponding function is recommended to employ.

Model of Genetic regulatory networks

Although metabolic reactions determine anabolism and catabolism, the regulation of metabolism is mainly based on the regulation of gene expression because this determines whether a protein is present to carry out its particular metabolic reaction. Gene regulatory networks are the on-off switches and rheostats of a cell operating at the gene level. Based on interactions between genes and proteins and reactions of genes and proteins, they dynamically orchestrate the level of expression for each gene in the genome by controlling whether and how vigorously that gene will be transcribed into RNA. Each RNA transcript then functions as the template for synthesis of a specific protein by the process of translation. Process of gene regulatory networks is not restricted to the level of transcription, but also may be carried out at the levels of translation [Pyronnet et al., 1996], splicing [Yao et al., 1996], posttranslational protein degradation [Hochstrasser, 1996], active membrane transport [Weissmuller and Bisch, 1993], and other processes. In addition, such networks often include dynamic feedback loops that provide for further regulation of network architecture and output.

Building complete kinetic models of genetic regulatory systems requires detailed knowledge on reaction mechanisms, often, the following steps ad hoc are considered:

  1. The gene (DNA) is transcribed into RNA by the enzyme RNA polymerase.

  2. RNA transcripts are subjected to post-transcriptional modification and control: rRNA transcript cut into appropriate size classes and initial assembly in nuclear organizer; tRNA transcript folds into shape; mRNA transcripts are modified, noncoding sequences (introns) removed from interior of transcript; in eukaryotes, all RNA types must move to the cytoplasm via the nuclear membrane pores.

  3. Then mRNA molecules are translated by ribosomes (rRNA + ribosomal proteins) that match the 3-base codons of the mRNA to the 3-base anticodons of the appropriate tRNA molecules.

  4. Finally, newly synthesized proteins are often modified after translation (post-translation) before carrying out its function, which may be transporting oxygen, catalyzing reactions or responding to extracellular signals, or even directly or indirectly binding to DNA to perform transcriptional regulation and thus forming a closed feedback loop of gene regulation.

However, at present time, the information of the bioprocesses from genes to the gene-encoded products is still unclear or unavailable. So far, we can regard the unknown part as a black box of transition (one transition that can be visualized as the representation of a part of Petri net) and simplify the whole procedure as a higher level of abstraction (Fig. 5):

Figure 5: Petri net model simplification.

This modification does not involve changing the structure of the complete net and any modification to this subnet is reflected in the behavior of the original transition. Therefore, Petri net models are extensible and ways of plug-in concept, they can be extended without significant deviation from the existing structure.

As to model gene regulatory networks quantitatively, the state equations of the following form are used to model bioprocesses such as activation of proteins, binding of proteins to genes, binding of RNA polymerase and so on.

If state[i](condition), then  = state[i](consequence)

For example, the concentration of the gene product is state[i]. The condition contains regulatory terms for this gene and describes whether the gene is being expressed or not. It depends on the state of the cell, and may contain models for promoters, enhancers, other proteins, nucleic acid, etc. The consequence then describes the result of condition changing, here, the rate of gene expression.

So the differential mass balances describing the concentration of mRNA and of the encoded protein can be given as:

If ( (Gene, transcriptional factor(s), RNA nucleotides, binding of RNA polymerase, etc.) not (Repressors, etc.))
Then (transcription is initiated and mRNA is produced,
If (Modified mRNA, tRNA, initiation factor(s), amino acid, binding of ribosome, etc.)
Then (the gene-encoded protein is synthesized,

Where kts and ktl are the rates of transcription and translation respectively, kd is the rate of degradation and kr is the rate of consumption of biochemical reaction. GPC denotes the concentration of the binding complex of gene, TFs, RNA ploymerase etc. DNA is a stable molecule, but mRNA and proteins are constantly being degraded by cellular machinery and then to be recycled. Specifically, mRNA is degraded by a ribonuclease (RNase), which competes with ribosomes to bind to mRNA. If a ribosome binds, the mRNA will be translated, if the RNase binds, the mRNA will be degraded. Proteins are degraded by cellular machinery including proteasomes signaled by ubiquitin tagging. Protein degradation is regulated by a variety of more specific enzymes (which may differ from one protein target to another). In practice, the first-order rate constant of degradation kd often is replaced by a half life H, and the degradation rate is expressed as , where H = 0.693/kd. mRNAs have specific half-lives ranging from hours to days.

Regarding the model of binding procedures which also are common phenomena in signal transduction, say,

a general model , where Kb is the binding constant, is presented for systems consisting of one subject Ai binding with other subjects.

As in many situations, the information of gene regulatory pathway and mechanism is not available and one needs to take recourse to more approximate models. In this sense, the discrete model will be favorable.

Diffusion transportation

Most of the models deal with the amount of metabolites in a cell. In the simplest case, we may be able to assume that the cell is a "well-mixed pool", i. e., the amount of metabolites, enzymes, etc. is uniform across the cell. In many situations, however, concentration gradients exist which will affect the local rate of biochemical reactions, in particular for large systems and different compartments, we must consider explicitly the effect of diffusion or transportation.

In general, if concentration gradients exist within the spatial scale of interest it is highly likely that diffusion will have an impact on the modeling results, unless the gradients change so slowly that they can be considered stationary compared to the timescale of interest. A growing number of modeling studies [Markram et al., 1998; Naraghi and Neher, 1997] have emphasized the important effects of diffusion on molecular interactions. Moreover, many bioprocesses take place in different compartments in a cell. For instance, glycolysis conducts in cyotoplasma while TCA in mitochondria. Membranes play an important role to separate these bioprocesses and meanwhile maintain the normal transportation of metabolites inside and outside of them. In addition, signal transduction also occurs across the membranes.

So far, in order to model a metabolic network, not only all effect of metabolites and reaction behaviors but different compartments should be considered. Diffusion, facilitated diffusion and active transport could be the very important physical effects in the models. We will focus on the membrane transportation. The rate of penetration of a metabolite across a membrane is related to the concentration gradient by Fick's Law of Diffusion:

Rate of penetration

where [S]out and [S]in are concentrations of metabolite outside and inside the membrane, respectively; D denotes the diffusion coefficient (D decreases with the size of the metabolite); A is the area of membrane (the greater it is, the more metabolite that can pass); is the partition coefficient ( increases with increasing solubility) and dx is the membrane thickness (the greater the thickness, the slower the rate). Usually, is called the permeability constant, a constant for a given substance moving through a given membrane.

In carried systems, the carrier exhibit saturation kinetics, so that "Michaelis-Menten equation" formula might be used to describe such a process where low Km means a high rate of "affinity and transport", and high Km a low "affinity and transport" rate. Some metabolites and/or signals (hormones) may modify carriers and change Km. Vmax is related to "carrier mobility", the total number of carriers present.

Petri net modeling algorithm

Modeling algorithm and analysis of hybrid Petri nets can be done by the following procedures:

Draft network construction
Normally, a Petri net model is built manually by drawing places, transitions and arcs with mouse events. Fortunately, The XML based Petri net interchange format standardization which consists of a Petri net markup language (PNML) and a set of document type definitions (DTD) or XSL Schema is coming into being and intended to be applied. Several Petri net tools such as PNK, Renew, Design/CPN and GON are currently being equipped with an XML based file format exchange. We have developed an environment to extract data of metabolic networks from KEGG, BRENDA and RegulonDB and transform them into XML based files that can be used by PNK and Renew to display the Petri net models automatically [Chen, 2002; Chen et al., 2002b].
Data preparation
The main feature of metabolic processes is that the concentration of metabolites will influence the reaction activity of bioprocesses. Therefore, the actual concentration of any metabolite is an important component of the model. Although some data nowadays are available to public via the Internet and other sources, some other data may be not complete. It requires a time consumption searching throughout the literature. In case of the unavailability of concentration values, additional experiments might be required. However, some computational prediction and experiential data also might help. Assignment of initial value of places is made after data gathering.
Determination of the kinetics
A series of predefined types of kinetics which are frequently used in biochemical reaction models are collected by SBML [Hucka et al., 2001]. However there are some circumstances in which the kinetic types are not yet defined. Then the kinetic modeling strategy is to be applied to ouline the kinetic properties of each reaction. All relationships and influences of metabolites are to be fixed by introducing the conresponding variables into the self-defined functions.
Stoichiometry check
Many metabolic pathway schemes contain mass conservation relations that must be taken into account in order to carry out the simulation. To check the mass conservation relations of a model we can go to the original reaction data from databases or the literature, or calculate them manually. In fact, we construct the model with the identification of reaction stoichiometry. Otherwise, it will lose something when the simulation is carried out because in continuous Petri nets, the weight of arcs is disable, so that all components involved in the reaction are changed with same rate which is defined by the transitions function. However, in the reaction 2*A+B=>C, the change of A should be twice than that of reaction rate. In VON++, unfortunately, we have to add more transition from A with the same function in order to obey the mass conservation law.
Parameter tuning and simulation
To build a model precisely requires as more as possible the variables and parameters involved in the metabolic network. The values of variables and parameters are determined either by experimental methods or deduced from other related values. However, it is impossible or sometimes unnecessary to put all variables and parameters into a model. The model is plausible when main influences are included. On the other hand, because of different purposes and situations, most data from laboratory do not fit the model very well and vice versa. We have to compare and tune the differences in order to find suitable ones. Then the effects of various parameters on the gene regulated metabolic networks and their relations can be understood. The key enzymes/proteins and intermediates related in the metabolic pathway can be determined, which will provide the necessary information to solve the metabolic bottlenecks.

Case study

Inborn errors of metabolism

Inborn errors of metabolism are characterized by a block in a metabolic pathway, a deficiency of a transport protein or a defect in a storage mechanism caused by a gene defect. The defect gene leads to an absent or wrong production of essential proteins, especially enzymes that enable, disable or catalyze the biochemical reactions of metabolic networks. Thus, these disorders of the metabolism result in a threatening deficiency or accumulation of intermediate metabolites in the human organism and their following corresponding symptoms.

For inborn errors of metabolism a lot of data is available in different databases accessible via the Internet. Most inborn errors of metabolism are included in OMIM (Online Mendelian Inheritance in Man Aside from the major molecular biological databases, e. g., GenBank, EMBL, TRANSFAC, KEGG and BRENDA, our group's Metabolic Diseases Data Base (MD-Cave) has been developed to simplify the collect and persistent storage of knowledge about inborn metabolic diseases [Hofestaedt et al., 2000b; Kauert et al., 2001]. Although the amount of this electronically available knowledge of genes, enzymes, metabolic pathways and metabolic diseases increases rapidly, they give only highly specialized views of the biological systems. It is clearly that the next task is to integrate all this knowledge and make it biotechnologically and medically applicable. The MD-Cave is developed as such a bioinformatics system for representing, modeling and simulating genetic effects on gene regulation and metabolic processes in human cells. In the following section, we will present a case study that emphasizes on the modelling and simulation of the gene regulated urea cycle network by using the hybrid Petri net.

Urea Cycle and its regulation

In human cells, excess nitrogen is removed either by excretion of NH4+ (of which only a little happens) or by excretion of urea. Urea is largely produced in the liver by the urea cycle, a series of biochemical reactions that are distributed between the mitochondrial matrix and the cytosol (Figure 6). The cycle centers around the formation of carbamyl phosphate in hepatocyte mitochondria to pick up NH4+ and incorporate it into ornithine to make citrulline. Citrulline is then transported to the cytosol where aspartate is added. As urea is removed it is converted back to ornithine that goes back into the mitochondria to start over.

Figure 6: Key enzymes in regulation of urea cycle in cells. CPS1: Carbamyl phosphate synthetase, EC; OTC: Ornithine transcarbamylase, EC; ASS: Argininosuccinate synthetase, EC; ASL: Argininosuccinate lyase, EC; ARG: Arginase, EC

Deficiencies in the urea cycle enzymes lead to excessive NH4+ and accumulation of its intermediates, resulting in neurological disorders. Any of five enzymes of the urea cycle may be deficient and lead to carbamyl phosphate synthetase (CPS) deficiency, ornithine transcarbamylase (OTC) deficiency, citrullinemia, argininosuccinic aciduria and argininemia.

Although the urea cycle was discovered by Hans A. Krebs early in 1930's, the modeling and simulation of the urea cycle so far have never been systematically explored. This case study therefore will attempt to model and simulate it. A model will show the interrelations of main metabolites invovled in the urea cycle. The simulation will be used to test the physico-chemical limitations and feasibility of certain proposed reactions as well as disease occurrences.

Figure 7 shows a graphical representation of the urea cycle metabolic network using the objects presented above for describing entities and interactions. It shows an intricate network that links entities and interactions. This network includes not only the succession of biochemical reactions that lead to the transformation of CO2 and NH4+ to urea, but also the regulation of gene expression and enzymatic activities. It furthermore displays (e. g. asparate, fumarate) the links to other pathways, which are not detailed on the graph to preserve clarity.

Figure 7: Schematic diagram of urea cycle metabolic network.
Data sources: Metabolic pathway (enzyme reactions) from KEGG and BRENDA; Gene regulatory: TRANSFAC; Drug information: Metabolic Diseases Drug Database (MDDrugDB) (; drawing by Ralf Kauert).

PN model

Based on the proposed modeling strategy, a hybrid Petri net model of the gene regulated urea cycle metabolic network is presented (Figure 8). The model of the intracellular urea cycle is made of the composition of gene regulatory networks and metabolic pathways. It comprises 155 Petri net elements, 14 kinetic blocks, 39 dynamic variables, and 22 reaction constants. Experimental data, partially listed in Table 2, are used for the initial evaluation of certain parameters of enzymatic reactions within the system. The values of the model parameters lacking in the literature are verified through numerical experiments or modifed from several references.

Table 2: Some kinetic parameters of enzyme reactions in human cells
Enzyme Substrate (mean concentration, mM) Compartment Km (mM) Kcat · E0 Reference
CPS1 HCO3- (0.05), NH4+ (0.025) mitochondia HCO3-, 6.7
NH4+, 0.8
MgATP, 1.1
1.5 [Pierson and Brien, 1980]
OTC carbamyl phosphate (0.001), L-ornithine (0.05) mitochondia CP, 0.16
L-ornithine, 0.40
2 [Charles et al., 1997]
ASS L-citrulline (0.02), L-aspartate (0.325) cytoplasm L-citrulline, 0.03
L-aspartate, 0.03
0.15 [Charles et al., 1997]
ASL argininosuccinate (0.034) cytoplasm argininosuccinate, 0.017 0.2 [Charles et al., 1997; Pierson and Brien, 1980]
ARG L-arginine (0.06) cytoplasm L-arginine, 10 1.7 [Charles et al., 1997]

The dynamic behavior of the model system, such as the metabolite fluxes, NH4+ input and urea output are well described with continuous elements; while control of gene expression are outlined with discrete ones due to the insufficiency of explicit expression data. Nevertheless, when explicate knowledge about expression levels of the enzymes are available; it is possible to exploit our model of gene regulatory network to handle realistic gene expression data with the state equations. The initial values of variables were assigned and tuned so that the model system behavior would comply maximally with available experimetnal data on the dynamic characteristics of the system's behavior, based on the following considerations:

The availability of ammonia or amino acides (denoted as NH3) is ingested continuously from plasma into mitochondria with a stable speed, i. e. the changes of ammonia concentration due to the rate of protein metabolism are not taken into account. The concentration of nitrogen excreted (urea) in plasma ranges from 3 mmol/L to 8 mmol/L and thus is regarded to be discharged with a certain rate. The degradation rates of an enzyme is 0.001 times of its concentration. The places of the main metabolites are directly linked to the transitions. Reaction rates assigned to these transitions are interpreted with differentail equations. However, from reality point of view, these transitions involve more than one variable that are presented in differential functions. In order to get a better understanding of these relationship, serval test arcs are used, e. g. the test arc between asparate and transition of ASS. There are no real input and output within these arcs, but the places linked are exploited by the transition firing speed.

In the model, inhibitor arcs are also used to present negative effects of repressors and/or inhibitors to gene expression. On the biochemical reaction level, negative effects of metabolites are expressed as enzyme inhibitions that include competitive inhibition, noncompetitive inhibition, irreversible inhibition and feedback inhibition. Sequentially, the regulation of urea cycle enzyme activities can be brought in these two ways: First, the gene expression that is regulated by activators and inhibitors controls the enzyme synthesis, while enzyme synthesis and degradation determine the amount of the enzyme. Second, the activity of the enzyme can be altered during the metabolic catalysis.

Figure 8: Petri net model of the gene regulated urea cycle and its dynamic simulation result.

The formalization of the urea cycle model allows the quantitative simulation of metabolic pathways. Dynamics of the main components on the model regulating the urea cycle are shown in Figure 8, too. Moreover, several tests on interfering the fluxes intentionally are conducted and results are observed (Table 3).

Table 3: Interfering tests on the urea cycle Petri net model
Interfering test Value of metabolites Urea cycle defect
NH4+ Citrulline (plasma) Argininosuccinate (plasma) Ornithine(plasma) Arginine
CPS1 blockade Carbamylphosphate synthase deficiency
OTC blockade Ornithine transcarbamylase deficiency
ASS blockade Argininosuccinate synthase deficiency
ASL blockade Argininosuccinase deficiency
ARG blockade Arginase deficiency
Membrane transportation blockade _ HHH syndrome

As the urea cycle operates only to eliminate excess nitrogen. High concentration level of ammonia in the cell results in hyperammonemia which is a typical fatal event, coma and death ever been reported. Laboratory studies can reveal elevated arginine levels, mild hyperammonemia, and a mild increase in urine orotic acid. The diagnosis now can be confirmed by enzymatic analysis in the model. On high-protein diets or under starvation state, proteins are degraded and amino acid carbon skeletons are used to provide energy, thus increasing the quantity of nitrogen. But the amino nitrogen must be excreted. To facilitate this process, enzymes of the urea cycle are controlled at the gene level to enhance the concentrations of enzymes. As the urea cycle takes place both in mitochondria and cytoplasma, the effects involved also come from the membrane transportation. Some mitochondrial membrane diseases, e. g. ornithine transporter deficiency, surely effect the transportation of ornithine into matrix and results in high concentration of ornithine accumulation in plasma, which gets a feedback to the transition of arginine into urea and finally hyperammonnemia. From the model we know the treatment for defects of the enzymes in the urea cycle could be either limiting the input of ammonia (limiting protein intake) or replacing the missing intermediates from the cycle (supplementing with arginine or citrulline). Patients with OTC deficiency benefit from citrulline supplementation because citrulline can accept ammonia to form arginine.


The case study shows that the Petri net allows easy incorporation of qualitative insights into a pure mathematical model and adaptive identification and optimization of key parameters to fit system behaviors observed in gene regulated metabolic networks. The advantages of applying hybrid Petri nets (HPN) to model and simulate are: (i) The HPN model has a user-friendly graphical interface that allows an easy design, simulation and visualization. (ii) With the discrete and continuous events, HPN can easily handle gene regulatory and metabolic reactions. (iii) The inhibitor arc is useful for mechanistic studies to learn how enzymes interact with their substrates, to know the role of inhibitors in enzyme regulation and gene expression. Moreover, powered with mathematical equations, simulation is executable and dynamic results are visible.

As in the cell, there are usually hundreds of interconnected metabolic pathways and gene regulatory networks and control of these presents more complex features. It is feasible to extend the Petri net model in a plug-in way. A large complex network model can be handled with the same set of structural and behavioral properties. When HPN models are applied to such a large network, the hierarchical concept makes it possible to develop a generalized variant of HPN at a global level. On the other hand, the subnet of the Petri net model provides us the basic model of which we already know its inner behavior and functions. Then we can construct a system by plugging together sub-models in order to understand the higher-level system and to predict its behavior.

With rapid development of PN, many important extensions to the above general Petri nets classification appear. In the past few years, a number of Petri net tools have been used to model and simulate metabolic pathways and gene regulatory networks, e. g. UltraSAN, Design/CPN and VON++. More tools can be found at However, different tools have their characteristics and cannot embed various functions. Almost all Petri net tools were intended to model manufacture, distribution, and communication systems rather than biological systems. It requires the collaboration of biologists and Petri net researchers to construct a specific Petri net tool that contains all necessary features. Fortunately, Matsuno et al. [2001] adapted the VON++ to the GON program, and have made a significant progress.

Provided with rich information about biochemical reactions and gene regulation, availability of various biological databases, building an integrative model of the whole cell (virtual cell modelling) that incorporates gene regulation, metabolic reactions and signal transductions is becoming a promising field in the post-genomic era. Several projects have been established. The challenge created with Petri nets is to understand how all the cellular proteins work collectively as a living system. Using powerful Petri nets and computer techniques, data of metabolic pathways, gene regulatory networks and signalling pathways could be converted for Petri net applications. Thus, a Petri net based virtual cell model could be implemented, and the attempt to understand the logic of the cell could be accomplished.


The work was partly supported by the Ministry of Science and Art of the Government of Sachsen-Anhalt, and by the German Research Foundation (DFG) graduate program "Bioinformatics" in the Uni-Bielefeld.