Petri net-based approach to modeling and analysis of selected aspects of the molecular regulation of angiogenesis.

The functioning of both normal and pathological tissues depends on an adequate supply of oxygen through the blood vessels. A process called angiogenesis, in which new endothelial cells and smooth muscles interact with each other, forming new blood vessels either from the existing ones or from a primary vascular plexus, is particularly important and interesting, due to new therapeutic possibilities it offers. This is a multi-step and very complex process, so an accurate understanding of the underlying mechanisms is a significant task, especially in recent years, with the constantly increasing amount of new data that must be taken into account. A systems approach is necessary for these studies because it is not sufficient to analyze the properties of the building blocks separately and an analysis of the whole network of interactions is essential. This approach is based on building a mathematical model of the system, while the model is expressed in the formal language of a mathematical theory. Recently, the theory of Petri nets was shown to be especially promising for the modeling and analysis of biological phenomena. This analysis, based mainly on t-invariants, has led to a particularly important finding that a direct link (close connection) exist between transforming growth factor β1 (TGF-β1), endothelial nitric oxide synthase (eNOS), nitric oxide (NO), and hypoxia-inducible factor 1, the molecules that play a crucial roles during angiogenesis. We have shown that TGF-β1 may participate in the inhibition of angiogenesis through the upregulation of eNOS expression, which is responsible for catalyzing NO production. The results obtained in the previous studies, concerning the effects of NO on angiogenesis, have not been conclusive, and therefore, our study may contribute to a better understanding of this phenomenon.

Introduction

The rapidly increasing amount of data concerning the molecular aspects of biological processes led to a conclusion that living organisms represent complex systems, composed of basic building blocks connected by a dense interaction network. The structure of this network determines the structure and functionality of the system. To understand the nature of living organisms, it is not sufficient to analyze the properties of their basic building blocks separately. The studies of the entire network of relations among them are necessary, and the methods suitable for the analysis of complex systems should be applied in these studies [1-3]. Many such methods have already been developed, especially in the area of technical sciences, and some of them can be used for the analyses of biological systems. However, in many cases, the existing methods of systems analysis should be adapted to biological systems, or new methods should be developed.

The basis of systems analysis is a formal model of the studied system, which can be expressed in a formal language of a mathematical theory. Differential equations are frequently used (and not only for the modeling of biological systems), since they represent very powerful mathematical tools. These equations require the precise values of the parameters that correspond to certain quantitative properties of the system, which are often difficult to determine in the case of biological systems, and this makes the construction of a differential equation-based model a difficult (or sometimes impossible) task.

Different approaches have been used for the modeling of biological systems, and especially promising one is based on Petri net theory, whose basic concepts were formulated in the early 1960s by Carl A. Petri, in the context of computer systems [4] (cf. [5-7]). For many years, these types of systems were the main area of application of Petri nets. However, in the mid-1990s, they have been applied for the modeling and analysis of biological phenomena as well, especially the metabolic networks, and this led to further investigations and the development of methods based on Petri nets for the modeling and analysis of complex biological systems [8, 9].

One of the advantages of Petri nets is their intuitive graphical representation, which is very helpful in the phases of the development of a biological system model, and the simulation of system behavior. Petri nets are mathematical objects, therefore they can be analyzed using strict mathematical methods, and many software tools for the simulation and analysis of this type of nets exist.

Building a precise model of a biological system represents a challenging task because such systems are probably among the most complex systems, and, contrary to, for example, the technical ones, they are not constructed by humans. To develop a formal model of a biological system, its structure has to be determined first, which may prove to be difficult.

In our study, the process of the regulation of angiogenesis has been modeled and analyzed. Angiogenesis is a complex physiological process, in which a coordinated sequence of gene expression, protein activity and distribution, as well as cell proliferation and migration, results in the development of new vessels from preexisting ones or from a primary vascular plexus. The development of new vasculature, particularly the formation of new capillaries from endothelial cells, is important in a number of pathological and homeostatic processes. In particular, angiogenesis is a well-described feature of the atherogenic process, in both coronary and carotid disease [10], and this process is necessary for the development and progression of cancer as well. Petri net model of angiogenesis has been already proposed in [11, 12], where the Stochastic Petri net and the reduction techniques have been discussed. These papers focus rather on the simplification of the models that describe complex biological systems and angiogenesis has been seen as a starting point for this concept. The studies presented in these articles relate mainly to VEGF and signaling pathways associated with it. Our research includes a larger number of processes but at a more general level.

Biological background of angiogenesis

To date, several factors that promote angiogenesis have been described, including physical factors, such as shear stress of blood flow or hypoxia, and biological factors, such as the secretion of various proangiogenic factors. However, hypoxic environment, often associated with tissue inflammation, is considered as the strongest activator of this process and it is among the most studied ones (reviewed in [13, 14]). The regulation of hypoxia-inducible factor 1 (HIF-1) has drawn a particular attention in the modeling of angiogenesis. The regulatory processes differ in the conditions of low oxygen status in cells (hypoxia) and the normal oxygen levels (normoxia).

Human cells require oxygen for the production of the adequate amounts of ATP necessary for metabolic activities, and they have developed the mechanisms of adaptation to hypoxic conditions or oxygen deprivation. Oxygen level changes can lead to the activation or repression of certain homeostatic regulatory genes, allowing the survival of tissues and cells despite the variable environmental conditions. Normal response to ischemia (restriction in blood supply to tissues) encompasses reparative mechanisms summarized by the term neovascularization, which includes 3 processes: angiogenesis, arteriogenesis, and vasculogenesis [15].

All the important compounds and reactions of the analyzed system will be described in the following paragraphs. In the model they will be represented by the so called places and transitions of a Petri net, while here a biological description is given.

The formation of new blood vessels can be divided into four stages, depending on the state of endothelial cells (ECs), which is the main indicator of neovascularization. In the first step, ECs are activated, usually by local hypoxic conditions. The strongest stimulator of hypoxia-driven angiogenesis is HIF-1, which plays an important role in the maintenance of cellular oxygen homeostasis, since it is responsible for the induction of expression of about dozen genes that encode proteins responsible for cell adaptation to hypoxia.

HIF-1 consists of two subunits, α and β, and HIF family consists of three α and three β isoforms [16], but we have considered the importance of HIF-1 α and β isoforms only, and the impact of HIF-2 α isoform on erythropoiesis in our model. The role of HIF-3 within the cell remains unclear, and this molecule was not included in the model.

The availability of HIF-1 is mainly determined by the α subunit, which is regulated in an oxygen-sensitive manner (controlled by ubiquitin-mediated degradation) under normoxic condition), in contrast to the β subunit, which is constitutively expressed [17]. The stability and activity of the α subunit is regulated by the post-translational modifications, such as hydroxylation, ubiquitination, acetylation, phosphorylation, and S-nitrosylation.

Prolyl hydroxylase enzymes (PHDs) exert tight control over HIF-1α degradation, and their activity depends on the concentration of oxygen. There are three PHD isoforms, PHD1, PHD2, and PHD3 [18], which have the potential to hydroxylate HIF-1α. PHD2 was shown to be the key limiting enzyme [19], and we have included only this isoform in our model. In normoxia, proline and asparagine residues of HIF-1α are hydroxylated by PHDs and by factor inhibiting HIF (FIH-1), respectively, in an oxygen and Fe 2+-dependent manner.

Both Fe(II) and cysteine residues of PHD2 catalytic domain can react with nitric oxide (NO). The most direct mechanism of the modulation of hydroxylase activity by NO is through competition with oxygen for the active-site Fe(II) [20]. NO can inhibit PHD and FIH activity by interacting with the enzyme-bound Fe(II), but this interaction between NO and PHD2 may be much more complex [20].

Hydroxylated HIF-1α proline residues serve as markers for the proteasomal degradation of HIF-1α, by binding to von Hippel-Lindau tumor suppressor, E3 ubiquitin protein ligase (VHL). VHL is transported from perinuclear granules to the nucleus or cytoplasm by a chaperone protein, VHL binding protein-1 (VBP-1) [21]. Following the VHL binding, HIF-1α is polyubiquitinated and rapidly degraded by the 26S-proteasomal system, keeping its expression levels low under normoxic conditions.

The modulation of HIF domains N-TAD and C-TAD is an additional control mechanism of its activity. These domains recruit transcriptional coactivators, such as CREB-binding protein (CBP)/p300 [19]. FIH-1 inhibits the interactions between HIF-1 and its coactivators through the hydroxylation of HIF-1α asparagine residues, which serve as a scaffold, linking various transcription factors to basal transcription machinery (BTM), making their role in the activation of HIF-1 especially important [22]. The subsequent binding of HIF-1 to p300/CBP facilitates the adaptation and survival of cells in an environment that changes from normoxia (21% O2) to hypoxia (1% O2) [19]. During hypoxic conditions, the lack of oxygen leads to strong endothelial activation and inhibits the activity of PHDs, resulting in the suppression of the degradation of α subunit, which promotes protein stabilization. Under these conditions, HIF-1α migrates from the cytoplasm to nucleus, binding to HIF-1β and undergoing dimerization, and forming a transcriptionally active HIF complex.

To date, approximately 100 genes involved in angiogenesis, metabolic adaptation, apoptosis, and metastasis have been identified as direct targets of HIF-1. Some of them, such as erythropoietin (EPO), vascular endothelial growth factor (VEGF), and glucose transporter type 1 (GLUT-1) have been included in our model.

EPO is regulated not only by HIF-1 but also by HIF-2α in sirtuin-1 (sirt-1)-dependent manner. The activation of sirt-1 may facilitate HIF-2α-directed production of EPO in cellular niches characterized by hypoxic stress, and possibly other environmental stresses [23].

VEGF is the most investigated angiogenic factor, whose activity leads to the production of matrix metalloproteinases (MMPs) and endothelium-derived NO. NO can stabilize and induce the transcriptional activity of HIF-1 through S-nitrosylation, where nitroso group is attached to the sulfur atom of cysteine [24, 25], preventing its degradation and increasing the interaction of HIF-1 complex with its co-activators, p300 and CBP.

NO synthesis in vascular ECs is catalyzed by endothelial nitric oxide synthase (eNOS), which can be activated through the hypoxia-stimulated erythropoiesis (through EPO and EPO receptor (EPOR)), the stimulation of tetrahydrobiopterin (BH4), essential cofactor of eNOS in ECs, and the cleavage from calmodulin (CaM) complex in caveolin-1 structures. eNOS is inactivated by CaM dissociation, allowing caveolin-1 to rebind, which is accompanied by renitrosylation of the enzyme and dephosphorylation of eNOS. Increased NO concentration is critical for the following stages of angiogenesis because it induces vessel permeability and the expression of glycolytic genes (GLUT-1) [26].

In the second stage of angiogenesis, MMPs degrade the basement membrane of vessels and extracellular matrix (ECM) components, which allows EC migration and proliferation. ECM represents a scaffold for growth factors that are crucial for the control of angiogenesis, by providing mechanical stability and allowing cell adhesion to ECM. It regulates the key stages of blood vessel morphogenesis and maturation [27].

MMPs may be involved in the recruitment of pericytes as well. The results obtained by Lethi et al. [28] suggest a dominant role of membrane-type 1 matrix metalloproteinase (MT1-MMP) in platelet-derived growth factor receptor-β (PDGFR-β) axis, leading to the recruitment and migration of mural cells (pericytes and vascular smooth muscles cells (VSMC)). Additionally, transforming growth factor β1 (TGF-β1) modulates MMPs. Activated TGF-β1 binds to the cell-surface type II receptor (TGF-βRII), inducing the activation of transmembrane serine/threonine kinases, activin receptor-like kinase-1 (ALK1) and activin receptor-like kinase-5 (ALK5). The activation of ECs is regulated in this way through two opposing pathways. ALK1 induces the signaling pathway that leads to an increase in EC proliferation and migration, while ALK5 inhibits these processes.

Following this, the third stage of angiogenesis is characterized by EC migration and proliferation, through the tip/stalk cell mechanism. During the formation of new vessels, tip cells form a leading head of newly proliferated cells, while stalk cells follow them. Recently, studies have shown that there is continuous competition between the cells of the new outgrowth to become tip cells, through a high expression of vascular endothelial growth factor receptor-2 (VEGFR-2), a major effector of VEGF signaling in angiogenesis, and low expression of vascular endothelial growth factor receptor-1 (VEGFR-1). Through the binding of vascular endothelial growth factor (VEGF) to VEGFR-2 tip cells are able to sense the VEGF gradient and follow it, and this binding promotes EC differentiation, proliferation, and sprouting.

In the last stage of angiogenesis, new vessels are stabilized through the recruitment of mural cells. The interaction between ECs and VSMC-like mural cells is essential for the formation of mature vascular structures. Growing vascular sprouts generate a concentration gradient of platelet-derived growth factor B (PDGFB). High levels of PDGFB in tip cells promote the recruitment of pericytes, which express PDGF receptor β (PDGFR-β). Biologically active form of PDGF is a dimer consisting of A and B chains. PDGF activity depends on different dimeric forms (AA, AB, or BB). PDGFR is a dimer as well, and it is formed by different combinations of α and β chains (α-α, α-β, β-β). PDGF-BB and basic fibroblast growth factor (bFGF) have stimulatory effects on ECs. bFGF/fibroblast growth factor receptor (FGFR) signaling pathway plays a critical role in the robust angiogenic response and the upregulation of PDGFR-α and PDGFR-β, which induce the angiogenic activity of PDGF-BB and PDGF-AB [29]. PDGFR-β activity may involve the cooperation with a G-protein coupled receptor for sphingosine-1-phosphate (EDG-1) that binds sphingosine-1-phosphate (SPP), a platelet-derived bioactive sphingolipid secreted by ECs. This induces the production of ECM proteins, promoting the migration of pericytes to the vessel outgrowth [30]. Furthermore, membrane type-1 matrix metalloproteinase (MT1-MMP), the prototypical member of MMP family subset, cooperates with SPP, stimulating EC migration and morphogenic differentiation into capillary-like structures.

Mural cell migration is facilitated by the binding of angiopoietin-1 (Ang-1) on mural cells to endothelial tyrosine kinase receptor 2 (Tie-2 receptor) on the EC surface. Ang-1 induces the expression of heparin-binding epidermal growth factor-like growth factor (HB-EGF) in ECs, and influences the potential of ECs to stimulate VSMC migration, suggesting an indirect mechanism by which Ang-1 recruits VSMCs [31]. Recent studies suggest that the stimulation of EPO leads to an increase in Ang-1 levels, indicating that EPO may regulate angiogenesis, at least partially, by modulating Ang-1 expression.

Angiopoietin-2 (Ang-2), which is produced and stored in Weibel-Palade bodies in ECs, generally functions as an Ang-1 antagonist, and it is upregulated in the hypoxic conditions, in a HIF-dependent manner [32]. Ang-2 is expressed only at the vascular remodeling sites, and plays a crucial role in the destabilization of vessels during normal or pathological angiogenesis [33].

Newly formed vessels are surrounded by the basement membrane, which consists of several types of laminins. It has been hypothesized that interaction between laminins and ECs stops the alterations in actin production, resulting in stationary morphology of ECs and terminating the angiogenic process (reviewed in [27]).

As a summary for the biological explanation of our angiogenesis model, Table 1 containing crucial processes and the literature data is given below.

Table 1

Main processes with corresponding literature references.

Process	References
HIF-1 regulatory pathway	[16, 17, 22, 24, 25]
HIF-2 regulatory pathway	[23]
The impact of the factors responsible for the adaptation of the cells to hypoxia	[23, 26, 28, 29, 31–33]

Because of the complex interactions between a large number of molecules, angiogenesis is a complex process, requiring a systems approach to completely understand and elucidate the underlying mechanisms.

In this study, a Petri net-based model of angiogenesis that relies on human cells, is presented and formally analyzed. This model is qualitative, describing the structure of the analyzed system. However, in the case of many biological systems, and angiogenesis as well, the structure is crucial for their functions, making it possible to draw some interesting biological conclusions based on the analysis of the presented model.

Methods

Petri nets

Petri nets are mathematical objects, with a structure of a directed bipartite graph G = (V, A). A set of vertices of a graph of this type can be divided into two disjoint subsets, V1 and V2, such that ending vertices of each arc (v, v) ∈ A of this graph belong to different subsets, i.e., v ∈ V1 ∧ v ∈ V2 or v ∈ V2 ∧ v ∈ V1. In Petri nets, vertices that are the elements of one of these subsets are called places, and the ones belonging to the other subset are called transitions. The arcs of this net are labelled by positive integers, called weights.

A directed bipartite graph corresponds to the structure of a Petri net, but this net is not static. One of the fundamental properties of Petri nets is dynamics, which is based on an existence of another net components called tokens. They reside in places and flow from one place to another through transitions (the direction of this flow is determined by the structure of the net, i.e., by the bipartite graph) (compare [5, 9, 34, 35]).

When a Petri net represents a model of a system, places correspond to its passive components, while transitions represent active components, and the flow of tokens models the flow of information through the system. A distribution of tokens over a set of places is called marking, and corresponds to a state of the modeled system.

Formally, a Petri net is a 5-tuple Q = (P, T, F, W, M0), where:

P = {p1, p2,…, p} is a finite set of places,

T = {t1, t2,…, t} is a finite set of transitions,

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

is a weight function,

is an initial marking,

P ∩ T = ∅ ∧ P ∪ T ≠ ∅ [5].

Every transition can have a set of pre-places, i.e., the ones who are its immediate predecessors. Analogously, a transition also can have a set of post-places as its immediate successors. Similarly, every place can have sets of pre- and post-transitions. The flow of tokens through the net is governed by the transition firing rule, where a transition is enabled if the number of tokens residing in every of its pre-places is at least equal to the weight of the arc connecting this place with the transition. An enabled transition can fire, meaning that the tokens can flow from their pre-places to their post-places. The number of tokens flowing between the place and the transition is equal to the weight of the arc connecting these two vertices.

Petri nets are mathematical objects that have very intuitive graphical representation, which facilitates the understanding of the structure and behavior of the net. In this representation, the transitions are shown as rectangles or bars, places as circles, tokens as dots or numbers residing in the places, and arcs as arrows. The arcs are labeled with weights, except when the weight equals one (i.e., usually, the weights equal to one are not explicitly shown in the graphical representation).

This graphical representation of Petri net, although intuitive and useful, is not appropriate for the analysis of net properties. A different representation, called incidence matrix, is used for this purpose. In this matrix, A = [a] rows correspond to places and columns correspond to transitions. The matrix entries are integers, and entry a, i = 1, 2,…,n, j = 1, 2,…,m is equal to the difference between the number of tokens residing in place p before and after firing transition t.

t-invariants

If a Petri net is a model of a metabolic network, invariants of the net and an analysis of the relations among them are of particular importance [6, 9, 36–39]. Two kinds of invariants exist here, t-invariants and p-invariants. The former is vector , which is a solution of the equation

Similarly, a p-invariant is vector , representing a solution to the equation

Each of these invariants corresponds to a set of transitions or places. More precisely, a support supp(x) is associated with each t-invariant x, representing a set of transitions that correspond to positive entries in x, i.e., supp(x) = {t : x > 0, j = 1, 2,…,m}. Every p-invariant y has an associated support supp(y) = {p : y > 0, i = 1, 2,…,n}.

Especially important are the minimal invariants, i.e., those whose supports do not contain a support of any other invariant. The consideration of only minimal invariants is sufficient because every invariant can be obtained as a linear combination of the minimal ones. A support of a t-invariant contains transitions whose firing a proper number of times (i.e., the number of times being the invariant entry which corresponds to a given transition) reproduces the marking of a net (i.e., the state of the system modeled by the net does not change). A support of a p-invariant contains places whose weighted sum of tokens is constant (the weights are the positive entries of the invariant).

With the t-invariants a question of their realisability arises [6]. In short, the analysis of t-invariants should base on the feasible ones. The existence of non-feasible t-invariants depends on two factors, i.e., insufficient number of tokens in each of the minimal p-invariants in the initial marking and/or read arcs (i.e., arcs going in two directions between a place and a transition). From the further analysis of our model properties it will be clear that such factors are not an issue for the Petri net-based model presented in this paper, therefore all generated t-invariants can be considered as feasible.

Each transition in a Petri net should belong to a support of a feasible t-invariant, i.e., the net should be covered by t-invariants. In that case every modeled elementary biological process contributes to the behavior of the biological system. Every t-invariant corresponds to a subprocess and every transition corresponds to an elementary subprocess. Since every t-invariant has an associated subset of transitions (a support), and the supports can have non-empty intersections, the subprocesses, corresponding to various t-invariants, may be composed of common elementary subprocesses, i.e., modeled by transitions that belong to an intersection of the corresponding supports. Therefore, the analysis of the relations of this type among t-invariants may lead to the elucidation of interactions and dependencies between the subprocesses of the modeled biological system.

t-clusters

To find relations among t-invariants, they are grouped into sets called t-clusters [40]. The goal is to acquire different clusters in such a way that they will group together t-invariants modeling subprocesses which possibly interact with each other. On the other hand, the differences between the clusters should be significant enough to distinguish and name such different sets of processes that have been previously identified as t-invariants in the modeled biological system. This is neither a simple nor an automatic process. There are three main “parameters” in this task and changing any one of them will result in a different clustering, i.e., a set of clusters for the further biological analysis. These are the distance metric, the linkage (joining) algorithm and the desired number of clusters. In the presented work we have decided to use hierarchical clustering and we have tested different combinations of the above parameters using scripts written in R language, available as S1 Script. Libraries used in our scripts provide different distance metrics and linkage algorithms. The distance metrics are referred to as Binary, Canberra, Pearson’s Correlation, Euclidean, Manhattan, Maximum, Minkowski, and Uncentered Pearson. The clustering algorithms are as follows: Unweighted Pair Group Method with Arithmetic mean (UPGMA), McQuitty linkage, Median linkage (also known as Weighted Pair-Group Method using Centroids, WPGMC), Single linkage, Complete linkage, Centroid linkage (also known as Unweighted Pair-Group Method using Centroids, UPGMC), and Ward’s linkage method [41]. As it has been already stated, the right combination of the distance measure and the clustering method is not an obvious choice. Moreover, the number of clusters is also an important, yet unknown variable. To evaluate the results, the Mean Split Silhouette index (MSS) [42, 43] has been used. Scripts have been prepared in such a way that for each selected distance measure and linkage algorithm, the different clusterings are generated for the selected range of the numbers of clusters, and their content is evaluated using MSS index. The results are presented as PDF files, also for each clustering a dendrogram is generated. Choosing the proper clustering is a challenging task, because the MSS evaluation can help, but alone is not enough. Generally, the most interesting results are obtained by the clustering methods that give the smallest number of single-invariant clusters and the highest MSS value. As it will be discussed later, an initial biological examination of promising clusterings is required, because the evaluation done on the basis of MSS value only may be not sufficient to draw conclusions interesting from the biological point of view.

MCT and ADT sets

A support of a t-invariant defines a minimal self-contained subnet, representing a single subprocess within the modeled biological system. All the subprocesses defined by the t-invariants can be analyzed easier if some common subsets of transitions in their supports could be identified. Such subsets are called Maximal Common Transition sets (MCT sets for short), providing an aid in the analysis of biological systems modeled by Petri nets [6, 37]. A set of this type contains transitions which are elements of supports of exactly the same t-invariants. They can be computed using, e.g., a simple algorithm described in [37].

In contrast to t-clusters, which usually induce overlapping subnetworks, MCT sets partition the net into disjunctive subnets and can be considered as the smallest biologically meaningful units [6]. They can have various sizes and also a single transition can form a (trivial) MCT set on its own. For the biological examination of such functional units we consider only non-trivial MCT sets. It should be noted that an MCT set can form a disconnected subnet, i.e., within such a set there can be at least one pair of transitions such that in every path connecting them there is at least one transition not belonging to this set. In [44, 45] different name for such sets has been proposed—maximal Abstract Dependent Transition (ADT) sets. It should be stressed that MCT sets are exactly the same as maximal ADT sets (cf. [6] and [44]). For the latter, a further decomposition has been proposed in [45]—into ADT sets which are not maximal, but induce connected subnets. One of the possible reasons for such a decomposition is to create a hierarchical structure of a net, where connected subnets induced by ADT sets are hidden in the so called macro transitions. Such a hierarchical structure may allow identification of important net components that connect various ADT sets, e.g., interface places for which additional studies can be performed [45]. In our approach we do not use hierarchical net structure, therefore we will use MCT sets. In the presented Petri net based model there is only one MCT set which induces a disconnected subnet and it consists of only two transitions. We have decided its further decomposition would not contribute to the presented analysis.

At the end a difference between t-clusters and MCT sets should be discussed. Both analytical approaches are based on t-invariants. MCT sets consist of transitions, and as it has been already mentioned, they represent a subnet of a net structure, not necessarily a connected one. t-clusters on the other hand group similar t-invariants. The supports of the t-invariants consists of transitions, therefore it is theoretically possible to represent a t-cluster as a multiset of transitions. The same transition appears in such a multiset as many times as it is present in all the supports of t-invariants of a given t-cluster. Such a representation can provide addition help in the understanding of a function of the cluster, because one can highlight a whole region of the net which is involved in all subprocesses which t-invariants of the t-cluster represent. This is just another way of looking on the t-clusters, but it should be noted that such a representation loses the information about the distinct subprocesses represented by the cluster t-invariants. Therefore, assigning a biological function to the t-cluster should first be performed by the analysis of group of its t-invariants. Analysis of the subnet representing the t-cluster can be auxiliary in such an analysis.

Knockout analysis

Further analysis has been performed using knockout analysis, which has been previously presented on the model of Duchenne muscular dystrophy [40]. The knockout analysis allows investigating which subprocesses will be influenced (i.e., knocked-out) if some elementary functional unit is disabled. For this task MonaLisa software [46] has been used. Each activity is represented by an MCT set, either a multi-element or a trivial one. Since transitions that belong to the same MCT set always appear together in the same t-invariants, the knockout impact of the members of a given MCT set is always the same. For every activity knocked out, the percent of affected t-invariants have been determined, thus leading to a ranking of the functional units within the model that carry the most and the least impact for the biological system processes.

Such a knockout analysis is based on the structural properties of the model, in this case on the t-invariants. It only gives the information about the percent of all processes that are dependent on the chosen (and disabled) functional unit. In fact, disabling functionality of any given transition should influence the rest of the net in some way that cannot be measured using the described method. Therefore, we have also studied the knockout behavior on the basis of the net dynamic, i.e., by observing tokens fluctuations and changing of transition firing frequencies when some functional units are being disabled. Data describing the normal model behavior (when nothing is disabled) are gathered, they can be compared with the data sets obtained, when a single transition or a non-trivial MCT set has been disabled.

Such a different type of knockout analysis involves a simulation approach. It is performed using token game simulation, in a mode similar to “synchronous” as described in [47], yet in our approach every active transition has 50% chance of firing. Each simulation involves 10,000 steps from the same initial marking and the simulations are repeated 4,000 times. For every step and every repetition data concerning transition firing and tokens distribution are gathered. Finally, a data package is created, containing average firing chance for each transition, as well as the information about the average tokens sum in every place. A standard deviation is calculated for transition firings and for the tokens accumulation from such 4,000 separate simulations to evaluate how repetitive is the behavior of the model. First data package gathered in that way contains the information about the model behavior when nothing is disabled. Then, for every transition a similar simulation procedure is repeated, but this time involving a knockout of that given transition. It means that a selected transition is always considered to be inactive by the simulation engine in every step in every repetition, no matter of the number of tokens in its pre-places. Then the change of firing chances of active transitions, when some important selected transitions, have been marked as knocked out, can be compared.

The net presented in this paper has been created using Snoopy software [48], while the analysis of its formal properties has been performed with INA [49]. Export to INA Petri net format as well as into many other PN file formats is possible and quite easy using Snoopy. Knockout analysis based on t-invariants has been performed using MonaLisa [46]. We have used our own, not yet published tool to acquire data for the described knockout analysis based on token game simulation. Computations involving cluster analysis have been performed using scripts written in R language. They required two additional libraries in R: cluster (v2.0.3) and amap (Another Multidimensional Analysis Package, v0.8.14). Scripts and the t-invariants input file (in fact a simple CSV format) are given in S1 Script.

Results and discussion

Analysis of the proposed model

The Petri net-based model of the angiogenic process developed and analyzed in this study (shown in Fig 1, also given as S1 File) contains 74 transitions and 63 places. Some places are identically named, but are shown as two concentric circles to improve the model readability, as these are the logical places, representing independent graphical elements that correspond to the same vertex of the net. The names of places and transitions are listed in Tables 2 and 3, respectively. Places correspond to the passive elements of the modeled system, i.e., chemical compounds, reaction products, and substrates, while transitions correspond to the elementary processes occurring in the system.

Fig 1

The Petri net model of angiogenesis with marked non-trivial MCT-sets.

Table 2

The list of places.

Place	Corresponding molecules	Place	Corresponding molecules
p0	HIF-1α-SH stable form not degraded in proteasome	p32	Tie-2 receptor for Ang-1 and Ang-2
p1	Active transcriptional factors	p33	Angiopoietin and Tie receptor complexes
p2	VEGF	p34	SPP
p3	EPO	p35	EDG-1
p4	VEGFR-2 on the endothelial surface	p36	SPP and EDG-1 complex
p5	VEGF and VEGFR-2 complex	p37	HB-EGF
p6	GLUT-1	p38	EGFR
p7	PHD-2, activated	p39	HB-EGF and EGFR complex
p8	FIH-1, activated	p40	TGF-1β
p9	HIF-1 β, constitutively expressed	p41	TGFβRII
p10	HIF-1α-OH, inactivated	p42	TGF-1β-TGFβRII complex
p11	HIF-1α, degraded	p43	Oxygen
p12	VHL and VBP-1 complex	p44	CaM and Ca2+ complex
p13	CBP and p300	p45	PDGFRs in the absence of ligands
p14	HIF-1 hydroxylases, inactivated	p46	MMPs
p15	eNOS, activated	p47	Hypoxia
p16	NO	p48	EPOR
p17	bFGF	p49	EPO and EPOR complex
p18	FGFRs, activated	p50	HIF-1α and HIF-1 β complex, complete (HIF molecule)
p19	PDGFRs, upregulated	p51	Tyrosine kinase signaling cascade
p20	PDGF-AB and PDGF-BB, secreted	p52	CaM and EPOR complex
p21	PDGF-AB and PDGF-BB, and PDGFR complex	p53	CaM free of Ca2+
p22	Endothelium activated by hypoxia	p54	Increased intracellular Ca2+ level
p23	Recruited pericytes	p55	HIF-2α
p24	ECs, upon proliferation and migration	p56	Sirt-1
p25	Pericytes and smooth muscle cells, upon proliferation and migration	p57	Fe(II)
p26	Stable vessels	p58	HIF-1α, hydroxylated (HIF-1α-OH-OH)
p27	CaM-Ca2+-caveolin 1 complex	p59	HIF-1α-OH-OH-VHL-VBP1 complex
p28	Caveolin-1	p60	Ubiqiutin
p29	Increased shear stress	p61	BTM
p30	Ang-1	p62	Normoxia
p31	Ang-2

Table 3

The list of transitions.

Transition	Corresponding reactions	Transition	Corresponding reactions
t0	HIF-1α-SH and HIF-1β dimerization	t37	Constitutive expression of Ang-1
t1	The genes responsible for cell adaptation to hypoxia	t38	Tie expression
t2	VEGF and VEGFR-2 binding	t39	The expression of angiopoietins
t3	VEGFR-2 expression	t40	The remodeling of vessels
t4	Tyrosine kinase signaling cascade initiation	t41	The binding of SPP and EDG-1
t5	GLUT-1 induced processes	t42	SPP synthesis
t6	The binding of EPO and EPOR	t43	EDG-1 synthesis
t7	Constitutive expression of HIF-1β	t44	The enhancement and increase of synthesis of ECM, and the release of HB-EGF
t8	Proline hydroxylation	t45	The binding of HB-EGF and EGFR
t9	Conservative hydroxylation of asparagine	t46	The binding of TGF-β1 and TGFβRII
t10	The binding of HIF-1-OH-OH and VHL-VBP-1	t47	ALK-1 or ALK-5 activation
t11	The lack of HIF-1α and HIF-1β dimerization	t48	The expression of TGFβRII
t12	The processes leading to low oxygen concentrations in an organism	t49	TGF-β1 synthesis
t13	HIF-1 hydroxylase inactivation	t50	Complex binding to HREs hypoxia responsive elements
t14	The recruitment of transcriptional coactivators, such as CBP/p300	t51	The synthesis of factors that stimulate vessel permeability in various ways
t15	The lack of HIF-1α degradation	t52	The synthesis of factors that stimulate proliferation and survival by bFGF
t16	The normalization of oxygen status in the organism	t53	The synthesis of factors that variously stimulate migration through MMPs
t17	VHL synthesis	t54	The increase in intracellular Ca2+
t18	NO synthesis	t55	CaM and EPOR binding
t19	The recruitment and invasion of pericytes by MMPs	t56	The upregulation of EPOR signaling pathways
t20	The secretion of PDGF-AB and PDGF-BB by pericytes	t57	The binding of CaM and Ca2+
t21	FGFRs activation	t58	CaM synthesis
t22	PDGFR upregulation	t59	HIF-2α influence on EPO
t23	The binding of PDGF-AB and PDGF-BB with the activated PDGFRs	t60	The influence of sirt-1 on HIF-2α
t24	The removal of ligands	t61	Sirt-1 synthesis
t25	Pericyte and smooth muscle cell migration	t62	The increase of Fe(II)
t26	Strong endothelial activation	t63	HIF-1α asparagine hydroxylase activation
t27	Endothelial cells proliferation and migration	t64	HIF-1α proline hydroxylase activation
t28	The binding process	t65	Oxygenation
t29	The use of new stable vessels	t66	The activation of proteasome degradation
t30	The activation of eNOS through the interactions with CaM and Ca2+ complex	t67	Ubiquitin synthesis
t31	The process of cleavage	t68	The regulation of gene expression through BTM
t32	Caveolin usage	t69	S-nitrosylation
t33	The processes that increase shear stress	t70	PHD activity influenced by NO
t34	CaM and caveolin-1 binding	t71	The processes that lead to the increase in NO
t35	Caveolin-1 synthesis in ER	t72	EPOR synthesis induced by hypoxia
t36	The binding of angiopoietins with Tie receptor on ECs	t73	Normal cellular state

The presented model is a discrete Petri net with the weights of all arcs being equal to one. The model contains neither information about the speed of reactions, nor substrate or product quantities, as this is a qualitative model, containing only information about the structure of the analyzed system. Functions of a biological system generally depend on its structure, and the analysis of this structure may lead to the deeper understanding of the system properties.

The net is ordinary, meaning that all weights of arcs are equal to one. For each place, all outgoing arcs have the same weight (no quantitative data is considered), so the model is homogeneous. The net is connected but not in the strong sense, i.e., there are undirected paths between any two places, but there may be no directed path between all pairs of places. Connectivity means that there are no independent processes within the model. The net is not structurally conflict-free because it contains places with two or more outgoing arcs (for instance transitions t7, t63 and t64 have the same preceding place p62). Additionally, the model is pure, i.e., it does not contain any pairs of oppositely directed arcs (read arcs). There are no p-invariants, while there are 48 minimal t-invariants, which cover the net.

When modeling a biological system, the most important part of the analysis is based on t-invariants, and the full list of these invariants is given in Table 4. MCT sets have been determined based on t-invariants and they are listed in Table 5, together with their biological interpretations. It should be noted that only MCT 9 corresponds to a disconnected subnet. It is however one of the smallest sets (2-element one), therefore it has been decided that its further decomposition into smaller sets (each one in that case would be a trivial, i.e., 1-transition MCT set) would not contribute in a significant way to the model analysis.

Table 4

The list of t-invariants.

The second and third columns show the elements of the support of t-invariants listed in the first column.

t-invariant	MCT sets	Single transitions
x1	m10
x2	m7	t49
x3	m9, m11	t40
x4	m6, m11	t59
x5	m3	t64, t65, t73
x6	m2, m8	t12, t18, t30, t51, t69, t73
x7	m3	t12, t16, t64, t65
x8	m2, m8	t12, t16, t18, t30, t51, t65, t69
x9	m3	t65, t70, t71, t73
x10	m2, m6	t12, t18, t51, t69, t71, t73
x11	m1, m2, m6, m7, m9, m11	t12, t27, t69, t71, t73
x12	m1, m2, m4, m6, m7, m8, m9, m11	t12, t69, t71, t73
x13	m2, m5, m6	t12, t18, t51, t71, t73
x14	m1, m2, m5, m6, m7, m9, m11	t12, t27, t71, t73
x15	m1, m2, m4, m5, m6, m7, m8, m9, m11	t12, t71, t73
x16	m1, m2, m7, m8, m9, m11	t12, t18, t27, t30, t69, t71, t73
x17	m1, m2, m4, m7, m8, m9, m11	t12, t18, t30, t69, t71, t73
x18	m1, m2, m5, m7, m8, m9, m11	t12, t18, t27, t30, t71, t73
x19	m1, m2, m4, m5, m7, m8, m9, m11	t12, t18, t30, t71, t73
x20	m3	t12, t16, t65, t70, t71
x21	m2, m6	t12, t16, t18, t51, t65, t69, t71
x22	m1, m2, m6, m7, m9, m11	t12, t16, t27, t65, t69, t71
x23	m1, m2, m4, m6, m7, m8, m9, m11	t12, t16, t65, t69, t71
x24	m2, m5, m6	t12, t16, t18, t51, t65, t71
x25	m1, m2, m5, m6, m7, m9, m11	t12, t16, t27, t65, t71
x26	m1, m2, m4, m5, m6, m7, m8, m9, m11	t12, t16, t65, t71
x27	m1, m2, m7, m8, m9, m11	t12, t16, t18, t27, t30, t65, t69, t71
x28	m1, m2, m4, m7, m8, m9, m11	t12, t16, t18, t30, t65, t69, t71
x29	m1, m2, m5, m7, m8, m9, m11	t12, t16, t18, t27, t30, t65, t71
x30	m1, m2, m4, m5, m7, m8, m9, m11	t12, t16, t18, t30, t65, t71
x31	m3, m8, m11	t18, t30, t59, t65, t70, t73
x32	m1, m2, m7, m8, m9, m11	t12, t18, t27, t30, t59, t69, t73
x33	m1, m2, m4, m7, m8, m9, m11	t12, t18, t30, t59, t69, t73
x34	m1, m2, m5, m7, m8, m9, m11	t12, t18, t27, t30, t59, t73
x35	m1, m2, m4, m5, m7, m8, m9, m11	t12, t18, t30, t59, t73
x36	m3, m8, m11	t12, t16, t18, t30, t59, t65, t70
x37	m1, m2, m7, m8, m9, m11	t12, t16, t18, t27, t30, t59, t65, t69
x38	m1, m2, m4, m7, m8, m9, m11	t12, t16, t18, t30, t59, t65, t69
x39	m1, m2, m5, m7, m8, m9, m11	t12, t16, t18, t27, t30, t59, t65
x40	m1, m2, m4, m5, m7, m8, m9, m11	t12, t16, t18, t30, t59, t65
x41	m2, m3, m5, m8	t12, t18, t30, t51, t65, t70, t73
x42	m2, m5, m6, m8	t12, t18, t30, t51, t73
x43	m1, m2, m5, m7, m8, m9, m11	t12, t18, t27, t30, t51, t73
x44	m1, m2, m4, m5, m7, m8, m9, m11	t12, t18, t30, t51, t73
x45	m2, m3, m5, m8	t12, t16, t18, t30, t51, t65, t70
x46	m2, m5, m6, m8	t12, t16, t18, t30, t51, t65
x47	m1, m2, m5, m7, m8, m9, m11	t12, t16, t18, t27, t30, t51, t65
x48	m1, m2, m4, m5, m7, m8, m9, m11	t12, t16, t18, t30, t51, t65

Table 5

The list of non-trivial MCT sets.

MCT set	Contained transitions	Biological interpretation
m1	t19, t21, t22, t23, t24, t25, t28, t36, t37, t38, t41, t42, t43, t44, t45, t52, t53	The stimulation of neovascularization caused by the selected factors with mitogenic activity in endothelial cells
m2	t0, t1, t2, t3, t4, t5, t7, t14, t20, t50, t68	The activation of the signaling pathways involving mitogen-activated kinases
m3	t8, t9, t10, t11, t17, t62, t63, t66, t67	The degradation of HIF-1α under normoxic conditions
m4	t6, t55, t56, t72	EPO pathway regulation
m5	t13, t15, t26	HIF-1α accumulation under hypoxic conditions
m6	t31, t33, t34	Shear stress-mediated CaM-caveolin-1 cleavage
m7	t46, t47, t48	The activation of TGF-β signaling pathway
m8	t54, t57, t58	The activation of eNOS modulated by CaM
m9	t29, t39	The stabilization of new vessels by angiopoietins
m10	t32, t35	The regulation of caveolin-1 expression in caveolae
m11	t60, t61	Sirt-1 dependent HIF-2 regulation

t-invariants, MCT sets, and clusters-based analyses

The analysis of the model is primarily based on t-invariants. On the basis of these transitions have been grouped into 11 MCT sets, while t-invariants has been grouped into t-clusters. The biological significance of these objects has been determined. t-clusters contain t-invariants corresponding to biological subprocesses that interact with each other in some way and new understanding of biological processes can be obtained by a proper construction of t-clusters and the analysis of their contents. This may lead to the discovery of dependencies among the biological subprocesses. For simplicity, in further text, we will sometimes write that a t-invariant contains some MCT sets or transitions to say that these MCT sets and transitions are elements of a support of this invariant.

One part of a table that contains the evaluation results of the clustering methods is presented in Fig 2. The complete data are given as an Excel file in S1 Table. Each result is described by two values; the first, narrow column represents the number of single-invariant clusters, while the second column contains the MSS index value for the entire clustering (obtained using a given method for a given number of clusters). The pairs of columns (except the first one) in the Table correspond to the clustering algorithms (linkage methods), while the rows contain the results for given numbers of clusters, for each distance metric.

Fig 2

A part of the MSS evaluation table of the Pearson’s Correlation distance and different joining algorithms.

For the biological analysis 13 clusters has been chosen, obtained using Correlated Person metric and UMPGA (average) linkage method. It requires a further explanation and points out the main issue with the cluster analysis in general. One can see that, for example, for the aforementioned clustering the MSS value is 0.5088, while e.g., for 11 clusters obtained using Binary metric and UPGMA linkage algorithm, MSS value of the entire clustering is higher, i.e., 0.5856 (see S1 Table). However, this clustering contains only seven clusters composed of more than one t-invariant. A detailed analysis of them, shown in Table 6, demonstrates that cluster number 8 contains 28 t-invariants, more than one-half of all t-invariants. This type of cluster corresponds to most of the transitions of the analyzed Petri net and it is very difficult to analyze. In other words, it is quite hard to give such a subset of processes a relevant, common biological function. Therefore, the results presented in Table 6 represent a first step when choosing of the best clustering. The second stage of the analysis, taking into account the composition of each clustering with one of the highest global MSS value is necessary. A careful choosing of distance metric has been also very important step. It seems that correlated Person metric is quite good in measuring the distances for the invariants (i.e., properly distinguishing their similarity) for the clusters computation.

Table 6

Two clusterings obtained with different distance metrics and joining algorithms.

Mean Split Silhouette (MSS) values are given for every cluster.

Correlation-UPGMA			Binary-UPGMA
Cluster number	Number of t-invariants	MSS value	Cluster number	Number of t-invariants	MSS value
1	1	0.00	1	1	0.00
2	1	0.00	2	1	0.00
3	1	0.00	3	1	0.00
4	1	0.00	4	1	0.00
5	2	0.38	5	4	0.46
6	10	0.51	6	2	0.56
7	2	0.49	7	2	0.72
8	4	0.64	8	28	0.69
9	8	0.29	9	4	0.33
10	4	0.66	10	2	0.85
11	4	0.59	11	2	0.63
12	2	0.82
13	8	0.75

Table 7 shows the chosen clusters and their biological interpretation, while Fig 3 represents the dendrogram of the chosen clustering. Red horizontal line represents the threshold level that distinguishes the 13 clusters. The number of t-invariants for each cluster is given below the red line. The obtained clustering served as the basis for the investigations of the relationships between biological subprocesses that correspond to t-invariants.

Table 7

The list of clusters.

Cluster	Contained invariants	Biological significance
c1	x1	The regulation of caveolin-1 expression caveolae
c2	x2	The activation of TGF-β signaling pathway
c3	x3	The remodeling of vessels and their stabilization, influenced by angiopoietins and HIF-2
c4	x4	The regulation of erythropoiesis via HIF-2 and caveolin-1
c5	x5, x9	The degradation of HIF-1α under normoxic conditions, influenced by processes that stimulate NO signaling pathway
c6	x6, x8, x41, x42, x43, x44, x45, x46, x47, x48	The relationships between metabolic degradation-related pathways and the formation of HIF-1 particles with no effect on remodeling and the regulation of expression of caveolin-1 in caveolae
c7	x7, x20	The degradation of HIF-1α under normoxic conditions, influenced by processes that stimulate NO signaling pathway and processes that lead to transition from hypoxia to normoxia, terminating angiogenesis
c8	x10, x13, x21, x24	The impact of oxidative stress on the regulation of HIF-1α and HIF-1β dimerization, together with the formation of active transcriptional factors, responsible for cell adaptation to hypoxia, with particular emphasis on VEGF signaling pathway
c9	x11, x12, x16, x17, x22, x23, x27, x28	The effect of S-nitrosylation on the dimerization of HIF-1α molecules, together with the formation of active transcriptional factors responsible for cell adaptation to hypoxia, excluding the influence of the endothelium
c10	x14, x15, x25, x26	The effect of S-nitrosylation on the dimerization of HIF-1α molecules, together with the formation of active transcriptional factors responsible for cell adaptation to hypoxia
c11	x18, x19, x29, x30	The effect of S-nitrosylation on the dimerization of HIF-1α molecules, with particular emphasis on TGFβ signaling pathway, when switching from normoxia to hypoxia
c12	x31, x36	The degradation of HIF-1α under normoxic conditions, influenced by CaM-Ca and HIF-2α-sirt-1 pathways
c13	x32, x33, x34, x35, x37, x38, x39, x40	The dimerization of HIF-1α molecules, together with the formation of active transcriptional factors responsible for cell adaptation to hypoxia, excluding the binding between CaM and ceveolin-1

Fig 3

Cluster dendrogram.

The analysis of t-invariants showed that the most important process condition affecting almost all of the subprocesses in the presented model is hypoxia and processes that contribute to it are involved in the response to hypoxia. Transition t12, corresponding to hypoxia, was present in 42 of 48 t-invariants. Six other t-invariants were related either to the regulatory processes, such as the: regulation of eNOS via caveolae (x1), regulation of the state of endothelium via TGF-β1 (x2), regulation of vascular remodeling by Ang-2 (x3), regulation of the vascular shear stress response via caveolin-1 and sirt-1 (x4), or to processes that occur under normoxic conditions, such as the degradation of HIF-1α by the ubiquitin-proteasome system (x5, and the regulation of HIF-1α via by NO (x31)). The last t-invariant, i.e., x31, corresponds to subprocess comprising phenomena such as proteosomal degradation of HIF-1α in normoxia, EPO-mediated activation of eNOS, the increase in cellular calcium levels influenced by EPO, and the pathway sirt-1/HIF-2α/EPO pathway. These phenomena are found in one t-invariant; hence, they are closely related to each other. What is especially interesting is that they reflect the crosstalk between NO and HIF signaling pathways. These results are reflected in the content of cluster c12, which is formed from two t-invariants, i.e., x31 and x36, which links HIF and sirt-1 together. By demonstrating that sirt-1 can regulate the activity of HIF proteins, this cluster links HIF and sirt-1 together.

The results of the studies that investigated whether there is a direct interaction between HIF-1α and sirt-1 have been conflicting. Dioum et al. [23] have reported that sirt-1 does not target HIF-1α, but it deacetylates HIF-2α and their interaction promotes HIF-2 transcriptional activity. Different studies suggest that sirt-1 positively regulates the transcriptional activity of HIF-1 [50]. Our results indicate that there is a relationship between HIF-1α and sirt-1, and HIF-2α and sirt-1.

Forty-two of the remaining t-invariants include subprocesses in which one transition always corresponds to hypoxia. Analysis of these t-invariants showed that none of the supports of t-invariants contain a pair of MCT sets m3 (degradation of HIF-1α under normoxic conditions) and m10 (regulation of the expression of caveolin-1 in caveolae). This is consistent with the results obtained in the previous studies. Recent studies have highlighted the significance of caveolin-1-dependent control of eNOS, stressing the relationship between these two proteins [51]. Under the baseline conditions, caveolin-1 binds eNOS and inhibits the production of NO through the caveolin-1 scaffolding domain. This implies that the association of eNOS with caveolin-1 maintains eNOS in an inactive state, occupying the calcium/calmodulin binding site of eNOS. Our model showed that, under the normoxic conditions, chronic exposure to NO destabilizes HIF-1α and promotes its proteosomal degradation. Therefore, a pair of MCT sets m3 and m10 existing side by side would suggest that two contradictory processes coexist. In our model, the inconsistencies of this type do not occur, indicating the correctness of the model.

During hypoxia, HIF-1 may be stabilized in several ways and S-nitrosylation (t69) plays an important role in this process. It is regulated by the inhibitors of the caveolin-1 peptide scaffolding activity toward eNOS. In the supports of many of t-invariants, S-nitrosylation coincides with the formation of caveolin-1 in hypoxia, showing that the process of HIF-1 stabilization, crucial for cellular adaptation to hypoxic stress, is strictly controlled in mammalian cells to maintain proper oxygen homeostasis.

The analysis of the key phenomena in our model showed that the formation of NO and S-nitrosylation, together with hypoxia, are extremely important processes. The increased levels of NO inhibit angiogenic process [52], but low doses of NO stimulate this process, and this molecule is crucial for the initial steps of angiogenesis. This leads to the conclusion that ECs may have an additional anti-angiogenic mechanism based on NO concentration switch.

Our investigations were able to lead to the identification of the interactions between TGF-β signaling pathway activation (m7), NO synthesis (t18), the processes that stimulate NO synthesis (t71), and S-nitrosylation (t69). These relationships have been determined, among others, in t-invariants x16, x17, x27, and x28. Additionally, a number of t-invariants contains relationships such as m7, t18, and t71 (x18, x19, x29, x30) or m7, t18, and t69 (x37 and x38). An interesting summary of these observations is cluster c11, were t-invariants x18, x19, x29, and x30 coexist. Supports of these invariants contain transitions that correspond to the synthesis of TGF-β1, the increase in NO levels, and the reactions of S-nitrosylation, in a condition switching from normoxia to hypoxia. In conclusion, TGF-β1 may participate in the inhibition of angiogenesis, through the upregulation of eNOS expression. Furthermore, the analysis of the contents of x7, x9, x10, and x11 clusters showed that there is a distinct relationship between S-nitrosylation and HIF-1 degradation, leading to a conclusion that NO affects HIF-1 stabilization under hypoxic conditions, which is consistent with the results obtained by Callapina et al. [53].

Additionally, in supports of many t-invariants (x7, x8, x20-x30, x36-x41, x45-x48), the transitions that represent opposing processes coexist, e.g., hypoxia (t12) and oxidation (t65). This may be considered as a paradox, it may also be explained as an adaptive mechanism of the cells to hypoxic conditions.

Endothelial NOS, a critical regulator of cardiovascular homeostasis, whose dysregulation leads to different vascular pathologies, requires L-arginine and oxygen as substrates for NO synthesis. Our suggestion that, in hypoxic conditions, oxygen demand in the mitochondria of ECs is lower, allowing eNOS to produce higher quantities of NO, which in turn improves vasodilatation and oxygen uptake, may explain the adaptive mechanism of the cells. This can be sufficient only for short-term oxygen flow arrest because a prolonged lack of oxygen leads to cell injury and death. In future, the investigations should focus on how endothelial cells react in normoxia, compared with anoxia (a total depletion of oxygen, an extreme form of hypoxia), rather than hypoxic conditions.

Knockout analyses

The results of these analyses, sorted according to the importance of the knocked-out MCT set or single transition, are presented in Table 8. This type of knockout is based on the analysis of t-invariants only.

Table 8

The impact of net element knockout depending on the affected t-invariants.

MCT-set / transition	Biological function	Affected t-invariants
t12	The processes leading to low oxygen concentrations in an organism	85.42%
m2	The activation of the signaling pathway involving mitogen-activated kinases	79.17%
t18	NO synthesis	66.67%
m8	The activation of eNOS modulated by CaM	66.67%
m11	Sirt-1 dependent HIF-2 regulation	66.67%
m9	The stabilization of new vessels by angiopoietins	60.42%
m7	The activation of TGFβ signaling pathway	60.42%
m1	The stimulation of neovascularization caused by the selected factors with mitogenic activity in endothelial cells	58.33%
t30	The activation of eNOS through the interactions with CaM and Ca2+ complex	58.33%
t65	Oxygenation	54.17%
m5	HIF-1α accumulation under hypoxic conditions	45.83%
t16	The normalization of oxygen status in the organism	45.83%
t71	The processes that lead to the increase in NO	45.83%
t73	Normal cellular state	45.83%
t69	S-nitrosylation	33.33%
m6	Shear stress-mediated CaM-caveolin-1 cleavage	31.25%
t51	The synthesis of factors that stimulate vessel permeability in various ways	29.17%
m4	EPO pathway regulation	29.17%
t27	Endothelial cells proliferation and migration	29.17%
t59	HIF-2α influence on EPO	22.92%
m3	The degradation of HIF-1α under normoxic conditions	16.67%
t70	PHD activity influenced by NO	12.50%
t64	HIF-1α proline hydroxylase activation	4.17%
t49	TGF-β1 synthesis	2.08%
m10	The regulation of caveolin-1 expression in caveole	2.08%
t40	The remodeling of vessels	2.08%

The second kind of knockout analysis, as described in Methods section, is based on simulation of transition firings in different conditions. The results presented in S2 and S3 Tables have been obtained using our own tool, and they show the impact of manually disabled transitions. The knockout impact of a transition on the net can be presented graphically in the net structure. The graphical representation of the results found in S2 and S3 Tables is presented in Figs 4 and 5, respectively. The transitions with a small thunder icon have been manually disabled (they are considered as permanently inactive by the simulator, as explained earlier), the transitions marked with a red circle (but without the icon) have been knocked out as a result. The colors and a degree of filling of transitions and places represent the average chance of firing for a transition or an average number of tokens accumulated in places during the simulation. For the places this is represented in a separate bars located on the upper right side of each place, while for the transitions there is additional fraction, indicating the average chance of firing. For example, if a graphical representation of a transition is mostly painted green and the fraction is 0.358705, it means that its average chance of firing during the simulation is about 35.87%. Yellow and red filings represent a much lower average chance of firing.

Fig 4

Graphical representation of the knockout results for the entire net, upon the disabling of transition t1.

Fig 5

Graphical representation of the knockout results for the entire net, upon the disabling of transitions t12, t18, and the entire MCT set m2.

The first column in the Tables represents the results obtained following the disabling of transition t1 (”OFFLINE” status in the second row and the second column of the table). All transitions with the status “KNOCKOUT” have been knocked-out (the firing was completely inhibited) as a direct result of t1 inability to fire during the simulation. In other words, the offline status means that a transition has been manually set as permanently inactive no matter the tokens in its pre-places, while the knockout status means that a transition could not fire because there were not enough tokens in its pre-places to activate it during all the repeated simulations. The values in the last three columns of the table represent the Reference set, Knockout set, and Difference column, as described in the Methods.

The average chance of firing for many transitions changes depending on the conditions, and for example, transition t0 had an average chance of firing equal to 14.98% in the Reference set, i.e., in simulations where there had not been any manual disabling of invariants. When t1 was disabled, t0 average chance of firing dropped by about 3.82%, leading to the average chance of firing that equaled 11.17% (Knockout set). S3 Table shows the results obtained when t12, t18, and the whole MCT set m2 have been disabled. The word “disabled” in the last column instead of the percent value helps in a quick identification of rows with transitions either being disabled manually (“offline”) or being knocked-out.

Knock-out analyses conclusively show that the most important processes during angiogenesis are the ones leading to lower oxygen concentration (t12), the activation of the signaling pathway involving mitogen-activated kinases (m2) and NO synthesis (t18). The blocking of the transition/MCT set that corresponds to these processes leads to the blockage of 85.42%, 79.17% and 66.67% of all t-invariants in the investigated network, respectively. Since t-invariants correspond to subprocesses involved in angiogenesis, their blockage prevents the formation of new vessels. As previously mentioned, a knockout of an MCT set means a knockout of any transition that represents the element of this set. Therefore, when a knockout of MCT set m2 influences 79.17% of t-invariants it means that every transition belonging to this set (t0, t1, t2, t3, t4, t5, t7, t14, t20, t50, and t68) has the same impact on the net.

Simulation-based knockout analyzes, where transition t1 or transitions t12 and t18, and an MCT set m2 have been knocked-out (see S2 and S3 Tables) have led to some interesting observations. If transition t1 is knocked-out, this leads to the reduction of NO synthesis, reduced activation of eNOS and the enhancement of the processes that lead to an increase in EPO levels. Knockout of the transitions that correspond to hypoxia (t12), NO synthesis (t18) and activation of the signaling pathway that involves mitogen-activated kinases (m2), a vast decrease in eNOS activation and processes involving calmodulin and calcium pathways has been observed. Additionally, many transitions were disabled, demonstrating the extent of the influence of the excluded processes on the process of angiogenesis.

Final summary

A particularly important result of the presented research is finding a direct link (close connection) between molecules, such as TGF-β1, eNOS, NO and HIF-1 that play a crucial roles during angiogenesis. We have shown that TGF-β1 may participate in the inhibition of angiogenesis through the upregulation of eNOS expression, which is responsible for catalyzing NO production. The results obtained in the previous studies, concerning the effects of NO on angiogenesis, are inconclusive. These studies have suggested that NO inhibits EC migration, which is an essential step in angiogenesis [54]. More recent studies show that NO stimulates angiogenesis in vitro and in vivo, influencing the migration and differentiation of capillary EC. However, the results presented here are simply a starting point for further studies. A wet lab experimental validation of the biological interpretations of mathematical analysis results is necessary. These experiments could provide additional data, which could be used for the improvement of our model. Another direction of the future research should be an extension of this model, through the inclusion of quantitative data. This may be done, for example, by adding weights on the arcs, which correspond to quantities of substances that are involved in the elementary subprocesses occurring in the system. Additionally, a variant of Petri net may be used (e.g., timed, continuous, or hybrid Petri net) to model time dependencies or reaction rates.

Scripts for the R language performing clustering according to the chosen distance, linkage and cluster number.

The results are given in PDF files containing MSS evaluation for clusters.

(PDF)

Click here for additional data file.

The combinations of seven clustering algorithms and eight distance metrics are presented, resulting in different MSS evaluation values for the resulting clusterings.

(XLS)

Click here for additional data file.

Knockout results for all transitions following the disabling of transition t1.

(PDF)

Click here for additional data file.

Knockout results for all transitions following the disabling of transitions t12 and t18, and MCT set m2.

(PDF)

Click here for additional data file.

File in Snoopy representing the discussed net.

(SPPED)

Click here for additional data file.