Topological analysis of metabolic networks based on Petri net theory

---

Abstract

Petri net concepts provide additional tools for the modelling of metabolic networks. Here, the similarities between the counterparts in traditional biochemical modelling and Petri net theory are discussed. For example the stoichiometry matrix of a metabolic network corresponds to the incidence matrix of the Petri net. The flux modes and conservation relations have the T-invariants, respectively, P-invariants as counterparts. We reveal the biological meaning of some notions specific to the Petri net framework (traps, siphons, deadlocks, liveness). We focus on the topological analysis rather than on the analysis of the dynamic behaviour. The treatment of external metabolites is discussed. Some simple theoretical examples are presented for illustration. Also the Petri nets corresponding to some biochemical networks are built to support our results. For example, the role of triose phosphate isomerase (TPI) in Trypanosoma brucei metabolism is evaluated by detecting siphons and traps. All Petri net properties treated in this contribution are exemplified on a system extracted from nucleotide metabolism.

Key words: Petri nets, elementary flux mode, metabolic networks, P-invariant, T-invariant, incidence matrix

---

Introduction

The aim of theoretical approaches is to build models, which should facilitate the study of real systems. A huge variety of models has been developed in theoretical biology. For example, biochemical reaction networks are the subject of extensive modelling studies [Fell and Wagner, 2000; Heinrich and Schuster, 1996; Jeong et al., 2000; Peleg et al., 2002; Schuster et al., 2002a, 2002b; Teusink et al., 1998]. This modelling has important applications in biotechnology [Klamt and Stelling, 2003; Liao et al., 1996; Schuster et al., 2000; Van Dien and Lidstrom, 2002] and genome research [Dandekar et al., 1999; Price et al., 2002; Förster et al., 2002]. Specific features of biochemical systems are that most reactions are catalyzed by enzymes and that many reactions utilize more than one substrate (reactant) and/or generate more than one product. Reaction systems that only involve isomerizations (that is, reactions with one substrate and one product) can be depicted as graphs and their properties can be studied from the point of view of graph theory. However, real biochemical networks cannot be represented as graphs due to bi- or multimolecular reactions. These cases would require that arcs linking three or more nodes exist. One, quite complicated, approach to coping with this situation is to consider the groups of substances on each side of reaction arrows as nodes (so-called complexes) [Clarke, 1980; Horn and Jackson, 1972].

A simpler solution is offered by Petri net theory [Reisig, 1985; Starke, 1990]. Two kinds of nodes are considered: places and transitions. The nodes and arcs between them represent the static structure, while some more elements/components such as tokens indicating time-dependent weights of places are used to describe the dynamics. Beside place/transition nets, also condition/event Petri nets have been proposed in the literature. Here, we only deal with the first type. Petri nets can be employed for the graphical description of processes. They allow us to understand more intuitively the temporal evolution of systems by considering flows of tokens through the nets. They offer also an appropriate formalism for the analysis of biochemical networks, as has been pointed out earlier by several authors [Genrich et al., 2001; Heiner et al., 2000; Heiner et al., 2001; Hofestädt, 1994; Küffner et al., 2000; Oliveira et al., 2003; Peleg et al., 2002; Reddy et al., 1996], while the above-mentioned modelling approaches [Fell and Wagner, 2000; Heinrich and Schuster, 1996; Jeong et al., 2000; Leiser and Blum, 1987; Schuster et al., 2002a, 2002b; Teusink, 1998] are independent of Petri net theory. In the present paper, we shall show the correspondence between concepts in both, Petri nets theory and traditional metabolic network analysis. Some examples will help the reader see the similarities and to exploit them.

In the present paper, we shall focus on the use of Petri nets for the topological analysis of biochemical networks rather than for the analysis of the dynamic behaviour. In particular, we shall deal with various invariants and other features in these nets such as boundedness and liveness and reveal their biochemical meaning. Moreover, we shall discuss the appropriate treatment of source and sink metabolites.

---

Similarities between Petri net theory and traditional biochemical modelling

In graphical representations of Petri nets, circles are used for places, while rectangles stand for transitions (Fig. 1). The correspondence place - substance (in biochemistry often called metabolite) and transition - reaction/enzyme is obvious. Metabolic networks have a static level - the stoichiometry, and a dynamic one, characterized by fluxes. The stoichiometric coefficients indicate how many molecules of a substance have to react to produce how many molecules of product. The stoichiometric coefficients are described by arc weights. Thus, the stoichiometry matrix containing these coefficients corresponds to the incidence matrix of a Petri net (see below).

A further object - the token - was introduced to describe the dynamics of a Petri net. It is denoted by a solid dot () inside the circles representing places. In ordinary Petri nets, the tokens are indistinguishable. They indicate the presence or absence of a condition, a signal, or a resource. In our case, the number of tokens in a place stands for the number of molecules of that metabolite existing at a given moment. Alternatively, tokens may correspond to any predefined unit measuring the amount of substance, such as mole, millimole etc. However, this brings about that non-integer token numbers should be admitted. This leads to continuous Petri nets, which are currently being developed [Alla and David, 1998; Matsuno et al., 2000].

The tokens that exist in the system at a given time describe the state of the system. This is called marking, M(P) [Reisig, 1985; Starke, 1990]. The system state changes when a transition fires. This can happen only if the transition is active/enabled, that means that every place from the input places set (Fig. 2 and Table 1) of the considered transition has at least as many tokens as the weight of the corresponding arc. The set of input and output places of a transition t is denoted by  t and t , respectively. This set corresponds to the metabolites that act as reactants and products, respectively, in the reaction t. The new state is obtained by subtracting from each input place of the considered transition a number of tokens equal to the weight of the corresponding arc and adding in each output place a number of tokens equal to the weight of the corresponding arc (Fig. 2 and Table 1). Formally, we can also speak about the input and output transitions set of a place p - p and p, respectively, which contain all the transitions which produce, respectively consume, the metabolite p.

It is useful to know between which places (P) and transitions (T) there exist arcs. For this purpose, two mappings describing weights were introduced (see also Table 1): pre: PT N and post: TP N (with N denoting the set of natural numbers). One can think about them also as matrices. The rows in pre correspond to places and the columns to transitions, while in the matrix post, the roles of rows and columns are transposed. The entries of these matrices have a nonzero value (equal to the weight of the arc), if an arc exists, and zero otherwise.

Further, the topological structure of a Petri net can be represented by an integer matrix, C, called an incidence or flow matrix. C is an nm matrix whose m columns correspond to the transitions and n rows correspond to the places of the net. The following relation holds true: C = post^T - pre. The mappings pre and post can be reconstructed from the matrix C in the following simple way: post(t_j, p_i) = max{C_ij, 0}, pre(p_i, t_j) = min{C_ij, 0}.

It is worth finding whether another state can be reached from a given state. This is related with the property of reachability. In metabolic networks, we can search all possible subsequent states, knowing the initial state of resources. Another interesting problem is to deduce all appropriate initial states from a desired later state.

Places can hold an arbitrary number of tokens or can be restricted by a given number - capacitated places. In Fig. 3, the unnamed places are considered external (with inexhaustible numbers of tokens). Transitions T₁ and T₂ are activated if there is at least one token in the place ATP (the currency of metabolic energy), respectively ADP. If the initial marking is [c, c, c, c, 1, 0], transition T₁ can fire and produce the marking [c, c, c, c, 0, 1]. Now, transition T₂ is enabled. It fires and we obtain again the initial marking. In this example, ATP + ADP = 1, independent of the system state. This is a conservation relation, which leads to boundedness of the capacity of all the internal places. Usually the places of biological systems are not considered to be limited because the limitation due to the finite size of living cells is not critical to most biochemical processes. There are only cases where the limitation comes from a conservation relation such as in the above case.

Another situation important in biology is the presence of inhibitors. The corresponding Petri net model can be extended by a special element, called inhibitory arc (Fig. 1i). The inhibitor is represented by a place. If there is a token at that place, the transition is not enabled, so it does not fire.

Note that the incidence matrix corresponds to the stoichiometric matrix [Érdi and Tóth, 1989; Heinrich and Schuster, 1996] for metabolic networks if the Petri nets are pure. That means that the networks do not involve self-loops (Fig. 4), because self-loops cannot be represented in the incidence matrix: a coefficient of 1 and a coefficient of -1 cancel each other to yield zero in the matrix, thus losing track of the self-loop. Thus, we should identify the situations that produce self-loops and the way to treat them without losing the biological meaning. One such situation occurs when enzymes are considered as normal substrates. There are algorithms for deleting/eliminating self-loops: First, one can treat the enzyme as a parameter (as usually done in biochemical modelling, see Heinrich and Schuster, 1996) rather than explicitly as a substrate. Second, for each loop, one may introduce another place and another transition. If all weights equal unity, the number of tokens in the new place is the difference between the capacity of the old place and the number of tokens that this old place contains already (Fig. 4b). By this construction, the new place corresponds to the enzyme-substrate complex, and the so-called overall reaction (old transition) has been decomposed into half-reactions (the new transitions). This construction can be applied also when the arcs have a multiplicity larger than 1, but care should be taken that the new arcs might have also such multiplicities. If p' and t' stand for the new place, respectively for the new transition, the new arcs have to respect the formulae: pre(p',t) = post(t',p') = pre(p,t') := pre(p,t), while the old arcs keep the same multiplicity. This situation can be nicely illustrated by the biochemical example of autocatalysis: A + B give 2B (Fig. 5). This reaction cannot simply be reduced to A gives B, because a small quantity of the product B is needed to start the reaction.

Second, a reversible reaction might be wrongly interpreted as a self-loop (Fig. 6). So, for each reversible reaction, we should consider only one flux direction and for the opposite one, we should introduce another transition. So, metabolic networks can be transformed into pure ones. Alternatively, one might think of defining reversible transitions. This has, however, not been dealt with so far in the Petri net literature.

Generally, Petri nets can be designed in different ways, called top-down and bottom-up. The first method supposes to start with a very generalized form of the system and then, to detail it as much as possible, until the basic units are reached. The second one starts with the "atoms", building modules which are then joint to model the real system. Sometimes it is advantageous to combine both approaches. The importance of modularity is expressed by the ancient saying "divide et impera".

---

Modelling of external metabolites

In metabolic networks, one needs to differentiate between internal and external metabolites. The former are totally produced and then consumed in the given network, while the external metabolites represent sources or sinks [Heinrich and Schuster, 1996]. Their amount is usually assumed to be constant, due to availability in large excess or well-tuned biological regulation. If one considers the given net as a part of a larger system, the external metabolites are a kind of boundary; or connection points with the remaining part, in which pathways producing or consuming these metabolites exist. An extension of the system so as to include those pathways is not useful, as the following example illustrates. Glycolysis (the well known pathway of sugar degradation [Stryer, 1995]) contains a sequence of reactions that transforms glucose into pyruvate, producing ATP. Glucose, pyruvate, ATP and also some other metabolites are usually considered "external" for this pathway. We might include, in the model, a reaction or pathway that consumes pyruvate, for example, for producing alanine. However, alanine then would be an external metabolite. The model needs to be delimited somewhere. Algebraically, the external metabolites can usually be identified also in the incidence matrix. Provided that each internal metabolite is both produced and consumed within the net, the external metabolites correspond to those rows in which all the coefficients have the same sign.

The modelling of external metabolites can be done in different ways. One of them is to fill all initial places with an inexhaustible number of tokens (modelled by infinity). For the sink places, one could allow them to accumulate tokens but has to take care in computing T-invariants (see below). If it is preferred to use finite token numbers for the initial places, one could redefine the firing rule for the transitions that have initial places in their preset (see Table 1) or final places in their postset, in such a way to not "consume" the input-place tokens and to not "produce" final-place tokens,

where I is the set of initial places and F is the set of final/terminal places. Note that the firing rule for the internal metabolites is

Another possibility is to connect sink places with source places by additional transitions so that a circular flow occurs [Heiner et al., 2000]. However, it is difficult to find which places have exactly to be connected with each other, because such transitions could impose unrealistic constraints on the flow ratios. For example, one cannot regenerate carbon atoms from outgoing nitrogen atoms. A solution can be to use coloured Petri nets, in which different atom groups can be modelled by tokens of different colour. Starke [1990] proposed, as another way of description, not to include the initial and final places in the net. Thus, the boundary is made up of transitions without presets or without postsets. The initial transitions do not need any tokens to fire. In the traditional modelling of metabolic networks, a similar description is indeed sometimes used for external metabolites that are of minor importance, such as inorganic phosphate, water, protons etc. However, applying this technique to all external metabolites has the drawback that they are not made explicit so that overall molar yields cannot be computed.

Here, we propose an alternative method. For each initial place, we add an arc feeding from the transition back to this place (Fig. 7) and use the firing rule (2) both for internal and external metabolites. This guarantees that the number of tokens in the initial places remains unaltered. For each final place, we add an arc feeding from this place back to the transition producing it. To guarantee that the transition can always fire, at least one token should be put in the final place at the beginning. However, one should be aware that this generates self-loops, so that the Petri net is no longer pure. Thus, as far as the external metabolites are concerned, the incidence matrix does not equal the stoichiometry matrix. This is no problem since the external metabolites are not usually included in the stoichiometry matrix [Heinrich and Schuster, 1996].

---

Invariants in Petri nets

When studying a system, it is always appropriate to begin with the study of its structural invariants. They help in analysing the system's behaviour and checking its logical properties. The same is true for Petri nets describing biochemical networks because the structural invariants do not depend on kinetic enzyme parameters, which vary due to external influences and internal fluctuations. Basically, there are two types of invariants in Petri nets: P-invariants and T-invariants [Reisig, 1985; Starke, 1990]. P-invariants (place invariants) are vectors, Y, with the property that multiplication of these vectors with any place marking reachable from a given initial marking yields the same result. If M₀ is the initial marking and m is some arbitrary marking, the relation Y^T·M = Y^T·M₀ describes a P-invariant and is called relation of marking conservation. Taking into account consecutive markings (that are obtained by firing of only one transition), it results that Y^T·col_t(C) = 0, for each transition t, where C is the incidence matrix. That means that, algebraically, these vectors are solutions of the equation

Invariants in Petri nets correspond to basic concepts in traditional biochemical modelling. In particular, P-invariants express conservation relations for metabolites, as becomes clear in the scheme shown in Fig. 3. This net has the P-invariant ATP + ADP = const. In general, eq. (3) is known, as for metabolic systems, as the general form of conservation relations [Horn and Jackson, 1972; Clarke, 1980, Heinrich and Schuster, 1996]. In most cases, these relations express the conservation of atom groups [Schuster and Hilgetag, 1995; Schuster and Höfer, 1991]. In the example in Fig. 3, the adenosine moiety is conserved.

In algebraic terms, invariants form a linear vector space. This implies that if I₁ and I₂ are invariants, also c₁I₁ + c₂I₂ with c₁, c₂ being real numbers, are invariants of the net [Reisig, 1985]. For example, if a biochemical net involves the P-invariants ATP + ADP = const. and NAD + NADH = const. [Stryer, 1995], then also ATP + ADP + 2 NAD + 2 NADH = const. is a P-invariant. Normally, one chooses invariants with the smallest integer coefficients and tries to decompose the invariants into the minimal terms (such as ATP + ADP = const.). This leads to the concept of minimal invariants (see below).

In order that conservation relations reflect the conservation of atom groups, the coefficients in these relations have to be non-negative. This leads to non-negative conservation relations [Schuster and Hilgetag, 1995; Schuster and Höfer, 1991]. They correspond to semi-positive P-invariants in Petri nets [Colom and Silva, 1990]. If all substances are involved in such conservation relations, the system is called conservative [Érdi and Tóth, 1989; Horn and Jackson, 1972]. This implies that a positive linear combination of all concentrations (token numbers in Petri nets) is constant in time,

where denotes time. can be, for example, the number of some sort of atoms. In a closed system, that is, a system without external metabolites, there is always one relation of the type (4) in which represents total mass. In addition, there may be further relations of the type (4). If there is a positive linear combination that increases in time, the system is called superconservative [Érdi and Tóth, 1989]:

If the sum in eq. (5) decreases in time until it reaches zero, the system is subconservative. The terms conservative, superconservative and subconservative have also been coined for Petri nets. The program INA developed by Starke and coworkers [www.informatik.hu-berlin.de/~starke/ina.html] determines whether a Petri net has one of these properties. Note that these three cases do not cover all networks. In fact, biochemical networks are usually open systems with a throughput of mass, as described by a flux between external metabolites. Therefore they may, depending on conditions, have a positive or negative mass balance, so that they usually belong to none of these classes.

P-invariants are useful in checking the property of mutual exclusion. Two transitions are called to be in mutual exclusion if there is no reachable marking that allows the two transitions to fire simultaneously. The first step is to identify the marking set that could characterise the simultaneous activation of the specified transitions. This marking is reachable only if it satisfies the conservation relation given by the P-invariants. For example, if four places need to have at least one token each to enable two transitions to fire simultaneously, while the conservation sum is three, mutual exclusion occurs. Such a case is irrelevant for metabolic networks provided that the molecule numbers are large enough. Alternatively, if the token numbers represent millimoles or the like, token numbers need not be integer, so that mutual exclusion is no problem either.

Often, two transitions are in mutual exclusion when they compete for the same input places set. If the tokens are indivisible and once a transition takes the existent tokens, the other transition cannot fire. If the tokens needed to reactivate the competing transitions are simultaneously regenerated for each conflict case, the net is called persistent. In this case, the two transitions do not deactivate each other. At first sight, some metabolic networks seam to be non-persistent because different enzyme reactions often compete for the same substance. However, in real metabolic nets, even if the quantity of the common resource is very small, the concurrent reactions share it, maybe in different percentage according to the various reactions rates. Until now, there are no techniques based on Petri nets that can model accurately this behaviour.

Biochemical networks often reach, after some initial transient, a stationary state. More concretely, this is the case when the kinetic properties of the networks are such that the stationary state is asymptotically stable, as is often the case [Clarke, 1980; Heinrich and Schuster, 1996]. At steady state, the following equation holds:

where V stands for the vector of net fluxes. They correspond to the flow of tokens per time in Petri nets. Special attention has to be paid to the involvement of external places in matrix C. If we connect outputs with inputs by additional transitions as explained in the previous section or additional arcs are added to create self-loops next to initial and final places, C can contain the coefficients both for internal and external places. Otherwise, it should only contain the coefficients for the internal places in order for eq. (6) to hold true.

A T-invariant (transition invariant) is a vector with the property that if each transition fires as many times as the value of the corresponding component of the vector indicates, the original marking is restored. Algebraically, these vectors are the solutions of eq. (6). Therefore, T-invariants correspond to flux distributions in steady state.

As Petri nets usually involve irreversible transitions only, all components of a T-invariant must be non-negative. T-invariants with this property are called true T-invariants. Frequently, the net direction of all biochemical reactions in a network is known, for example, because they are irreversible or have a defined biochemical function. In this case, the orientation of reactions can be chosen in such a way that all (net) fluxes are non-negative. Then, only steady-state flux distributions corresponding to true T-invariants are relevant.

A central concept in metabolic network analysis is that of elementary flux modes [Schuster and Hilgetag, 1994, Schuster et al., 2002a]. They stand for minimal sets of enzymes that could operate at steady state and are uniquely determined. That means that no other flux modes at steady state are proper subsets of the elementary modes. For a better understanding of the behaviour of biochemical systems, they can be decomposed into such simplest relevant routes. This has recently been demonstrated for sugar cane metabolism [Rohwer and Botha, 2001] and bacterial metabolism [Van Dien and Lidstrom, 2002]. In Petri net theory, elementary modes have, as counterpart, the minimal T-invariants [Starke, 1990]. The concept of elementary modes is, however, more general because reversible reactions are allowed. Colom and Silva [1990] developed an algorithm for computing minimal P-invariants. It can, by transposition of the incidence matrix, be used also for computing minimal T-invariants. It is based on row operations on the incidence matrix augmented with an identity matrix. In the course of calculation, care has to be taken to eliminate non-minimal and duplicate T-invariants. Colom and Silva [1990] propose two alternative tests to do so. A method for computing elementary flux modes based on convex analysis was proposed in [Pfeiffer et al., 1999; Schuster and Hilgetag, 1994; Schuster et al., 2000]. Although the latter algorithm was developed with different goals (Colom and Silva [1990] did not deal with metabolic networks) and completely independently of Petri net theory, the two algorithms show some similarities. However, they differ in that elementary modes can involve reversible reactions. This is taken into account by partitioning the stoichiometry matrix into "reversible" and "irreversible" submatrices. Moreover, the test for eliminating non-minimal and duplicate T-invariants (flux modes) is slightly different. For a more detailed comparison of the algorithms, see [Schuster et al., 2002a].

The T-invariants are helpful in studying several properties of Petri nets, such as consistency. This property means that there exists an initial marking and a corresponding firing sequence that regenerates the initial state and contains each transition at least once. As can be seen in the system shown in Fig. 8 with either reactions 3 and 4 completely inhibited or reactions 5 and 6 completely inhibited, not every metabolic system is consistent according to this definition. However, reactions 5 and 6 in the former case and reactions 3 and 4 in the latter case are not covered by true T-invariants. If we consider only a subnet that is covered by true T-invariants, such as transitions t₁ and t₂, it is consistent. This is because, once the minimal T-invariants (elementary flux modes) are identified, appropriate initial markings enabling these invariants to operate can be linearly combined and a new initial marking is obtained. The system can fire each above-mentioned T-invariant consecutively (the necessary resources exist due to the "construction" of the initial marking), each transition is used at least once and the initial marking is always regenerated.

A further property studied for Petri nets, reversibility, means that for every marking m that can be reached from M₀, M₀ can also be reached from M. It holds for metabolic networks, if some constraints are fulfilled. One constraint is that the network is covered by true T-invariants. The second constraint is that all external metabolites have enough tokens to operate all true T-invariants. The arguments read as follows: Let us denote the number of times the transitions ti have to fire in order to reach a marking m from M₀, by w_i. The numbers w_i are gathered in a vector W. Note that W need not fulfil eq. (6). As the net is covered by true T-invariants, we can find a vector V that does satisfy eq. (6) and a sufficiently large natural number such that V-W involves positive components only. This vector indicates how many times the transitions need to be fired to reach the initial marking again.

---

Siphons, traps, deadlocks and liveness

In Petri nets, special sets of places can be identified, for example, siphons, called also structural deadlocks, and traps [Reisig, 1985]. A siphon is a set of places that - once it is unmarked - remains so. A trap is a set of places that - once it is sufficiently marked - can never lose all its tokens. (It can happen that, if only some places of the trap are marked with a number of tokens smaller than a certain limit, the trap may lose all its tokens.) Clearly, any semi-positive P-invariant implies a trap because the total number of tokens is constant and can, hence, not reach zero. Moreover, superconservative subnets form traps, while subconservative subnets form siphons. The algorithms calculating siphons and traps [Schmidt, 1996a, 1996b; Yamauchi and Watanabe, 1999; Yamauchi et al., 1996; Tanimoto et al., 1996] in nets with specific properties are based on the following alternative definitions: A siphon is a set of places having the property that its input transitions set is contained in its output transitions set. A trap is a set of places for which its output transitions set is contained in its input transitions set. A Petri net N, having m₀ as initial marking, is said to be deadlock-free if for any reachable marking m, there is an enabled transition. A Petri net N, characterized by a current marking m is in deadlock if no transition is enabled to fire at marking m. Preventing deadlocks in an efficient way represents an intensely researched field [Kemper, 1993; Moody and Antsaklis, 1998; Huang et al., 2001; Chu and Xie, 1997; Varpaaniemi, 1993; Iordache et al., 2001].

Traps and structural deadlocks are interesting for biochemical modelling. Many biochemical networks have the function to produce storage substances in certain periods and consume these substances in other periods. For example, the potato plant produces starch and accumulates it in the potato tubers during growth, while starch is consumed after the tubers are deposited after the harvest. The starch and several of its precursors then form traps in the reaction net during growth, while starch and possible intermediates of degradation form siphons after the harvest. Consider the simple reaction system shown in Fig. 8. Transition t₁ is always activated while t₂ only fires if at least one token exists in S₁. Importantly, we consider t₃ and t₄ to be inoperative if t₅ and t₆ are operative in the system and vice versa. For example, the system could describe the production and degradation of starch. The internal metabolites then would be: S₁, glucose-1-phosphate, S₂, UDP-glucose, S₃, starch [Stryer, 1995]. In the starch example, it is not necessary to consider an intermediate S₄, while for other storage metabolites, it may be. In most cells containing starch, either the branch producing starch or the branch degrading it is functional. This is realized by complete inhibition of the appropriate enzymes. It can be easily observed that S₂ and S₃ form a trap when reactions 3 and 4 are operative. Once a token arrives in S₃, no transition able to fire exists in the system to consume this token, so it remains there independently of the later evolution of the system. In the other case, once the last tokens were extracted from S₃ and S₄, no transition able to generate a new token in these places exists, so they remain empty. This means that S₃ and S₄ form a siphon.

Current computer programs for simulating metabolic networks deal only partially with siphons and traps. For example, the program GEPASI developed by Pedro Mendes [1997] (http://gepasi.dbs.aber.ac.uk/softw/gepasi.html) detects all reactions that are at equilibrium in any steady state. For the example system shown in Fig. 8 with transitions 5 and 6 blocked, GEPASI would detect reactions 3 and 4 to have this status. Here, these reactions 3 and 4 are irreversible, so that no steady state can be reached. However, if the reactions were reversible, they would indeed attain thermodynamic equilibrium in any steady state. The Program METATOOL [Pfeiffer et al., 1999] (http://www.bioinf.mdc-berlin.de/projects/metabolic/metatool/) indicates metabolite S₃ as a "non-balanced internal metabolite", while it does not say anything about metabolites S₂ or S₄. It is worth including, in future refinements of simulation packages, routines for detecting all metabolites involved in traps or siphons.

Another important concept is liveness [Reisig, 1985]. A transition t is said to be live if, for any marking m reachable from m₀, there is a marking m' reachable from m, such that t is enabled by m'. A transition t is dead at marking m if no marking m' reachable from m enables t. A Petri net (N; m₀) is said to be live if every transition is live.

Importantly, deadlock-freeness and liveness are two different notions. Liveness means that all the system's transitions may be repeated infinitely often, while deadlock-freeness only implies that at least a subset of these transitions may be repeated, but not necessarily all of these. It is interesting how the concepts of liveness, deadlocks, siphons and traps are connected with each other. A net satisfies the deadlock-trap property if each non-empty siphon includes a trap and the maximal trap in each minimal deadlock is sufficiently marked. In this case, no dead marking is reachable. So, the net is deadlock-free. A further special class of Petri nets is made up of the free-choice nets. In these nets, each place has at most one output transition or the input places of the output transitions of p consist only of p for any place p belonging to the net. A free-choice net is live if and only if every non-empty siphon includes an initially marked trap. This property is also known as Commoner's theorem [Commoner, 1972]. Siphon and trap are dual notions. A siphon in a Petri net N is a trap of the net N' obtained by reversing direction of all edges of N. Therefore, the properties satisfied by siphons have counterparts for traps.

Liveness and deadlock-freeness are structural properties, in which the initial marking plays, however, an essential role. An example for a net that is in deadlock is the above example A + B gives 2 B (Fig. 5) if the initial token number of B is zero. The input and output transition sets of B coincide. Therefore, {B} is simultaneously a siphon and a trap. {B} = {A,B}  {B} (where {B} denotes the input places of the output transitions of B, that is, the co-substrates of B), but |A| = |{t}| = 1 and |B| = |{t}| = 1, so "A + B gives 2 B" is a free-choice net. Due to Commoner' theorem, this net cannot be live if the token number of A is not infinite (if A is not what we called external metabolite in Section 2) and trap {B} included in siphon {B} has not at least one token. Without tokens in B, this autocatalytic reaction then does not start proceeding. In chemistry, such a situation is known as false equilibrium [Othmer, 1981]. A larger example is glycolysis, which requires 2 moles of ATP in its upper part and produces 4 moles of ATP in its lower part. If no ATP is present at the beginning, the glycolytic pathway cannot proceed. Therefore, this pathway has been said to have a turbo design [Teusink, 1998]. A test for liveness and deadlock-freeness of the net can thus help us decide whether the metabolic system can attain a situation where it is blocked. The detection of siphons and traps is instrumental for this purpose.

---

Evaluating the role of TPI in Trypanosoma brucei metabolism by detecting siphons and traps

T. brucei is a unicellular, extracellular, eukaryotic parasite of the blood and tissue fluids of mammals. It is transmitted by tsetse flies and causes sleeping sickness in humans. The infections are lethal unless treated, but the few existing drugs have severe side-effects. Many studies [e. g. Bakker et al., 1999; Helfert et al., 2001] were focused on the carbon and free energy metabolism of this organism, which depends entirely on glycolysis. Accordingly to Scheme 1 in Helfert et al. [2001], glucose is imported into the glycosome and then converted into F-1,6-P. The two consecutive enzymes (HXK and PFK) are contracted in one step which consumes two ATP and produces two ADP. ALD converts F-1,6-P into DHAP and GA3P. These two substances are isomerised into each other by a reversible enzyme, TPI (triose-phosphate isomerase). DHAP is transformed into Gly3P by GPDH, with consumption of NADH and production of NAD⁺. GAPDH1 uses NAD⁺ to transform GA3P into BPGA and NADH. Gly3P can be either converted into glycerol by GLYK with consumption of ADP and production of ATP, or transported into the cytosol, where GPO oxidizes it to DHAP and H2O, DHAP being transported back to the glycosome. PGK uses an ADP molecule to convert BPGA in 3-PGA and ATP. 3-PGA is transported into the cytosol and converted, via 2-PGA, into PEP, which gives pyruvate and ATP, consuming one ADP. A similar system was treated by Overkamp et al. [2002]. They maximized the yield of glycerol in Saccharomyces cerevisiae using metabolic engineering.

At the beginning, Helfert et al. [2001] supposed that glycolysis could proceed without TPI, producing glycerol and pyruvate in the same amount. This corresponds to an elementary mode described in Schuster et al. [2000]. However, TPI knockout mutants turned out to be inviable. Thus, they built a kinetic model to explain the unexpected result that all system fluxes (PYR, GLYCEROL) decrease. We now give a structural explanation, thus ignoring the kinetics. In Fig. 9, the Petri net model corresponding to Scheme 1 in Helfert et al. [2001] is given. It should be noted that this network is not a free-choice net because, for example, |GLY3P| = |{t2, GLYK}|=2 and {GLY3P} = {GLY3P, ADP}  {GLY3P}. The first property is known in Petri net theory as conflict, because two transitions {t2, GLYK} compete for the same resources (the tokens from place B). But in metabolic networks, the token number is large enough that the transitions in competition will simply "agree" on the tokens distribution depending on their reaction rate. Taking this aspect into account, we do not need to know the reaction rates, but we only assume that the flux through transition T₂ in Fig. 9 is always greater than zero. Another important knowledge that we use is that ALD is reversible and therefore inhibited by its products [Stryer, 1995].

Let us consider the case when TPI is knocked out. t₁ = {NADH, NAD⁺} forms a siphon and a trap at the same time. Its input transitions set {GPDH, GAPDH} coincides with its output transitions. This means that once this set of places is sufficiently marked, it keeps its tokens. Moreover, {GPDH, GAPDH} forms also a P-invariant, their tokens sum remaining constant during the whole process. Another trap (T2) consists of {DHAPc, DHAPg, GLY3Pc, GLY3Pg, Gly} because its input transitions set {ALD, GPDH, GPO, GLYK, t₁, t₂} includes its output transitions set {GPDH, GPO, GLYK, t₁, t₂}. If the flux through t₂ were equal to zero, Gly would accumulate, but because the flux through t₂ is greater than zero, GLY3P is partially transformed back into DHAPg. Let us start proceeding with the marking m0. Following the firing sequence {e₃, ALD, GPDH, GAPDH, t₂, GPO, t₁, GPDH, t₂, GPO, t₁, PGK, t₃, e₁, e₂} as Table 2 illustrates, the network reaches a dead marking m1, because DHAP accumulates - no NADH being available for further converting it. This is because GPO is draining the flux, consuming NADH faster than GAPDH can produce it. This continues until the product inhibition of ALD is so strong that ALD ceases operating. Therefore, the whole system is dead, no more Gly and Pyr being produced. For deriving this result, it is informative that after deleting TPI, the three transitions t₁, t₂ and GPO are not involved in any elementary mode (minimal T-invariants) anymore.

The importance of TPI can be seen if its corresponding transitions are added in the model. TPI being a reversible enzyme, two transitions, one acting forwards and one - backwards, have to be added. Of course, the T-invariant {TPI1, TPI2} has to be ignored, having no biological significance. In this new context, T₂ is not any more a trap. Another minimal trap occurs in t₃ = {DHAPC, DHAPg, GLY3Pc, GLY3Pg, GLY, GA3P, BPGA, 3PGK, 2PGK, PEP, PYR} corresponding to the accumulation of GLY and PYR. Whenever DHAPg tends to accumulate, due to NADH lack, TPI1 converts a part of DHAPg tokens into GA3P. Due to the sufficient amount of NAD⁺, GAPDH fires with production of the necessary NADH, which gives GPDH the possibility to fire. Also if NAD⁺ is deficient, but GA3P is sufficient, TPI2 converts GA3P into DHAPg, GPDH fires and produces the required NAD⁺. Taking into account almost only structural properties of the given network, especially the presence of traps, we could evaluate the role of TPI in glycolysis and glycerol production. In the next section an example taken from nucleotide metabolism will be used to illustrate the notions presented above. The program INA (www.informatik.hu-berlin.de/~starke/ina.html) is utilized to facilitate the calculations.

---

Illustration of Petri net properties on a system extracted from nucleotide metabolism

Let us now consider the biochemical system depicted in Fig. 7A. It represents part of nucleotide metabolism, as it occurs, for example, in human liver [Stryer, 1995]. We have translated it in terms of a Petri net (Fig. 7B) and then, analysed it using the program INA. For modelling the external metabolites, we have chosen to introduce a self-loop for each source and each sink and to keep the same firing rule (2) independently of the metabolite's type. If we do not impose capacities for the internal metabolites, the net is unbounded because uridine, for example, can accumulate more and more tokens if Cdd keeps firing while Urk1 is not. Accordingly, the reachability tree is infinite. As INA has reported, the net is strongly connected, not pure (which is obviously due to the introduced self-loops), and not (sub-)conservative. There is no P-invariant, except the external metabolites on their own. Again due to the self-loops next to the external metabolites, the number of tokens in each of these places remains constant.

INA reports four minimal semi-positive T-invariants. We give them here by indicating the transitions with positive components in the vectors representing these T-invariants:

1. Urk2, Kcy2
2. Cdd, Urk1, Kcy1
3. 2 Kcy1, KPR, 2 UPP, APT
4. Kcy1, KAD, KPR, UPP.

Note that some transitions (such as Kcy1 and UPP in the third invariant) have to fire twice, but not necessarily successively. One can see that firing all the activated transitions that belong to an invariant regenerates the initial marking. As each enzyme occurs in at least one minimal T-invariant, the net is covered by these invariants. Therefore, the net is persistent and live. For simplicity's sake, although in biological organisms the reactions KAD and KPR are reversible, we considered them irreversible. If they are treated as reversible, care has to be taken that extra, irrelevant T-invariants, containing only KAD and KAD', and KPR and KPR' respectively, result (where the primed symbols denote the reverse reactions). They have to be discarded.

The biochemical meaning of the minimal T-invariants can be explained as follows: (1) production of cytidine-diphosphate (CDP) from cytidine, (2) production of uridine-diphosphate (UDP) from cytidine, (3, 4) two invariants producing uridine-diphosphate (UDP) from uracil in different ways. In the invariant (3), one mole of adenine per two moles of UDP produced is formed as a by-product. This is because ATP is used as a source of the ribose moiety, which is necessary for forming UDP from uracil. Note that this invariant is not easy to determine by inspection. Moreover, it can be seen that the molar yield with respect to ATP is different for the pathways (3) and (4). While invariant (3) consumes 3 moles of ATP per mole of UDP produced, invariant (4) uses 3 moles of ATP per two moles of UDP [Schuster et al., 2002a]. All of these T-invariants correspond to the so-called salvage pathways, which serve to save nucleotides from leaving the cell and redirect them to nucleotide phosphates [Stryer, 1995].

Let us now assume that ATP and ADP are internal metabolites and that the two enzymes KPR and KAD are not expressed in a certain cell type. If we modify the network by eliminating the transitions that stand for these enzymes (and also the arcs that connect them with their neighbouring places), we do obtain a P-invariant. It can be translated in terms of a conservation relation: ATP + ADP = const. The constant would be 2 if we define the initial token numbers of ATP and ADP to be 1 each. With any conservation sum less than four, the four remaining transitions consuming ATP (Urk1, Urk2, Kcy1, and Kcy2) are in mutual exclusion. Since the net does not include transitions producing ATP or (in the second case) AMP, these substances are eventually running out, so that the places standing for ATP and AMP are deadlocks, while in the complete net, there is neither a trap nor a deadlock. To maintain a steady state, nucleotide metabolism requires permanent production of ATP, for example, by glycolysis.

---

Discussion

Petri nets provide a special formalism to describe processes in networks. In particular, they are suitable to model biochemical networks. Here, we have shown that several concepts from Petri net theory have a significance for this modelling. However, there are alternative formalisms, and it is difficult to decide which formalism is best suited. To implement the calculations on computer, one usually translates Petri nets into matrices. So one may argue that the networks could be modelled by matrices from the very beginning. Indeed, Petri nets have the advantage to provide a means of visualisation. On the other hand, biochemists use a special way of visualisation for decades [Stryer, 1995; Kanehisa and Goto, 2000] (and chemists already for centuries). Multimolecular reactions such as "A + B gives C + D + E" are represented by an arrow that has two upper ends and three lower ends. This arrow can be represented, in a formal language, as a pair of n-tuples: ((A,B), (C,D,E)). In contrast, in a normal graph as used in graph theory, edges correspond to simple pairs of nodes. In Petri nets, the representation is "disentangled" by introducing additional nodes and arcs. The above reaction would then be represented by five place nodes and one transition node (T) linked by five arcs. The arcs correspond to the following pairs of nodes: (A,T), (B,T), (T,C), (T,D) and (T,E). It is a matter of taste which representation is preferred - one pair of n-tuples or several pairs.

Many concepts from Petri net theory have counterparts in traditional biochemical modelling, for example, P-invariants (conservation relations), T-invariants (flux modes), and minimal T-invariants (elementary flux modes). In metabolism, minimal T-invariants can be interpreted as biochemical pathways. Detection of these in complex networks is often not straightforward. It is helpful in determining maximal conversion yields [Rohwer and Botha, 2001; Schuster et al., 2000; Van Dien and Lidstrom, 2002].

The concepts of trap, siphon, deadlock, and liveness, among others, have not been considered in biochemical modelling so far. Here, we have shown that these are helpful to characterize special properties of metabolic networks. For example, the test for deadlock-freeness helps to determine whether a biochemical pathway can attain a false equilibrium, where it is blocked. From another point of view, this situation has been referred to as the danger of a turbo design of pathways [Teusink, 1998]. The liveness of a system indicates that all transitions are able to fire infinitely often, and the processes are not eventually restricted to a subsystem. Traps can correspond to storage metabolites that are produced during growth of an organism and steadily increase in their concentrations.

We have here analysed the example of the energy metabolism in t. brucei. If the accumulations in the trap exceed a certain amount, this can cause product inhibition of some transitions (the aldolase reaction in the example), forcing the system to stop working. This result is of interest for elementary-modes analysis. It has been argued that this analysis can help assert the effects of enzyme deficiencies and knockout mutations [Klamt and Stelling, 2003; Schuster et al., 2000]. The example analysed here shows that in a deficient system, the remaining elementary modes may not be functional because of occurrence of a trap. Therefore, pathway analysis should be refined by considering traps, siphons, deadlock-freeness and liveness.

Siphons can correspond to storage substances when they are gradually depleted during starvation. An analysis of traps and siphons appears to be promising in studying diseases such as obesity and hypercholesterolemia, which are related to over-accumulation of storage substances. It will be worth including the analysis of traps, siphons, deadlocks, and liveness in metabolic simulation packages.

So far, in Petri net theory, transitions are always considered to be unidirectional. However, many biochemical reactions such as all isomerases are known to be reversible in that their net flow can change sign depending on the physiological state. If such a reaction is described by two oppositely directed transitions, meaningless T-invariants arise. For example, in the scheme shown in Fig. 6, the T-invariant {T₁, t₂} occurs. In order to avoid the cancellation of such T-invariants after their computation, it will be worthwhile extending Petri net theory by allowing for reversible transitions.

A property that can be checked for Petri nets is boundedness. As biochemical networks are open systems, they are not usually covered by semi-positive P-invariants; that is, they are not conservative. Nevertheless, subnets are often covered by such invariants and are, therefore, bounded. For example, if the conservation relation ATP + ADP = const. holds, one can deduce that the energy currency metabolite ATP cannot exceed a certain limit. If negative coefficients exist in the conservation relation, boundedness cannot be guaranteed even for the corresponding subnet. Beside conservative subnets, there may be superconservative subnets. Obviously, they imply unboundedness. First, there may be metabolites that are only produced by irreversible reactions but not consumed by any reaction (Fig. 8). Second, if consuming reactions exist, the catalysing enzymes may have such a low maximal velocity (saturation level) that the rate of production is higher than the rate of consumption.

Here, we have focussed on topological analysis, which deals with the properties that occur from the static construction of the network. For many biological applications, such as the assignment of the metabolic function to an enzyme gene (functional genomics) [Dandekar et al., 1999; Förster et al., 2002; Selkov et al., 1997], it is sufficient to analyse these properties rather than the dynamics. The structural properties are the most representative features that one should look for. Compared to kinetic parameters of enzymes, they are constant in time and often much better known. Thus, reaction stoichiometries are easier to get from databases [Kanehisa and Goto, 2000; Selkov et al., 1996]. Topological analysis, (in particular, the computation of invariants) constitutes the basis for the simulation of the dynamics of the system.