 Methodology Article
 Open Access
 Published:
Manatee invariants reveal functional pathways in signaling networks
BMC Systems Biology volume 11, Article number: 72 (2017)
Abstract
Background
Signal transduction pathways are important cellular processes to maintain the cell’s integrity. Their imbalance can cause severe pathologies. As signal transduction pathways feature complex regulations, they form intertwined networks. Mathematical models aim to capture their regulatory logic and allow an unbiased analysis of robustness and vulnerability of the signaling network. Pathway detection is yet a challenge for the analysis of signaling networks in the field of systems biology. A rigorous mathematical formalism is lacking to identify all possible signal flows in a network model.
Results
In this paper, we introduce the concept of Manatee invariants for the analysis of signal transduction networks. We present an algorithm for the characterization of the combinatorial diversity of signal flows, e.g., from signal reception to cellular response. We demonstrate the concept for a small model of the TNFR1mediated NF κB signaling pathway. Manatee invariants reveal all possible signal flows in the network. Further, we show the application of Manatee invariants for in silico knockout experiments. Here, we illustrate the biological relevance of the concept.
Conclusions
The proposed mathematical framework reveals the entire variety of signal flows in models of signaling systems, including cyclic regulations. Thereby, Manatee invariants allow for the analysis of robustness and vulnerability of signaling networks. The application to further analyses such as for in silico knockout was shown. The new framework of Manatee invariants contributes to an advanced examination of signaling systems.
Background
Living cells interact with their environment to adapt to changes and perturbations. The cell needs to control pivotal decisions such as cell differentiation, proliferation or cell death to provide tissue homeostasis. Signal transduction processes mediate these cellular responses to changes of environmental and intracellular conditions. The appropriate response is orchestrated and processed by a highly intertwined and complex network of signal transduction pathways. Alterations of these processes can result in severe diseases and pathologies [1].
Mathematical models are useful frameworks to capture signaling systems and to gain insights of the systemwide processes. To elucidate and understand the complex regulatory mechanisms of signal transduction networks, various computational models have been applied. The models range from small, but detailed, kinetic models to abstract, largescale models, which exhibit a coarsegrained view of the cellular processes. Predominantly, ordinary differential equation (ODE)based kinetic modeling has been applied to quantitatively describe the changes of the species concentrations in time and the resulting substance flow. As the availability of kinetic parameters is often limited, ODEbased modeling is restricted to small and wellcharacterized processes, such as the I κBNF κB signaling module [2]. The majority of experiments on signaling pathways provides mainly qualitative information about the interrelations between the diverse molecular components of a cell. Qualitative data as, e.g., from knockout experiments encourage the application of alternative, topologybased modeling strategies. The behavior of a signaling network can be examined on an abstract level by the application of methods such as Boolean networks [3], interaction graphs [4] and Petri nets (PN) [5]. A Boolean network describes the interdependencies of the entities of a system by a set of Boolean functions. An interaction graph analyzes the network’s topology to reveal the architecture and regulatory principles of the signaling system. The Petri net formalism has been successfully applied to model signal transduction systems, for example, in the human iron homeostasis process [6], in the gene regulation of the Duchenne muscular dystrophy [7], in the Von HippelLindau tumor suppressor interaction network [8] and in the insulin receptor recycling [9, 10].
To investigate the robustness and vulnerability of a signaling network, all possible pathways, i.e., signal flows from signal reception to cellular response, need to be determined. The wellestablished concept of elementary modes defines subprocesses in metabolic networks that occur under steadystate conditions [11]. Within the PN formalism elementary modes coincide with minimal, semipositive transition invariants (TI). However, in the analysis of steadystate processes for signaling networks, there are crucial differences that need to be considered for the application of elementary modes or TI. Metabolic systems occur in homeostasis, while signaling pathways display a rather transient and timedependent behavior. Klamt et al. [4], Behre and Schuster [12] and Schuster and Junker [13] have successfully applied the steadystate assumption to signaling networks as well.
A necessary condition for signal transduction models is to reproduce experimentally measured knockout behaviour. In silico knockout experiments can predict the model behavior to perturbations, such as gene knockouts or inhibitions for drug treatment. A common concept to examine in silico knockouts for mathematical models is based on the computation of TI to reveal correct pathway dependencies. The in silico knockout matrix provides a valuable visualization of the effects of knockouts for signaling systems [14].
Signaling systems exhibit specific characteristics of signal propagation, such as inhibitory relations, crosstalks, feedback loops and signal amplifications. Signal amplification is an important mechanism, since the signal strength is enhanced in signaling cascades. For example, a kinase may have several substrates and thereby catalyzes phosphorylation reactions repeatedly until it is deactivated again. Another pivotal mechanism of signaling is the feedback loop, which influences upstream signaling processes either in a positive or negative way. To gain insights of the biological system, these specific regulatory features need to be captured in the mathematical model. As a consequence, the network topology contains directed cycles that hamper straightforward application of TI for the detection of complete signaling pathways.
To explore the basic behavior and perturbations of signaling systems, many approaches have been proposed to address the challenging issue of pathway analysis. Various adaptations of TI and elementary mode analysis have been introduced. A standard approach for PN is to substitute bidirectional read arcs by unidirectional arcs to represent the main direction of the signal flow [15, 16]. Bidirectional arcs represent cycles in the network and thereby affect the TI analysis. The modification simplifies the description of the molecular processes and may ignore important regulatory features. Sackmann et al. [17] have introduced the notion of feasible TI to find complete pathways in signaling systems. Feasible TI are specific linear combinations of TI determined via read arcs adjacent to place invariants (PI). The concept of feasible TI is applicable to reestablish a read arc linkage. Nevertheless, an algorithm for the computation of feasible TI is lacking. Behre and Schuster [12] have adapted the concept of elementary flux modes to signaling routes in enzyme cascades, in particular for systems consisting of phosphorylation and dephosphorylation cascades. Klamt et al. [4] have introduced interaction graphs and logical interaction hypergraphs to analyze the structure of signaling and regulatory networks. For interaction graphs, the detection of paths between pairs of species is equivalent to the computation of elementary modes.
Furthermore, alternative, graphbased approaches to analyze signaling pathways have been proposed as well, such as the theoretical framework for detecting signal transfer routes by ZevedeiOancea and Schuster [18]. This framework applies graphtheoretical breadthfirst search to reveal routes in a network from a specific initial factor to a determined target and vice versa. The routes detect dependencies of specific factors on targets. A related approach for signaling networks, which are represented as signed hypergraphs and consider composite nodes, has been presented by Wang and Albert [19]. The introduced method of elementary signaling modes is an extension of the simple path analysis and represents a counterpart to elementary modes. All previous concepts fall short to reveal the complete combinatorial diversity of signal flows at steady state in models of signaling systems, which usually exhibit feedback loops or amplification cycles.
In this article, we adapt the notion of feasible TI [17] and introduce the concept of Manatee invariants (MI) to detect complete pathways from signal reception to cellular response. For an introduction to the terminology of the PN formalism, we refer to the methods section. In the results section, we define the concept of MI to obtain feasibility for TI and give the formal definition of MI. Further, we apply the theoretical framework to the TNFR1mediated NF κB pathway and demonstrate the applicability of the concept for in silico knockout experiments. In the supplement, we describe an algorithm for the construction of MI.
Methods
In this section, we introduce all terms of the Petri net formalism used in the study based on [20–22].
Petri nets
A Petri net (PN) is a quintuple N=(P,T,F,W,m _{0}) with:

P and T are finite and disjunct sets of places and transitions, respectively.

F⊆(P×T)∪(T×P) is a set of arcs.

\(W: F \rightarrow \mathbb {N} \) defines the weight of each arc.

\(m_{0}: P \rightarrow \mathbb {N}_{0}\) is the initial marking.
The directed arcs define preplaces and postplaces for each transition. For a transition t∈T, the set F t={p∈P∣(p,t)∈F} denotes its preplaces and the set t F={p∈P∣(t,p)∈F} its postplaces. The sets F p={t∈T∣(t,p)∈F} and p F={t∈T∣(p,t)∈F} define pretransitions and posttransitions of a place p∈P, respectively. A transition is defined as an output transition if the transition has only preplaces and no postplaces. Analogously, a transition is defined as an input transition if it has only postplaces and no preplaces.
A read arc is a bidirectional arc and denotes a relation of a transition to a preplace that is also its postplace. A PN is pure if no transition exists, for which a preplace is also a postplace, i.e., without read arcs. A PN is ordinary if the weight of each arc is one.
Tokens are movable objects, which can be assigned to places. The token distribution is the marking, m, of a PN and defines a state of the system. The initial state of a system is given by the initial marking, m _{0}. For every state of a system, a place p carries an amount of tokens m(p)≥0. A transition t∈T has concession, i.e., it is enabled or activated, in a marking m if m(p)≥w(p,t) applies for each preplace, p∈F t, and the corresponding weights of the arcs, w(p,t)∈W. An enabled transition t _{ j } fires by moving tokens from its pre to postplaces according to m _{ new }(p _{ i })=m _{ old }(p _{ i })+c _{ ij }. In the incidence matrix C, the element c _{ ij } indicates the rearrangement of the tokens on a place p _{ i } if transition t _{ j } fires:
Feasibility
A firing of a sequence of transitions σ=(t _{1},t _{2},…,t _{ n }) results in a shift of tokens \(\Delta m:P\to \mathbb {Z}_{0}\) with
The vector \(x : T \to \mathbb {N}_{0}\) gives the number of occurrences of transition t _{ k } in the firing sequence σ by the component x _{ k }=# t _{ k }. In PN, the vector x is the Parikh vector of the sequence σ and is denoted by \(\overline {\sigma }\).
A sequence of transitions σ=(t _{1},t _{2},…,t _{ n }) is called to originate from a marking, m, if the marking enables the firing of each transition of the sequence. The pth transition, t _{ p }, of the sequence, σ, must have concession after firing the transitions of the prefix (t _{1},t _{2},…,t _{ p−1}) of sequence σ. This condition has to be true for each p=1,2,…,n. The Parikh vector, \(\overline {\sigma }\), of a sequence σ that originates from the initial marking m _{0} is called feasible in the marking m _{0}. The set L _{ N }(m) collects all firing sequences that originate from a marking m of a PN, N. The Parikh vector of each sequence σ=L _{ N }(m) is feasible in the marking m. A marking m ^{′} is reachable from a marking m _{0} if a sequence σ∈L _{ N }(m _{0}) exists with \(m_{0}\overset {\sigma }{\rightarrow }m'\). A Parikh vector x is called realizable in the marking m if a marking m ^{′} is reachable that makes x feasible, i.e., a sequence σ∈L _{ N }(m ^{′}) exists such that \(x=\overline {\sigma }\). A Parikh vector x that is not feasible in a marking m _{0} may become feasible for a reachable marking m ^{′}. In this case, x is realizable, but not feasible, for the marking m _{0}. A Parikh vector that is feasible in the marking m _{0} is always realizable in the marking m _{0} as well.
Invariants
A transition invariant (TI) of a PN is defined as a Parikh vector \(x:T\to \mathbb {N}_{0}\) that fulfills the equation
A TI is denoted as true or semipositive if x has no negative component, x≥0. The set of transitions, whose corresponding components in x are positive, is called support of x and is denoted by s u p p(x). If the greatest common divisor of its nonnull elements is one, a TI is called canonical. The set of canonical TI of a PN can be infinite. To obtain finite sets, we introduce the set of minimal TI. A TI, x, is minimal if no other TI, x ^{′}, exists such that s u p p(x ^{′})⊆s u p p(x). Any TI can be generated from linear combinations of minimal TI. The set of minimal TI is unique and finite. In the following, we consider minimal, semipositive TI and call them TI. A PN is covered by TI (CTI) if each transition occurs in at least one TI. For a Parikh vector x of a sequence σ that fulfills Eq. (2), we define realizable TI and feasible TI. A Parikh vector x is a realizable TI in N if a sequence σ∈L _{ N }(m ^{′}) exists with \(\overline {\sigma } = x\) and a reachable marking m ^{′} with \(m'\overset {\sigma }{\rightarrow }m'\) [21]. The term feasible TI is a special case of realizable TI [17]. A Parikh vector x is a feasible TI in N if a sequence σ∈L _{ N }(m _{0}) exists with \(\overline {\sigma } = x\) and \(m_{0}\overset {\sigma }{\rightarrow }m_{0}\).
Equivalently, the equation
defines minimal place invariants (PI). A PI, \({y\in \mathbb {N}_{0}^{P}, y \geq 0}\), characterizes a token conservation rule for a set of places. The set of places, whose corresponding components in y are positive, is called support of y and is denoted by s u p p(y). For further definitions and applications of PN theory, we refer to [5, 20–23] and for algorithms for the computation of invariants and the reachability problem, see [24–26].
Results and discussion
Concept of Manatee invariants
Signal transduction systems exhibit cyclic regulations, such as feedback loops, which cause cyclic structures in network models. A common example of a cycle in signaling cascades emerges in amplifications such as enzymecatalyzed reactions, i.e., reactions that restore the enzyme’s activity. The catalytic reaction describes a recycling of the enzyme and may be represented by a reaction motif like the MichaelisMenten reaction scheme. Such a recycling step has an impact on the biological relevance of TI. The PN in Fig. 1a describes an enzymecatalyzed reaction motif. The network is free of PI (PIfree) and CTI, T I _{1} =(syn S, bin, rel, deg P) and T I _{2}=(syn E, deg E). T I _{1} describes the inflow of substrate S, its binding to the enzyme E, forming the intermediate E:S complex and the production of P under release of E. T I _{2} captures the synthesis and degradation of enzyme E.
The TI of an enzymecatalyzed reaction such as in Fig. 1a covers, on the one hand, the recycling of the enzyme, T I _{1}, and, on the other hand, its synthesis and degradation, T I _{2}. Transition syn E, which synthesizes enzyme E, is not part of T I _{1} because the coupling of transition syn E to T I _{1} would result in an accumulation of the enzyme. An accumulation contradicts the definition of TI to describe processes at steady state. In signaling systems, synthesis or activation of enzymes is usually regulated by upstream signaling events. The enzyme may need to be activated by upstream processes in the first place, before it is able to catalyze any downstream reactions. Both the upstream processes of enzyme activation and the catalyzed downstream reaction may be covered by TI. However, due to the condition of minimality inherent to TI, the TI that covers the upstream processes and the TI of the downstream reaction remain uncoupled. This is also the case for the small example in Fig. 1a. T I _{2}, which captures the upstream process of the synthesis of E, remains uncoupled from the downstream process covered by T I _{1}, which is dependent on E.
The small network in Fig. 1a demonstrates an interrelation of TI that cannot be ignored for the analysis of signal transduction networks. Biologically, both TI are interrelated, since the functionality of T I _{1} depends on the synthesis of E in T I _{2}. These findings for TI of enzymecatalyzed reactions also occur in other cyclic regulations, such as feedback loops. Consequently, for networks with cyclic structures, TI are unable to discover all pathways in terms of sequences of reactions in a network, which correspond to independent processes of signal propagation in a signaling system. The coupling of biologically interrelated TI that represent complete pathways motivated the concept of Manatee invariants (MI).
We adapted the notion of feasible TI of Sackmann et al. [17]. TI that meet the property to be feasible in the initial marking describe an entire biological pathway best assuming that the initial marking corresponds to the physiological state of an unstimulated cell. A TI is feasible in a given input condition or initial marking if all transitions of the TI can fire in sequential order. For a detailed definition of feasibility, see the Methods section. We propose MI as linear combinations of TI to attain feasibility.
The concept of MI is best explained for the small PN in Fig. 1a. The catalyzed conversion of the substrate S to the product P in T I _{1} is blocked when no enzyme E is present in the system. Consequently, T I _{1} is only feasible if at least one token is assigned to the place E in the initial marking. Alternatively, the enzyme E can be produced by firing of transition syn E of T I _{2}. The linear combination of T I _{1} and T I _{2} is therefore also feasible in the zero initial marking, i.e., no tokens are assigned to any places in the initial marking.
PI have an important impact on the feasibility of TI, since the associated places need to be provided with a sufficient amount of tokens. Let us consider the network in Fig. 1b. This network represents the T I _{1}induced network of the PN in Fig. 1a. The TIinduced network covers all places and arcs between the transitions of the TI. Whereas the complete network in Fig. 1a is PIfree, the T I _{1}induced network in Fig. 1b exhibits a PI, P I _{1}=(E, E:S complex), which describes the recycling of the enzyme E. The linear combination y _{1}=T I _{1}+T I _{2} induces a network that matches the network in Fig. 1a and therefore is PIfree.
We define MI as linear combinations of TI that fulfill the following property. The MIinduced network, given by its transitions, all their preplaces and all arcs inbetween, is free of PI, except for the PI of the complete PN. We call the MI pure if the MIinduced network is PIfree and impure if the network contains a PI of the complete PN. Therefore, the pure MI of the PN model in Fig. 1a are y _{1}=T I _{1}+T I _{2} and y _{2}=T I _{2}. Note that, T I _{2} is feasible in any initial marking, and the T I _{2}induced network is PIfree.
The construction of MI is motivated by the correlation between the feasibility of TI and the existence of PI in TIinduced networks. We assume a PN with a minimal marking, i.e., tokens are assigned only to places belonging to a PI. For PI within the TIinduced network, no transition of the TI can provide tokens to the places of the PI. Thus, posttransitions of the places of the PI are only enabled if the transitions not belonging to the TI, provide tokens to the places of the PI. Consequently, the TI is not feasible in the minimal initial marking. Therefore, to attain feasibility, the TI need to be coupled to the processes that provide tokens to the PI in the TIinduced network. MI interrelate the corresponding TI and thus constitute a network that is PIfree. Evidently, a PI within the TIinduced network that is simultaneously a PI of the complete PN cannot be dissolved by any coupling to other TI. We assume that PI of the complete PN were assigned with tokens in the initial marking, but any PI within the TIinduced network needs to be provided with tokens by upstream processes of other TI.
Definition of Manatee invariants
We extended the concept of TI and defined the MI within the PN formalism. We assumed a PN model of the molecular processes of a signal transduction system. These processes should be able to operate in a stable equilibrium. The resources that are required for signal transduction should either be produced by the system itself or have to be provided by the environment. In terms of PN theory, the network should be CTI, pure and PIfree. The following definitions and methods were worthwhile also for other cases, but were motivated by the advantages of their application to PN that are CTI, pure and PIfree.
Definition 1:TIinduced network.Let X _{ TI } be the set of TI of the PN, N=(P,T,F,W,m _{0}). For Y⊆X _{ TI }, the TIinduced network is given by N _{ Y }=(P ^{′},T ^{′},F ^{′},W,m _{0}) with

\(T' = \bigcup _{x \in Y} supp(x)\),

\(P' = \bigcup _{t \in T'} Ft\) and

F ^{′}=((P ^{′}×T ^{′})∪(T ^{′}×P ^{′})) ∩F.
The nodes of the TIinduced network, N _{ Y }, comprise the transitions of the support, s u p p(x), of the TI, x, in Y and all preplaces, Ft, of these transitions. The arcs of N _{ Y } are all arcs between the transitions and places in N _{ Y } that are also arcs in N.
Definition 2: PIfree TIset (TIset).For a PN, N, let X _{ TI } and X _{ PI } denote the set of TI and the set of PI, respectively. Let Y⊆X _{ TI } be a (sub)set of TI of N. N _{ Y } is the TIinduced network of Y, and Y _{ PI } denotes the set of PI of N _{ Y }. We call the set of TI, Y, a PIfree TIset iff
In the following, we will use the term TIset as a synonym for PIfree TIset. Apart from the PI of the original PN, the TIinduced network is free of additional PI. For a PN that is CTI, the set of all TI is a TIset because condition (4) is fulfilled with Y _{ PI }=X _{ PI }.
Definition 3: pure/impure TIset.Let Y be a TIset of N, and Y _{ PI } denotes the set of PI of N _{ Y }. We call the TIset Y of N pure, iff
and impure otherwise.
For a PIfree PN, the conditions (4) and (5) become equivalent. The number of TIsets, n, of a given PN is finite with an upper bound of \(\phantom {\dot {i}\!}n\le 2^{X_{TI}}1\), where X _{ TI } is the cardinality of the set of TI of the PN. We introduce the term of a minimal TIset as the minimal basis to construct TIsets as unions of such minimal TIsets.
Definition 4: minimal TIset.Let M be the set of all TIsets of a PN. A TIset, Y∈M, is minimal, iff
Note that, the union of any two TIsets is a TIset, and hence, it is sufficient to consider unions of two TIsets to prove that a minimal TIset is unique and not a union of a finite number of other TIsets.
Definition 5: Manatee invariant (MI).Let Y be a minimal TIset of a PN. An integer linear combination
with \(c_{x} \in \mathbb {N}^{+} \) is a Manatee invariant (of Y). We call the Manatee invariant, y, pure if the TIset Y is pure and impure otherwise.
Because an MI is a linear combination of minimal TI, an MI is not necessarily a minimal invariant. The number of MI is infinite. We define a finite set of MI as minimal MI.
Definition 6: minimal Manatee invariant.Let Y be a minimal TIset of a PN. An MI, \(y =\sum _{x\in Y}\ c_{x}\ x\), is minimal if either

(type a) ∀x∈Y:c _{ x }=1, or

(type b) y is feasible in the initial marking, m _{0}=0, and no other MI \(y' =\sum _{x\in Y}\ c'_{x}\ x\) with c ^{′}≤c is feasible in the initial marking.
Preferably, we would like to compute MI of type b, i.e., MI that are feasible in any initial marking. However, the computation of MI of type b may not be possible because the PN either is not PIfree or contains a motif that restrains the feasibility. Exemplarily, Additional file 1: Figure S2 depicts a pure MI that is not feasible in m _{0}=0, i.e., the MI is of type a. In the following, we consider only minimal MI. For a description of an algorithm for the construction of MI, we refer to Additional file 1.
Case study: The NF κB signaling pathway
We applied the theoretical framework of MI to a model of the TNFR1mediated NF κB pathway. Figure 2 illustrates the biological processes of the NF κB pathway. The full protein names are given in the list of abbreviations in the declaration section.
TNFR1 is a membrane receptor, which mediates immune responses triggered by the cytokine TNF α. When TNFR1 is activated, various adaptor proteins, TRADD, RIP1, TRAF2 and E3 ligases, cIAP1, cIAP2, are recruited and form the receptor signaling complex (RSC). The E3 ligases catalyze the formation of ubiquitin chains on RIP1 to form a scaffold for the effector kinases, TAK1 and IKK. The activated kinases, TAK1 and IKK, induce the activation of the transcription factor, NF κB, by targeting the inhibitor of the transcription factor, I κB, for degradation. Released NF κB can promote gene expression of proinflammatory proteins. The signaling pathway is negatively regulated mainly via two feedback loops that eventually terminate the signal transduction. On the one hand, the inhibitor, I κB, gets restored to terminate NF κB activity and on the other hand, deubiquitinase A20 promotes the dissociation of the RSC to return to the physiological state of an unstimulated cell. For reviews of the signaling processes, we refer to [27–29].
PN model
Figure 3 illustrates the PN model of the NF κB signaling pathway. The model comprises the main signal flow of receptor activation, RSC formation and transcription factor activation. The regulation of transcription factor activity, along with initiation of gene expression and coordinated termination by feedback loops were also incorporated in the model. As this model serves as a case study, we restricted the model to the internal regulatory processes of the signaling pathway. Crosstalks or alternative cellular responses were not considered. Additional file 1: Tables S1 and S2 list all places and transitions of the PN model and their biological meaning as proteins or protein complexes and reactions, respectively. We developed and analyzed the model applying the opensource software MonaLisa [30]. We provide the file of the PN model in Additional file 2. The MI were computed by the algorithm described in Additional file 1.
The PN is composed of 31 transitions, 29 places and 69 arcs. The PN is pure and ordinary. The initial marking provides a token on each of the places of the gene of A20, Gen_A20, and gene of I κB, Gen_I κB. Note that, Gen_A20 and Gen_I κB are places belonging to P I _{1} and P I _{2}, respectively, see Additional file 1: Table S3. The PI concerns parts of gene expression and corresponds to the conservation of the genes of A20 and I κB. The PN has four TI and fulfills the CTI property, see Additional file 1: Table S4.
The trivial T I _{1} is highlighted blue in Additional file 1: Figure S3. T I _{1} consists of two transitions describing the process of NF κB synthesis and degradation. T I _{1} is feasible in the initial marking, since transition SynNF κB assigns tokens to the place, NF κB, which are consumed by transition DegNF κB. Therefore, all transitions of T I _{1} are enabled.
T I _{2} captures the process of TNFR1 activation and RSC formation coupled with gene expression of the deubiquitinase A20 and dissociation of the RSC promoted by A20. Thereby, T I _{2} describes the processes of the A20 feedback loop. T I _{2} is highlighted red in Additional file 1: Figure S3. T I _{2} is not feasible in the initial marking because gene expression of A20 is blocked without transcription factor NF κB located in the nucleus. The initial marking does not assign any tokens to the place of nuclear NF κB, NF κB_n, and thus, neither transition B i n8 nor transition SynA20 can ever be enabled. The T I _{2}induced network has two PI, P I _{1} and a PI comprising the places NF κB_n and NF κB_n :Gen_A20. Only P I _{1} of the complete net is provided with tokens in the initial marking.
T I _{3} captures the processes of the I κB feedback loop with I κB gene expression, binding of I κB to NF κB and subsequent activation of NF κB by the RSC complex targeting I κB for degradation. In Additional file 1: Figure S3, T I _{3} is colored green. T I _{3} is not feasible in the initial marking because, e.g., the gene expression of I κB is blocked without NF κB located in the nucleus. Since the place NF κB_n is not provided with any tokens in the initial marking, transitions B i n9 and S y n I κ B are never enabled. The T I _{3}induced network exhibits four PI, which correspond to the activation and inhibition of NF κB in the cytosol, the recycling of nuclear NF κB, the recycling of the gene of I κB during gene expression of I κB and the release of RSC during NF κB activation, respectively. The initial marking assigns tokens for the gene of I κB, whereas three PI remain insufficiently provided with tokens.
Additional file 1: Figure S4 depicts the PN and highlights T I _{4}. T I _{4} captures the processes of NF κB regulation, i.e., the activation of NF κB by kinases of the RSC complex and I κB degradation along with NF κBdependent initiation of the gene expression of I κB, which gets restored in the cytosol and translocates into the nucleus to shuttle NF κB back into the cytosol. T I _{4} is not feasible in the initial marking, since, e.g., the activation of NF κB is dependent on the kinases, TAK1 and IKK, of the RSC, but the place RSC_ub:TAK1:IKK carries no tokens in the initial marking. The T I _{4}induced network has three PI, which correspond to the conservation of NF κB in the cytosol and the nucleus, the recycling of the gene of I κB and the recycling of the activated RSC, respectively.
All three TI that are not feasible, T I _{2}, T I _{3} and T I _{4}, exhibit cycles in their induced networks. The cycles emerge from feedback loops or amplifications, such as the amplification during gene expression initiation or the repeated activation of NF κB by the kinases of the RSC.
Manatee invariant analysis
We applied the concept of MI to the NF κB pathway model. The computation yielded three MI, see Additional file 1: Table S5. The trivial TI of synthesis and degradation of NF κB, T I _{1}, is the smallest MI, M I _{1}. The M I _{1}induced network is free of any PI and hence, M I _{1} is pure and feasible in the initial marking.
M I _{2} is composed of all four TI, T I _{1},T I _{2},T I _{3} and T I _{4}, and covers the complete PN. M I _{2} describes the complete TNFR1mediated NF κB signaling pathway, from receptor ligation to transcription factor activation and initiation of feedback loops that eventually terminate signal transduction and restore the physiological state of an unstimulated cell. The M I _{2}induced network has two PI, P I _{1} and P I _{2}, of the complete PN, hence M I _{3} is impure and feasible in the initial marking.
M I _{3} combines the trivial TI, T I _{1}, with two nontrivial TI, T I _{2} and T I _{4}. Additional file 1: Figure S5 highlights M I _{3}. Aside from the direct binding of I κB to NF κB in the cytosol (transition Bin12), the M I _{2}induced network describes the complete PN. M I _{3} covers the processes of TNFR1 stimulation and subsequent formation of RSC, which leads to the activation of NF κB and NF κBdependent gene expression to terminate signal transduction via two feedback loops. Alternatively to the direct formation of the inhibitory complex of NF κB and I κB in the cytosol, free NF κB translocates into the nucleus to initiate the gene expression of I κB. Subsequently, restored I κB translocates into the nucleus, binds to NF κB and shuttles the inhibitory complex into the cytosol. The M I _{2}induced network has two PI, P I _{1} and P I _{2}, and therefore, M I _{3} is impure and feasible in the initial marking.
The MI reflected the complete combinatorial diversity of signal flows of the PN model. M I _{2} and M I _{3} cover complete pathways of TNFR1 signal transduction and subsequent signal termination along with restoration of the initial physiological state. The two MI distinguish for the two processes to form the inhibitory complex of NF κB and I κB and describe distinct, biologically reasonable processes of NF κB activation. While TI are suitable to determine feedback loops or other cyclic regulatory features in networks, MI are necessary to detect all the pathways, in which these feedback loops are involved.
In silico knockouts
For further analysis of the NF κB pathway model, we examined in silico knockouts based on MI. We studied the perturbation of the system for the knockouts of all proteins that constitute the signaling pathway. Figure 4 shows the resulting in silico knockout matrix. To study the influence of a knockout on certain proteins or protein complexes, the in silico knockout analysis adds output transitions to the respective species [14]. Additional file 1: Figure S6 shows the PN model of the NF κB pathway with the required additional output transitions grayed out. For a detailed description of in silico knockout analysis, we refer to Scheidel et al. [14]. Each row of the matrix represents the knockout of a protein, e.g., by deletion of an input transition that represents the synthesis of the protein. A column represents a protein or protein complex that may be affected by the respective knockout. The entries of the matrix are either green or red. A green entry denotes a place that has a pretransition in an MI, whereas a red entry denotes a place in the PN that has no pretransition in an MI. Biologically, a green entry indicates that, e.g., a protein complex can still be formed despite the knockout, while a red entry indicates that the complex cannot be formed.
The rows of the in silico knockout matrix shown in Fig. 4 exhibit two to eleven red entries. The knockout of the deubiquitinase A20 yields two red entries for protein A20 and A20 associated to the RSC complex, A20 and RSC_ub:TAK1:IKK:A20, respectively. All other complexes of the network are unaffected by a knockout of A20, which is in accordance to the known biological behavior, since A20 deubiquitinates the ubiquitin chains within the RSC [31]. The rows of TNF and TNFR1 have the highest number of eleven red entries, followed by TRADD and NF κB with ten entries each. Therefore, the knockout results indicate a high impact of TNF α, TNFR1, TRADD and NF κB. The receptor TNFR1 and ligand TNF α have the strongest impact on the system’s behaviour as the pathway is triggered by the receptor binding [32]. The knockout of TRADD, the adaptor protein of the stimulated receptor, has a high impact on the downstream system’s behavior as well. Also NF κB, the target transcription factor of the signaling pathway, has a strong influence on the functionality of the network. The knockout of NF κB affects the initiation of the feedback loops of A20 and I κB, while the formation of the RSC is unaffected. In the network, the synthesis of I κB and A20 is dependent on NF κB activity, therefore, the knockout results match the pathway’s behavior.
The matrix in Fig. 4 is a valuable representation for the effects of knockouts. The knockout results were intuitive and in accordance with expected signaling pathway dependencies. For comprehensive models, the influence of proteins and the detection of dependencies within the pathway may become more complex. In these cases, the in silico knockout can verify whether the model can reproduce experimentally measured knockout behavior, and predictions of the in silico knockout can be tested in further experiments.
The validity of the results of the in silico knockout analysis correlates with the property of MI to describe complete signal flows. The advantage of the application of MI can be conveniently determined in comparison to an in silico knockout analysis based on TI. The asterisks in Fig. 4 indicate red entries that turn green for the corresponding analysis based on TI, see Additional file 1: Figure S7. The in silico knockout matrix based on TI distinguished for 26 entries from the MIbased matrix. Surprisingly, TIbased analysis did not observe any effects for the NF κB knockout, even NF κB itself was unaffected by the knockout of its synthesis. The essential role of the target trascription factor, NF κB, is well known, and many processes in the system are controlled by NF κB activity [27]. Furthermore, the knockout matrix for TI indicated that the knockout of the proteins that constitute the RSC had no effect on the activation of NF κB as well, e.g., the knockout of IKK had no effect for the formation of complex RSC_ub:TAK1:IKK:NF κB:I κB. However, the phosphorylation of I κB can be triggered only by the activated kinase, IKK, bound to the ubiquitinated RSC [33]. Therefore in the model, all RSC proteins, including the kinases TAK and IKK, were necessary to initiate the phosphorylation and degradation of I κB. The TIbased analysis was unable to detect these pathway dependencies properly, and the results were misleading concerning their biological interpretation.
For in silico knockouts, the dependencies and relations of pathway components need to be captured correctly. The signal flows of a network reflect these dependencies best. TI or MI that are feasible in the initial marking define signal flows in the network and capture the pathway dependencies most appropriately. For TI that were not feasible, such as T I _{4}, we obtained misleading results for the in silico knockout. T I _{4} in Additional file 1: Figure S4 covers processes related to NF κB regulation and exhibits cyclic network structures of feedback loops and amplifications. Considering T I _{4}, a knockout of transition SynNF κB did not affect the TI, and the knockout matrix based on this TI observed no effects for NF κB knockout. Analogously, M I _{3} depicted in Additional file 1: Figure S5, comprising T I _{1}, T I _{2} and T I _{4}, would be affected by the knockout of transition SynNF κB. The valuable concept of MI attains the feasiblility by combination of interrelated TI to capture pathway dependencies of cyclic regulations to upstream or downstream processes. Thereby, only MI were able to derive the correct effects for in silico knockouts.
We proposed the concept of MI as a new method to compute functional pathways in PN models of signaling networks, even for signaling systems with amplification cycles or feedback loops. We showed that the application of MI as a precursor for further analysis like in silico knockouts is beneficial for the investigation of signaling networks. The detection of functional pathways in network models is elementary for many other rigorous network analysis approaches as well. For example, also the examination of crosstalks is dependent on a correct signal flow detection to determine shared processes of different signaling pathways. To evaluate the correctness of a signal transduction model and for a profound investigation that allow to postulate new hypotheses about its dynamic behavior, the computation of all possible signal flows is of great advantage.
Since the concept of MI is based on TI, it takes into account the steadystate assumption and causal relations of the components to determine signal flows in a network model. The proposed algorithm reveals all possible signal flows, which need to be considered for network analysis. An alternative simulation approach would be able to find the identical variety of signal flows in a model at least in the limit of an infinite number of simulation runs. However, the mathematical approach reveals the minimal solutions of all possible signal flows in a model. The mathematical concept of TI is an established precursor for further analyses, such as, e.g., maximal common transition sets [5], minimal cut sets [34] and Tclusters [35], which can be easily applied to MI as well.
A limitation of the rigorous analysis of all pathways in terms of MI is the complexity of the computational task. In worst case, the computation of TI requires exponential space [26]. Since the computation of MI requires the computation of TI, the computation of MI is also at least EXPSPACEhard. The computation of MI depends on size, structure and complexity of the network and may become infeasible for some network models due to the combinatorial explosion of the search space.
Conclusions
Characteristic and intrinsic regulation motifs of signal transduction like amplification reactions or feedback loops cause cycles in the topology of a network model. These cycles hamper the straightforward application of TI analysis for the detection of all possible signal flows, since cyclic, minimal TI usually do not reflect the entire pathways. In this article, we introduced the concept of MI, which aims to detect all signal flows from signal reception to cellular response including cyclic regulations. We adapted the concept of feasible TI. MI combine interrelated TI that are disconnected due to cyclic network structures with the objective to attain feasibility. Specific linear combinations of TI interrelate cyclic regulations to linked upstream or downstream processes, reflecting all signal flows from signal reception to cellular response. We presented an algorithm for the construction of MI to compute the combinatorial diversity of pathways from causal dependencies of reactions in a model.
We demonstrated the applicability of the concept of MI for a PN model of the TNFR1mediated NF κB signaling pathway. Exemplarily, we elucidated the benefit of MI application for in silico knockout studies. MIbased knockouts revealed correct effects for all protein knockouts of the network, whereas a TIbased analysis failed to detect essential interdependencies of network components. We suggest that other network analysis techniques can also benefit from the concept of MI to obtain biologically relevant conclusions. We presented MI as a straightforward approach for the detection of signal flows to advance modeling and functional pathway analyses, in particular of signal transduction networks.
Abbreviations
 A20:

Tumor necrosis factor αinduced protein 3
 cIAP1/2:

Cellular inhibitor of apoptosis protein 1/2
 CTI:

Covered by transition invariants
 IKK:

I κB kinase complex
 IκB:

Inhibitor of nuclear factor κB
 MI:

Manatee invariant
 NF κB:

Nuclear factor κB
 NLS:

Nuclear localization signal
 PI:

Place invariant
 PN:

Petri net
 RIP1:

Receptorinteracting protein 1
 RSC:

Receptor signaling complex
 TAK1:

TGF βactivated kinase 1
 TI:

Transition invariant
 TNF α :

Tumor necrosis factor α
 TNFR1:

Tumor necrosis factor receptor 1
 TRADD:

TNF receptorassociated protein with a death domain
 TRAF2:

TNF receptorassociated factor 2
 Bin:

Binding
 Deg:

Degradation
 Dis:

Dissociation
 Exp:

Nuclear export
 Imp:

Nuclear import
 n:

Nuclear
 p:

Phosphorylated
 Rel:

Release
 Syn:

Synthesis
 Ub:

Ubiquitinated
 X:Y:

Complex of protein X and protein Y
References
 1
Downward J. The ins and outs of signalling. Nature. 2001; 411(6839):759–62.
 2
Hoffmann A, Levchenko A, Scott M, Baltimore D. The I κBNF κB Signaling Module: Temporal Control and Selective Gene Activation. Sience. 2002; 298(5596):1241–5.
 3
Schlatter R, Schmich K, Avalos Vizcarra I, Scheurich P, Sauter T, Borner C, Ederer M, Merfort I, Sawodny O. ON/OFF and Beyond  A Boolean Model of Apoptosis. PLoS Comput Biol. 2009; 5(12):e1000595.
 4
Klamt S, SaezRodriguez J, Lindquist JA, Simeoni L, Gilles ED. A methodology for the structural and functional analysis of signaling and regulatory networks. BMC Bioinforma. 2006; 7(1):56.
 5
Koch I, Reisig W, Schreiber F. Modeling in Systems Biology: The Petri Net Approach. Computational Biology, vol. 16. London Dordrecht Heidelberg New York: Springer; 2011.
 6
Sackmann A, Formanowicz D, Formanowicz P, Koch I, Blazewicz J. An analysis of the Petri net based model of the human body iron homeostasis process. Comput Biol Chem. 2007; 31(1):1–10.
 7
Grunwald S, Speer A, Ackermann J, Koch I. Petri net modelling of gene regulation of the Duchenne muscular dystrophy. Biosystems. 2008; 92(2):189–205.
 8
Minervini G, Panizzoni E, Giollo M, Masiero A, Ferrari C, Tosatto SC. Design and Analysis of a Petri Net Model of the Von HippelLindau (VHL) Tumor Suppressor Interaction Network. PloS ONE. 2014; 9(6):96986.
 9
Balazki P, Lindauer K, Einloft J, Ackermann J, Koch I. MONALISA for stochastic simulations of Petri net models of biochemical systems. BMC Bioinforma. 2015; 16(1):215.
 10
Scheidel J, Lindauer K, Ackermann J, Koch I. QuasiSteadyState Analysis based on Structural Modules and Timed Petri Net Predict System’s Dynamics: The Life Cycle of the Insulin Receptor. Metabolites. 2015; 5(4):766–93.
 11
Schuster S, Fell DA, Dandekar T. A general definition of metabolic pathways useful for systematic organization and analysis of complex metabolic networks. Nat Biotechnol. 2000; 18(3):326–32.
 12
Behre J, Schuster S. Modeling Signal Transduction in Enzyme Cascades with the Concept of Elementary Flux Modes. J Comput Biol. 2009; 16(6):829–44.
 13
Schuster S, Junker BH. Topological Analysis of Metabolic and Regulatory Networks. In: Modeling in Systems Biology. London Dordrecht Heidelberg New York: Springer: 2011. p. 209–24.
 14
Scheidel J, Amstein L, Ackermann J, Dikic I, Koch I. In Silico Knockout Studies of Xenophagic Capturing of Salmonella. PLoS Comput Biol. 2016; 12(12):e1005200.
 15
Rodriguez EM, Rudy A, del Rosario RC, Vollmar AM, Mendoza ER. A discrete Petri net model for cephalostatininduced apoptosis in leukemic cells. Nat Comput. 2011; 10(3):993–1015.
 16
Heiner M, Koch I. Petri Net Based Model Validation in Systems Biology In: Cortadella J, Reisig W, editors. Applications and Theory of Petri Nets 2004. Berlin Heidelberg: Springer: 2004. p. 216–37.
 17
Sackmann A, Heiner M, Koch I. Application of Petri net based analysis techniques to signal transduction pathways. BMC Bioinforma. 2006; 7(1):482.
 18
ZevedeiOancea I, Schuster S. A theoretical framework for detecting signal transfer routes in signalling networks. Comput Chem Eng. 2005; 29(3):597–617.
 19
Wang RS, Albert R. Elementary signaling modes predict the essentiality of signal transduction network components. BMC Syst Biol. 2011; 5(1):44.
 20
Reisig W. Petri Nets: An Introduction. EATCS Monographs on Theoretical Computer Science, vol. 4. Berlin Heidelberg New York: Springer; 1985.
 21
Starke PH. Analyse von PetriNetzModellen. Leitfäden und Monographien der Informatik, vol. 6. Stuttgart: Teubner; 1990.
 22
Murata T. Petri nets: Properties, analysis and applications. Proc IEEE. 1989; 77(4):541–80.
 23
Goss PJ, Peccoud J. Quantitative modeling of stochastic systems in molecular biology by using stochastic Petri nets. Proc Natl Acad Sci. 1998; 95(12):6750–5.
 24
Colom JM, Silva M. Convex geometry and semiflows in P/T nets. A comparative study of algorithms for computation of minimal psemiflows. LNCS. 1991; 483:79–112.
 25
Ackermann J, Einloft J, Nöthen J, Koch I. Reduction techniques for network validation in systems biology. J Theor Biol. 2012; 315:71–80.
 26
Lipton RJ. The reachability problem requires exponential space. Research report 62, Dept. of Computer Science, Yale University. 1976.
 27
Wertz IE, Dixit VM. Signaling to NF κB: regulation by ubiquitination. Cold Spring Harb Perspect Biol. 2010; 2(3):a003350.
 28
Hinz M, Scheidereit C. The I κB kinase complex in NF κB regulation and beyond. EMBO Rep. 2013; 15(1):46–61.
 29
Ruland J. Return to homeostasis: downregulation of NF κB responses. Nat Immunol. 2011; 12(8):709–14.
 30
Einloft J, Ackermann J, Nöthen J, Koch I. MonaLisa–visualization and analysis of functional modules in biochemical networks. Bioinformatics. 2013; 29:1467–70.
 31
Wertz IE, O’rourke KM, Zhou H, Eby M, Aravind L, Seshagiri S, Wu P, Wiesmann C, Baker R, Boone DL, et al. Deubiquitination and ubiquitin ligase domains of A20 downregulate NF κB signalling. Nature. 2004; 430(7000):694–9.
 32
Wajant H, Scheurich P. TNFR1induced activation of the classical NF κB pathway. FEBS J. 2011; 278(6):862–76.
 33
Karin M, BenNeriah Y. Phosphorylation meets ubiquitination: the control of NF κB activity. Annu Rev Immunol. 2000; 18(1):621–63.
 34
Klamt S. Generalized concept of minimal cut sets in biochemical networks. Biosystems. 2006; 83(2–3):233–47.
 35
GrafahrendBelau E, Schreiber F, Heiner M, Sackmann A, Junker BH, Grunwald S, Speer A, Winder K, Koch I. Modularization of biochemical networks based on classification of Petri net tinvariants. BMC Bioinforma. 2008; 9(1):90.
Acknowledgements
Not applicable.
Funding
This work is supported by the Cluster of Excellence ’Macromolecular Complexes’ of the GoetheUniversity Frankfurt (3212070002/TP2) and by the LOEWE program Ubiquitin Networks (UbNet) of the State of Hesse (Germany)(20120712/B4). The founders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.
Availability of data and materials
The datasets supporting the conclusions of this article are included within the article and its additional files.
Author information
Affiliations
Contributions
LA conceptualized the study and wrote the manuscript. JA implemented the algorithm and was a major contributor in writing the manuscript. IK was involved in writing the paper. JS analyzed the data for the knockout analysis. SF, ID, IK conceptualized and supervised the study. All authors read and approved the final version of the manuscript.
Corresponding author
Correspondence to Ina Koch.
Ethics declarations
Ethics approval and consent to participate
Not applicable.
Consent for publication
Not applicable.
Competing interests
The authors declare that they have no competing interests.
Publisher’s Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Additional files
Additional file 1
Algorithm for the construction of Manatee invariants, Supplementary Figures S2S7 and Tables S1S5. (PDF 1116 kb)
Additional file 2
File of the NF κB Petri net model. (MLPROJECT 68 kb)
Rights and permissions
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
About this article
Received
Accepted
Published
DOI
Keywords
 Signaling pathway
 Mathematical model
 Petri net
 Transition invariant
 Feasibility
 Manatee invariant
 NF κB pathway