Petri nets and ODEs as complementary methods for comprehensive analysis on an example of the ATM–p53–NF-κ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\kappa$$\end{document}B signaling pathways

Intracellular processes are cascades of biochemical reactions, triggered in response to various types of stimuli. Mathematical models describing their dynamics have become increasingly popular in recent years, as tools supporting experimental work in analysis of pathways and regulatory networks. Not only do they provide insights into general properties of these systems, but also help in specific tasks, such as search for drug molecular targets or treatment protocols. Different tools and methods are used to model complex biological systems. In this work, we focus on ordinary differential equations (ODEs) and Petri nets. We consider specific methods of analysis of such models, i.e., sensitivity analysis (SA) and significance analysis. So far, they have been applied separately, with different goals. In this paper, we show that they can complement each other, combining the sensitivity of ODE models and the significance analysis of Petri nets. The former is used to find parameters, whose change results in the greatest quantitative and qualitative changes in the model response, while the latter is a structural analysis and allows indicating the most important subprocesses in terms of information flow in Petri net. Ultimately, both methods facilitate finding the essential processes in a given signaling pathway or regulatory network and may be used to support medical therapy development. In the paper, the use of dual modeling is illustrated with an example of ATM/p53/NF-κ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\kappa$$\end{document}B pathway. Each method was applied to analyze this system, resulting in finding different subsets of important processes that might be prospective targets for changing this system behavior. While some of the processes were indicated in each of the approaches, others were found by one method only and would be missed if only that method was applied. This leads to the conclusion about the complementarity of the methods under investigation. The dual modeling approach of comprehensive structural and parametric analysis yields results that would not be possible if these two modeling approaches were applied separately. The combined approach, proposed in this paper, facilitates finding not only key processes, with which significant parameters are associated, but also significant modules, corresponding to subsystems of regulatory networks. The results provide broader insight into therapy targets in diseases in which the natural control of intracellular processes is disturbed, leading to the development of more effective therapies in medicine.


Models.
In this study, we analyze two models, describing the same signaling pathway but implemented using two different methods: ODEs and Petri nets. The ODEs model describing the mechanism of the DNA damage response system based on the ATM/p53/NF-κ B pathways 12 was used as the basis for manual developing a corresponding model based on the Petri nets. Classical Petri nets enable creation of a qualitative model based on interactions between components, without any data on temporal system responses.
"Translation" of the ODE model into a Petri net. Diagrams depicting biochemical reaction networks, whose kinetics is the basis for creation of ODE models, resemble, to some extent, Petri net model representation. However, transformation of the former into the latter is not necessarily straightforward. An example, illustrating a part of the network analyzed in this paper, is shown in Fig. 1a. In the figure, text blocks represent molecules taking part in the processes, such as synthesized transcripts, proteins and their complexes. Solid arrows represent reactions, including transition of molecules into a different (e.g., active) form. Dashed arrows indicate that the molecule given in the associated text box affects a reaction represented by the solid arrow. An arrow ending with a sharp point represents stimulation of a reaction by a given compound, while a blunt point means inhibition of a reaction. Figure 1b shows a corresponding Petri net model, where text blocks from ODE schemas were translated into places, and reactions with their regulatory mechanisms into transitions. It should be noted that a model based on a classical Petri net is qualitative. Therefore, translation of ODE model into this kind of Petri net involves retaining the structure of the model but quantitative information, such as values of kinetic reaction rates or algebraic expressions describing kinetic rates, is not taken into account. In classical Petri nets, different types of arcs are not distinguished, contrary to the arcs in the diagrams used for ODE models. In a model based on Petri nets different types of reactions, like the transition of a given compound into an active form, forcing a reaction, inhibition of a reaction or degradation are marked by a classic directed arc (for example, there are no inhibitor arcs). Due to these limitations, the proposed Petri net model had to be adapted, in particular to allow for including terms of inhibiting specific reactions, degradation as well as enzymatic reactions. The first limitation related to inhibition reactions and degradation was solved in the way described by Gutowska et al. 32 . This solution is based on creating an additional transition that prevents the firing of a selected reaction (i.e., transition) by taking tokens from its pre-places. That way, a specific reaction is weakened (inhibited). A complete inhibition of a reaction is realized through knockout analysis, described farther in the text. The second limitation is associated with enzymatic reactions. Enzymes that take part in a reaction should not be removed but should be available for next reaction. Therefore, pools of particular components have been created to ensure their return to circulation. For example, a transition t 158 corresponds in a biological context to the creation of a pool of PIP3.
To illustrate how the ODE model was manually translated into the Petri net model, let us consider one of its equations, describing change of the level of nuclear phosphorylated p53 protein ( p53 pn ), depicted in Fig. 1a: In Eq. (1), subsequent expressions describe, respectively (in the order in which they appear in the equation): spontaneous p53 activation, ATM-dependent and Chk2-dependent p53 activation, Wip1-dependent inactivation, spontaneous and Mdm2-dependent degradation. They correspond to particular fragments of the the Petri net (Fig. 1b): • Spontaneous p53 activation. Three elements of the Petri net are needed for modeling this subprocess: a place (p53n) preceding a transition (spontaneous phosphorylation of p53n), the transition itself, and a place (p53pn) succeeding the transition (spontaneous phosphorylation of p53n). Due to the fact that it is a spontaneous reaction, no additional component (an additional preceding place) is needed to stimulate the reaction.
(1)  www.nature.com/scientificreports/ • ATM-dependent and Chk2-dependent p53 activation: Both reactions occur independently of each other, therefore they are modeled as two independent subprocesses. The following elements of the Petri net are required for modeling this subprocess: two places (p53n and ATM) preceding the transition (phosphorylation of p53n by ATM), the transition itself and a place (p53pn) succeeding the transition (phosphorylation of p53n by ATM). Two components are required to fire the transition (phosphorylation of p53n by ATM). Chk2-dependent p53 activation is realized in the same way. • Wip1-dependent p53 inactivation: Activation/inactivation (phosphorylation/dephosphorylation) is modeled similarly. The only difference is the direction of arcs (activation: p53n → activation → p53pn, inactivation: p53n ← inactivation ← p53pn). The following elements of the Petri net are needed for modeling this subprocess: two places (p53pn and Wip1n) preceding transition (dephosphorylation of p53pn by Wip1n), this transition and place (p53n) succeeding the transition (dephosphorylation of p53pn by Wip1n). • Spontaneous and Mdm2-dependent p53 degradation: A transition corresponding to degradation, unlike other reactions, does not have places succeeding it (such a transition can be treated as an output transition). Two elements of the Petri net are required for modeling of spontaneous degradation: a place (p53pn) preceding transition (degradation of p53pn), and the transition itself. In the case of Mdm2-dependent degradation two places (p53pn and Mdm2pn) preceding transition (degradation of p53pn by Mdm2pn), and this transition are needed. In this case, explicit output transitions are required to prevent further tokens flow. These output transitions, in addition to represent degradation, facilitate covering the Petri net model with t-invariants, which is a prerequisite for conducting the analysis. • p53 production: Similarly, as in the case of the ODEs model (which is presented in Fig. 1a) as a wiring diagram), this subprocess is symbolically represented in Fig. 1b) by a directed arc.
The processes of inactive p53 protein production (gene activation, transcription and translation) are included in other equations, which is why they are marked symbolically to the left of the block diagram ( Fig. 1a) as an arrow.
The ODE-based model. An extensive description of the ODE-based model used in this paper can be found in Jonak et al. 12 . In this model, ordinary differential equations were used to describe the dynamics of 61 components, such as proteins and transcripts, taking into account different possible forms of the same protein (e.g., phosphorylated and unphosphorylated). The model contains 135 parameters describing both the biochemical reaction rates and the physical properties of the simulated cell (e.g., the cytoplasmic to nuclear volume ratio or the number of alleles of a given gene). Simulations were run using ode23tb function in Matlab, with parameter values taken from Jonak et al. 12 .
Using the implemented model, it is possible to perform simulations with various combinations of input signals. The first input signal represents tumor necrosis factor α (abbreviated as TNF farther in the text), which stimulates the canonical NF-κ B signaling pathway. The second one represents ionizing radiation, which leads to DNA damage and subsequent activation of the p53 signaling pathway as well as NF-κ B signaling pathway.
We defined two simulation protocols for further analysis. These protocols differ only in the total irradiation dose. In each protocol, simulations were divided into three consecutive phases: 1. Waiting for equilibrium phase-24 h phase during which the model reaches steady state (no IR, no TNF). 2. Stimulation phase-the phase in which signaling pathways are stimulated by the simultaneous activation of TNF (10 ng/ml) and IR (1 h irradiation with a total dose of 4 or 10 Gy). 3. After-irradiation phase-240 h phase during which the cell attempts to repair DNA damage induced in phase 2 (no IR, TNF still present at 10 ng/ml).
Sample time courses corresponding to these simulation protocols are shown in Fig. 2.
Model based on Petri nets. The proposed model concentrate on ATM, p53, NF-κ B pathways, and excitations associated with the presence of TNF and IR. This model contains 170 transitions (elementary processes) and 89 places (biological components). Furthermore, it is covered by 541 t-invariants. All place and transition names are listed in Supplementary Table S1 and S2 (available online in Supplementary information section). Figure 3 shows the proposed model, taking into account division into particular modules: ATM (green), p53 (gray), NF-κ B (red), WIP1 (blue), CREB (purple) and extortions of the system (yellow). The Petri net model was created using a tool called Holmes 33 , facilitating manual creation of this type of models. The described model is available for download in two different formats, using links given in the section Supplementary information. The first one, .project, is a Holmes-specific format, while the second one, .spped, is a format specific to Snoopy (another tool used to model Petri nets) 34 . Holmes supports both of these formats, and the manual for Holmes software is available at the link given by Radom et al. 33 .
Adequacy of Petri net model and ODE model was confirmed by comparing crucial checkpoints, obtained for different inputs setup. For that aim, knockouts of modules associated with particular inputs were used, i.e., the knockout of a single input module (TNF or IR) and the knockout of both input modules (TNF and IR). The latter one corresponds to a steady state. That way, the system response to a TNF activation can be analyzed by a knockout analysis of the IR module. Similarly, the system response to an IR activation can be analyzed by a knockout analysis of the TNF module. In the case of TNF activation a significant activation of the NF-κ B pathway, and modules closely associated with it, should be expected, which is confirmed by the results from a knockout of IR module in the proposed model. For IR excitation, significant activation of p53 pathway is   Sensitivity analysis of ODE-based model. In this paper, we used the one-at-a-time (OAT) method based on the frequency distribution of a model transient response described in the work of Kardynska and Smieja 35 . Following the procedure proposed in Kardynska and Smieja 35 numerous simulations with randomly changed parameters were carried out. For each parameter in the model, 1000 simulations were carried out in which the value of this parameter was randomized within the range of ± 30% of its nominal value using a uniform distribution. Only the kinetic parameters of the model (such as the rates of biochemical reactions) were analyzed, as only these may become possible molecular targets for drugs. Moreover, only the parameters related to specific biochemical reactions have their counterparts in the Petri net as transitions. As a result, the number of analyzed parameters was reduced to 123 (from the initial 135). These parameters are listed in Supplementary  Table S4. The algorithm described in the work of Kardynska and Smieja 35 suggests selecting one model variable (e.g., the concentration of one of the proteins present in the model) as a representation of the model response. In the case of the model ATM/p53/NF-κ B, which combined multiple signaling pathways, finding such a single variable proved very difficult. To solve this problem, we decided to take a different approach and define a combined model response, based on time courses of all variables in the model. First, we created parameter rankings for each model variable according to the procedure described in Kardynska and Smieja 35 . Then we calculated the arithmetic mean of the sensitivity indices for parameters from all ranking (for all model variables). The resulting parameter ranking represents the cumulative effect of changing a given parameter value on all variables in the model. It should be noted that a high position of a given parameter in the ranking may suggest that a change of this parameter value results in a substantial change in few model variables or a moderate change in many model variables. Parameters ranking does not make distinguishing between the two cases possible.
We found that many parameters substantially affect the model response (Fig. 4). Considering the structure of the model, which combines several signaling pathways, results were discussed separately for groups of parameters associated with the p53 and NF-κ B module, the DNA damage detection module (related to ATM protein), Wip1/KSRP module and parameters related to the process of cell cycle arrest and apoptosis (Bax/p21 proteins).

Analysis of 10 ng/ml TNF and 4 Gy IR simulation results.
In the simulation with a low dose of radiation, sensitivity analysis showed that the parameter with the greatest impact on the model response is a DSB damage rate caused by IR. That finding is not surprising. Unfortunately, it does not lead to any conclusions concerning prospective molecular drug targets.
Among the parameters of the p53 signaling module, many are particularly important, including those directly affecting the level of p53 protein (synthesis of p53 mRNA, p53 transcript degradation, p53 translation rate) and its inhibitor, Mdm2 (synthesis of Mdm2 mRNA, Mdm2 transcript degradation, Mdm2 translation rate, Mdm2p protein spontaneous degradation, Mdm2 degradation by Chk2).
Surprisingly, the entire NF-κ B signaling module turned out to be very resistant to changes in parameters (compared to the remaining signaling pathways combined in the analyzed model) and none of its parameters substantially affects the model response. This may result from the fact that the NF-κ B responds to the TNF or IR very quickly, and its response is also quickly suppressed, resulting in low levels of the proteins involved. Unlike them, proteins associated with the p53 signaling pathway persist at elevated levels for several days following pathway excitation.
Analysis of 10 ng/ml TNF and 10 Gy IR simulation results. In the simulation with a high dose of radiation, sensitivity analysis showed that the parameter with the greatest impact on the model response is Mdm2 translation rate. This parameter appeared in the previous analysis but was not the most important one. In fact, all parameters mentioned in the previous subsection proved to significantly affect the model response in the simulation with 10 Gy. In addition to them, parameters of processes forming a positive feedback loop in the p53 control system (Mdm2 and PTEN gene inactivation) and p53 degradation by Mdm2 were crucial for the model response.
Substantial differences in the parameter rankings were observed for the ATM module. The only two parameters found to be relevant to the model response in the simulation with a high dose of radiation are the rates of DNA damage creation (DSB damage rate caused by IR) and repair (DSB repair rate).
No changes in the significance of the parameters for modules NF-κ B, Wip1/KSRP, Bax and p21 have been observed.
Analysis of Petri net-based model. The analysis for the model based on Petri nets is based on significance analysis at the level of individual transitions and significance analysis at the level of transition subsets (searching for a certain set of transitions that is included in a sufficiently large number of t-invariant supports).
The significance analysis finds subprocesses that are the most crucial for the functioning of the system under investigation, i.e. subprocesses leading to DNA damage and repair processes. The results showed that the most important elementary processes in the proposed model are associated with creation of DSB by IR, ATM, p53, and NF-κ B. Creation of DSB by IR is involved in over 77% of all modeled subprocesses. ATM and p53 subprocesses are involved in over 75% of all modeled subprocesses, while the NF-κ B module is involved in over 47% of them. The exact values for individual elementary processes are provided in Table 1, which contains only selected most relevant reactions. The results of the significance analysis for the remaining reactions (transitions) are given in Supplementary Table S3.
Structural analysis also includes searching for a set of transitions that is included in a sufficiently large number of t-invariant supports. To find that set, we used an algorithm described in the work by Gutowska and Formanowicz 36 . The result of this algorithm for the basic version of the model (541 t-invariants) allowed the determination of significant subprocess, which consists of the following transitions: t 1 , t 2 , t 13 , t 14 , t 16 , t 42 , t 43 , t 52 , t 70 , t 107 , t 147 (the names of the listed transitions are included in Supplementary Table S2). From biological perspective, such subset indicates key elements in the mechanism of the DNA damage response system based on ATM/p53/NF-κ B pathways. After the DNA damage, ATM detection module activates p53 and NF-κ B, which cooperate to determine the cell's response to the stimulus. Knockout of such subsets leads to exclusion of about 92% of all modeled subprocesses (what clearly emphasizes its importance). Table 1. Significance analysis for the proposed Petri net model with two excitations (541 t-inv).

Discussion
The independent analysis of the Petri net-based model, as well as the independent analysis of the ODE model, indicate significant components in the proposed model. In the case of the Petri net-based model analysis, we get information about crucial elementary processes (reactions), although we can also obtain such information for subprocesses and entire modules. In the case of ODE analysis, we obtain information about the parameters that have the greatest impact on the model response.
The significance analysis of the Petri net-based model is a structural analysis and allows indicating essential subprocesses in terms of information flow in such a net. While, the sensitivity analysis of the ODE-based model allows indicating the parameters whose change results in the greatest quantitative and qualitative changes in the model response. The results of the sensitivity analysis are influenced by both the values of parameters and inputs, therefore the analysis was performed for different simulation conditions. Some of the results for the significance and the sensitivity analysis are consistent (especially for the 4Gy dose; see Table 2). However, each method provides additional conclusions that another one is not able to yield. Therefore, their results indicate that the methods are complementary, in particular, in view of the following observations: • One of the most important elementary processes for the Petri net model concern ATM module and corresponding to them transitions t 42 , t 43 , t 52 , t 54 (see Supplementary Table S2)-from transcription, through translation, to ATM phosphorylation. This module is a key module for damage detection in response to IR radiation (it activates p53 and NF-κ B pathways). The parameters related to ATM module turned out to be important also in the ODE sensitivity analysis (parameters p 75 , p 76 , p 95 , 97 in the ODE model (see Supplementary Table S4)), but only for the simulation with 4Gy irradiation. Under extreme conditions (10Gy), the significance of these parameters is lower, which may be related to the existence of a strong, positive feedback loop in this regulatory system. • The sensitivity analysis showed a significant influence of parameters related to repair of double-stranded DNA damage (DSB repair), degradation of proteins, and degradation of ATM and Mdm2 transcripts. These elementary processes were not detected in the significance analysis for the Petri net model. www.nature.com/scientificreports/ The fact that degradation reactions have not been detected as essential may be associated with the structure of the Petri net model. In that framework, the larger number of subprocesses is dependent on a reaction, the more important this reaction is. A transition corresponding to degradation (called an output transition) does not affect other reactions (does not have places succeeding it) in any way other than limiting the amount of molecules being degraded (by taking tokens from a given place). The significance analysis revealed other important subprocesses, mainly related to transcription and translation.
In the case of DNA damage repair (parameter p 80 in the ODE model and corresponding transition t 0 in the Petri net model), the structure of the Petri net model also influences significance of this transition. If all repair modules were directly associated with this particular transition, its significance would be greater. In the analyzed model, damage repair is induced by p53pn, which is engaged in many other subprocesses. Therefore, though most transitions associated indirectly with DNA repair were found in the Petri net model to be important, the specific transition corresponding to damage repair by p53pn was not.
Interestingly, though in ODE model simulation for varying parameter p 75 (indicated as significant in Petri net model analysis, and less significant in ODE analysis), differences in the levels of the various forms of the ATM protein were noticed, no differences in settling time of the response were observed (see Fig. 5). Contrary to parameter p 75 , changes in parameter p 80 (indicated as not significant in Petri net model analysis, and significant in ODE analysis) do not result in differences in protein levels but result in different settling time of the response, see Fig. 6. Whether this is a coincidence, or a general rule, and how to explain these relationships remains an open question.

Conclusion
There are various methods of modeling complex biological systems, but each has its own limitations. Therefore, in this work, a dual approach, consisting in using two different modeling methods has been proposed. It facilitates a comprehensive structural and simulation analysis. Complementarity of sensitivity analysis (for ODEs model) and significance analysis (for Petri net model) has been confirmed by the results obtained for a sample largescale model of DNA damage response system based on the ATM/p53/NF-κ B pathways. Conclusions reached by application of only one of the methods would miss important information provided by another method, as described in sections "Results". An important, though unexpected finding is that changing parameter indicated as significant in the Petri net and less significant in ODE model leads to changing in protein levels, without changing of the settling time. On the other hand, changing a parameter rendered important in the ODE model, which is less important in the Petri net, leads to an opposite conclusion.
Our research shows that conducting a comprehensive analysis, both structural (for the Petri nets model) and parametric (for the ODE model), may not only allow for a better understanding of the studied phenomenon, but also allows to find essential parameters and significant modules which may prove important for the development of new, more effective therapies. Therefore, we show that the dual modeling approach can better supports the goals of medical therapy, and in the context of ATM/p53/NF-κ B pathways, it may lead to better outcomes in cancer therapy.

Methods
ODE-based approach. When the mathematical modeling of the intracellular signaling pathways is considered, usually one of the two main approaches is implemented: deterministic through ODEs or stochastic using Gillespie algorithm. Since parametric sensitivity analysis methods for the first type of models are well established, an ODE model is considered in this paper.
ODE models describe changes of concentration, or number of molecules, of molecular players, such as proteins, mRNAs, miRNAs, protein complexes etc., involved in the processes under investigation. They require parameter estimation which is based on the experimental results. Each process, taken into account, is represented by a single algebraic expression that appears in one or more equations and contains one or several parameters. A single parameter is associated with only one process. Therefore, high sensitivity of the model with respect to a given parameter indicates that system responses might be significantly altered by targeting the associated process. That, in turn, might be achieved by finding a molecule that would bind to a molecular player involved in this process, thus linking sensitivity analysis to search for molecular drug targets. Sensitivity analysis is the most important one taking into account the searching for the possible therapeutic target.
Sensitivity analysis became one of the necessary tools in in silico investigation of signaling pathways, as it provides information not only about dependence between parameter values and system behavior, but also about robustness of these systems 37 , a property that should characterize most of the pathways. This property is related to many biological phenomena such as homeostasis, stability, redundancy and plasticity 38 . SA methods can serve a number of useful purposes, e.g., uncover technical errors in the model, identify critical regions in the parameter space or establish priorities for research 39 . Among others, they also provide means to simplify high dimensional models that arise in systems biology 40 and can be used to indicate prospective molecular targets for new drugs against diseases associated with particular signaling pathways 41 .
Various methods of SA have been proposed, developed for either local or global sensitivity of systems under investigation. The local sensitivity analysis provides information about the effect of a small deviation of a single parameter from its nominal value on the system output. Global sensitivities, in turn, describe how the system output changes when multiple parameters change within a relatively wide range. Local methods are usually based on determining derivatives of model outputs with respect to parameters. Their undoubted advantage is simplicity and low computing power needed to perform the analysis, however, these methods do not guarantee correct results in the case of strongly non-linear models 39 . Global methods based on variance do not have such limitations, however, in the case of complex models of signaling pathways, they require enormous computing power, which may be a factor limiting their use. The compromise here may be variance-based methods from the OAT group, in which single parameters are randomly changed within a relatively wide range.
Petri net-based approach. Petri nets are mathematical objects with a structure of a directed weighted bipartite graph. A set of vertices of such a graph can be divided into two disjoint subsets in such a way that vertices belonging to the same subset are not connected by an arc. In the context of Petri nets, vertices being elements of one of these subsets are called places while elements of the other one are transitions. When a Petri net is a model of some system, places correspond to its elementary passive components (e.g., chemical compounds) and transitions are counterparts of elementary active components (e.g., chemical reactions). The bipartite graph describes the structure of the modeled system, but a Petri net is not just a graph. There is another type of component of such a net, i.e., tokens which bring a kind of dynamics into a net. They are located in places and flow from one place to another via transitions. This flow corresponds to a flow of substances, signals etc. through the modeled system. It is governed by the rule of firing transitions. According to this rule a transition is active if the number of tokens in all pre-places of this transition (i.e., the places which are its immediate predecessors) are equal to or greater than the weights of arcs connecting these places with the transition. An active transition can be fired what means that tokens can flow from its pre-places to its post-places (i.e., the places which are immediate successors of the transition). The number of flowing tokens is equal to the weights of the respective arcs.
Petri nets have a very intuitive graphical representation being helpful in understanding the structure and behavior of the modeled system. In this representation places are depicted as circles, transitions as rectangles, arcs www.nature.com/scientificreports/ as arrows and tokens as dots or numbers located in places. The weights, being positive integers, are represented as numbers associated with arcs (if a weight is equal to one, it is not shown).
Although the graphical representation of Petri nets is intuitive, it is not very well suited for a formal analysis of the proposed model. The used methods for analysis, are based on t-invariants, that is why for this purpose another representation, i.e., an incidence matrix, is used. In such matrix A = [a ij ] n×m , where n is the number of places and m is the number of transitions, rows correspond to places while columns correspond to transitions. Entry a ij , being an integer number, is equal to a difference between the numbers of tokens located in place p i after and before firing transition t j .
On the basis of matrix A invariants, which are especially important in the analysis of models of biological systems, can be calculated. There are two types of them, i.e., transition invariants (t-invariants) and place invariants (p-invariants). The former ones are vectors x ∈ N m being solutions of equation A · x = 0 , while the latter ones are vectors y ∈ N n satisfying equation y · A = 0.
In the analysis presented in this paper t-invariants are used. With every such an invariant a set of transitions, called its support, is associated. Support s(x) of t-invariant x contains those transitions which correspond to positive entries in vector x, i.e., s(x) = {t j : x j > 0, j = 1, 2, . . . , m} . If every transition t j from a support of t-invariant x is fired x j times, then the distribution of tokens over places of the whole net (also called a marking of the net) will not change, what means that the state of the modeled system is unchanged. Because a marking of a net corresponds to a state of the modeled system, it means that t-invariants correspond to subprocesses which do not change the state of the analyzed system and are crucial in the analysis of the Petri net-based models.
Methods of analysis of models based on Petri nets allow conducting structural analysis which includes significance analysis, knockout analysis, and searching for a significant subset of transitions. Such analyses are important to determine the significance of elementary processes or larger subprocesses.
The significance analysis is based on the frequency occurrence of each transition (elementary subprocess) in all supports of t-invariants 42 . On the basis of such frequency, the significance of transition expressed in percentage can be calculated, i.e., . Such an analysis allows to distinguish which subprocesses are more crucial for the functioning of the modeled system. Apart from determining the most significant elementary processes, it can be determined also the irrelevant ones. To complement the significance analysis a knockout analysis was performed (cf. 42 ).
The knockout analysis allows to estimate which subprocesses have been excluded in consequence of knockout of selected transition or sets of transitions. It allows to estimate how crucial selected elementary processes or whole subprocesses (sets of elementary processes) are, based on the number of excluded subprocesses. In such analysis, it is necessary to determine the number of t-invariants before and after knockout. Using these numbers of t-invariants, it is possible to estimate the percentage of excluded t-invariants as a result of a knockout and on this basis conclude about the significance of the knockouted elementary process or subprocess, i.e., (tinv b −tinv a )·100 An additional element of structural analysis is searching for a significant set of transitions. This approach is focused on finding a certain set of transitions that would be included in a sufficiently large number of t-invariant supports. We can distinguish two situations. The first one, a situation in which a certain set is included in almost all t-invariant supports (the greater number of t-invariant supports that include the searched set, the more important it is). Regardless of the expert knowledge about a given biological phenomenon, if a certain set occurs in more than 80% of all modeled subprocesses (supports of t-invariant), it will be considered as necessary for the system's functioning. A separate issue concerns the size of such a set. The smaller the set, the greater chances that it will be included in a bigger number of t-invariant supports, therefore, the larger sets of transitions (that can have assigned biological significance) are of interest. In the second situation, a set can occur in a smaller number of t-invariant supports, but if these t-invariants are crucial for the functioning of the model, then finding such set of transitions is also valuable. However, this case, unlike the first one, requires expert knowledge at every level of analysis. It should be realised which subprocesses are crucial for the functioning of the system in accordance with the state of biological knowledge. In both cases, a knockout of such a set of transitions results in turning off many or even the most of the modeled subprocesses (t-invariants). Transitions that are being elements of such a subset correspond to certain elementary processes. These elementary processes occur in many subprocesses, which may turn out to be essential to the functioning of a modeled system.