Interrelations between Iron and Vitamin A—Studied Using Systems Approach

A deficiency of vitamin A (VAD) and iron is the most common nutritional problem affecting people worldwide. Given the scale of the problem, the interactions between vitamin A and iron levels are widely studied. However, the exact mechanism of the impact of vitamin A on the regulation of iron metabolism remains unclear. An extremely significant issue becomes a better understanding of the nature of the studied biological phenomenon, which is possible by using a systems approach through developing and analyzing a mathematical model based on a Petri net. To study the considered system, the t-cluster analysis, the significance analysis, and the analysis of the average number of transition firings were performed. The used analyses have allowed distinguishing the most important mechanisms (both subprocesses and elementary processes) positively and negatively regulating an expression of hepcidin and allowed to distinguish elementary processes with a higher frequency of occurrence compared to others. The analysis also allowed to resolve doubts about the discrepancy in literature reports, where VAD leads to positive regulation of hepcidin expression or to negative regulation of hepcidin expression. The more detailed analyses have shown that VAD more frequently positively stimulates hepcidin expression and this mechanism is more significant than the mechanism inhibiting hepcidin expression indirectly by VAD.


Research Context
Considerations regarding the role of vitamin A deficiency (VAD) in the regulation of hepcidin expression and, consequently, VAD effect on iron metabolism were (and still are) undertaken by many researchers [1][2][3][4][5][6]. The first literature reports on the relationship between vitamin A deficiency and low iron levels come from studies conducted in the 1920s [7][8][9]. Additionally, vitamin A and iron deficiencies are the most common nutrient deficiencies affecting millions of people around the world [1,3,5,6]. However, the mechanism of vitamin A influence on iron homeostasis was not fully understood [2,6]. Moreover, the study of these relationships raises a lot of controversy among researchers, which can be observed by the discrepancies in the literature reports. With the purpose to understand the rules underlying this phenomenon, a systems approach based on Petri nets [10][11][12] was used and described in this paper. This approach allows the systematization of knowledge and a better understanding of the nature of the biological system through creating a mathematical model being a representation of the complex network of interactions [13].

Biological Background
Iron is a unique trace element that is biologically essential, but also potentially toxic for humans, hence, maintaining its homeostasis is an extremely important issue. The essence of intracellular iron homeostasis is providing adequate amounts of iron to maintain the continuity of vital biological processes, as well as limiting the toxicity of Fe 2+ ions involved in the Fenton reaction [14].
Iron differs from other minerals because, due to the lack of a physiological excretion mechanism, its balance in the body is tightly controlled at the level of absorption. Iron ions in the form of Fe 3+ are reduced to Fe 2+ by duodenal cytochrome b (Dcytb), and in the next stage, they are transported to the enterocyte by the divalent metal transporter 1 (DMT1). Then, Fe 2+ ions can be absorbed by ferritin (Ft), or they can be transported into the circulation via ferroportin (FPN) using hephaestin (Heph). Heph oxidizes Fe 2+ ions to Fe 3+ , which are then absorbed by mammalian cells [15]. Newly absorbed Fe 3+ binds to plasma transferrin (Tf). Iron-loaded Tf binds to transferrin receptor 1 (TfR1) on the surface of most of the cells and forms a complex. The Tf-TfR1 complex is internalized via clathrin-mediated endocytosis [15]. Next, after endocytosis of this complex, iron enters the cytoplasm via DMT1. If iron is needed it may traffic directly to the sites of utilization, like mitochondria. If there is no requirement for iron, it may be sequestered for later use within ferritin-the iron-storage protein. The cell can also get rid of iron by exporting via FPN. This export pathway can act as a safety valve if very large amounts of iron accumulate in the cell. In enterocytes, the efficiency of cellular iron export is enhanced by Heph, but in most cells of the body, this role is played by the circulating Heph homolog, ceruloplasmin [15]. An important element is the cytoplasmic labile iron pool (LIP), indicating a deficiency or excess of iron. LIP in response to abnormal iron levels, activates the pathways regulating iron homeostasis [16]. Iron metabolism is controlled by two regulatory mechanisms: the first one acts at the cellular level, while the other works at the systemic level.
The regulatory mechanisms at the cellular level are based on the iron regulatory protein (IRP)-iron responsive element (IRE) signaling pathway. They involve cytoplasmic regulatory proteins 1 and 2 (IRP1 and IRP2) that bind to specific RNA sequences, iron-reactive sequences (IREs), thus regulating the expression of proteins involved in the maintenance of homeostasis, like TfR1, DMT1, Ft and FPN [15]. Depending on the cellular iron levels the IRP/IRE mechanism controls the synthesis of these proteins and their action. For the detailed IRP/IRE signaling pathway translation modulation mechanisms see [17].
Iron deficiency induces the formation of apo-IRP1, which binds to the IRE sequence in the 5'UTR region of the Ft and Fpn mRNAs, leading to an inhibition of their translation. On the other hand, the binding of apo-IRP1 to the IRE sequence in the 3'UTR mRNA region of TfR1 and DMT1 increases the encoded proteins. Otherwise, excess iron will induce the formation of holo-IRP1, which has an opposite effect. It inhibits the translation of TfR1 and DMT1 and increases the level of Fpn and Ft proteins [18].
The regulatory mechanism functioning at the iron system homeostasis level depends mainly on hepcidin (a hormone produced in the liver), the master regulator of systemic iron homeostasis [19]. The major stimuli regulating hepcidin transcription include iron concentrations in the blood, liver iron stores, inflammation, and erythropoiesis [20]. Hepcidin binds to FPN (in the enterocyte membrane), and as a consequence, the Janus kinase 2 (Jak2) is attached. It leads to the phosphorylation of tyrosine (Tyr) residues in FPN, which signals the protein translocation from the membrane into the cytoplasm. As a consequence, FPN is dephosphorylated, ubiquitinated to be degraded eventually. Because of the FPN degradation, iron becomes blocked in enterocytes, decreasing the amount of iron transported to the blood [18].
Hepcidin production is stimulated by iron loading and inflammation. Iron-mediated hepcidin regulation occurs via the bone morphogenetic protein (BMP) and SMAD (BMP-SMAD) pathway, whereas inflammation-mediated regulation occurs in both the IL-6 and JAK and signal transducer and activator of transcription proteins (STAT) (IL-6/JAK/STAT) signaling axis and via the BMP-SMAD pathway [19].
One of the hepcidin regulation mechanisms involves Tf and hemochromatosis (HFE)proteins, which compete for the binding of the membrane TfR1 [1]. In normal iron states or excess serum iron, Tf binds with a high affinity to TfR1, thus preventing the binding of HFE to TfR1. Free HFE sends signals to the hepatocyte nucleus increasing hepcidin expression.
In turn, in the case of iron deficiency, HFE binds to TfR1, reducing the expression of hepcidin [19].
Hepcidin expression can be also modulated by the inflammatory processes [20,21]. Overproduction of cytokines such as IL1β, TNF-α and IL-6 by macrophages and IFN-γ by lymphocytes weakens erytropoietin (EPO) synthesis, reduces the response to erythropoiesis, increases hepcidin levels and may activate erythrophagocytosis, especially in an acute form [19]. IL-6 induces the expression of the Hamp gene through the JAK and the signal transducer and activator of the transcription 3 (STAT3) pathway (JAK-STAT3). IL-6 is bundled up with JAK protein triggering its autophosphorylation and activation, resulting in phosphorylation of the transcription factor STAT3. In turn, STAT3 is translocated to the cell nucleus, where leads to the stimulation of the promoter of the hepcidin gene, increasing the expression of hepcidin [19].
The hemojuvelin-bone morphogenetic protein 6-and signaling protein 4 (HJV-BMP6-SMAD4) pathway is another way to modulate hepcidin expression [22]. HJV is an essential cofactor for the binding of bone morphogenetic protein 6 (BMP6) to the BMP receptor (BMPRs), leading to the phosphorylation of the signaling proteins SMAD 1, 5, and 8 (nuclear signal transducers). Phosphorylated SMAD 1, 5, and 8 proteins form a complex to bind the SMAD4 protein in the next step, which is a signal to move to the nucleus and stimulate transcription of the hepcidin gene [23]. This mechanism is triggered in the event of an excess of iron.
In iron deficiency, the HJV-BMP6-SMAD4 pathway is inhibited by furin, specifically, by the action of the soluble form of hemojuvelin (sHJV). sHJV is released from the cell membrane by furin and then binds to the BMP6 protein, preventing it from binding to the receptor. Consequently, the HJV-BMP6-SMAD4 pathway is inhibited, and thus, the concentration of hepcidin is reduced [16].
Besides, hepcidin expression is modulated by VAD. It turns out that vitamin A may increase or decrease hepcidin expression. VAD is characterized by an increase in the level of the pro-inflammatory cytokine, such as IL-6 and IL-1β, responsible for increasing hepcidin expression [1]. VAD also increases the level of the BMP6 protein, which is essential in the HJV-BMP6-SMAD4 pathway that stimulates the transcription of the hepcidin gene [2]. On the other hand, the level of the BMP6 protein is reduced in the iron deficiency state [2], reducing the hepcidin expression by preventing the HJV-BMP6-SMAD4 pathway. An equally important mechanism for modulating hepcidin expression is the impairment of erythropoiesis, both in the case of iron deficiency and vitamin A deficiency [2]. VAD leads to impaired erythropoiesis by lowering erythropoietin expression. It leads to erythrocyte malformations and, consequently, accumulation of the heme group in the spleen [2].
The described biological background is graphically presented in Figures 1 and 2, which relate to two different situations. Figure 1 is focused on the interactions between modeled subprocesses in an iron overload situation, while Figure 2 concerns the interactions between the modeled subprocesses in an iron deficiency situation.

Petri Net-Based Model
The purpose of creating a model based on Petri nets is to systematize knowledge and better understand the interrelations between vitamin A deficiency and iron deficiency. Although vitamin A deficiency affects the regulation of iron homeostasis, these mechanisms remain unclear. This fact became a motivation to focus on this biological phenomenon. The proposed Petri net-based model was created using Snoopy (open-source software) [24], while the analysis was performed using Holmes (free, stand-alone Java application) [25]. Properties of the proposed model are as follows: 68 passive components (places), 89 active components (transitions), 402 t-invariants (subprocesses). The Petri net-based model was illustrated in Figure 3, which includes only names of places, while names of transitions are included in Table 1 (to increase the readability of the model). The model is available in Supplementary Materials.

Analysis
The analysis of the Petri net model is based on a search for similarities between t-invariants (subprocesses), leading to an identification of subprocesses that may interact with each other. Such knowledge allows for a better understanding of the modeled system, and may also allow for a discovery of interesting properties of the investigated biological system. The analysis of the proposed model includes a Maximal Common Transition (MCT) sets analysis, a t-cluster analysis, a significance analysis with the completion by a knockout analysis, and an analysis of an average number of firings of transitions (simulation). These methods are described in Section 4.
In the first stage of the analysis, the MCT sets analysis was performed. The proposed model contains 89 transitions (elementary processes), which were grouped into non-trivial MCT sets (see Section 4). In a biological context, MCT sets correspond to certain functional modules. Therefore, a biological meaning for each MCT set was assigned, what is described in Table 2. All the described modules are also marked in Figure 3 as colored frames.

MCT Set Contained Transitions Biological Interpretation
Iron ion Fe 3+ is reduced to ion Fe II by Dcytb and then is exported to enterocytes (to LIP) by DMT1. Hepcidin binds to Fpn and leads to degradation and internalization of Fpn. What results in prevention of iron release from enterocytes and macrophages.     Vitamin A deficiency and iron deficiency impair of erythropoiesis, which results in increased phagocytosis of malformed and undifferential erythrocytes. This mechanism in consequence leads to accumulation of iron in spleen and to decrease of Hamp mRNA level.

MCT Set Contained Transitions Biological Interpretation
Iron is engaged in Fenton reaction, which leads to increase of oxidative stress via ROS. m 10 t 29 , t 32 Iron absorption by mammals via endocytosis. To be precise, iron ion Fe 3+ binds to transferrin (Tf) and next binds to transferrin receptor protein 1 (TfR1), which results in iron absorption.    In the next stage, clustering of 402 t-invariants was performed. To assess the best clustering Mean Split Silhouette index (MSS) and Calinski-Harabasz (C-H) coefficient were used [22]. The first coefficient, MSS, evaluates the fit of each t-invariant to its cluster and the average quality of a given clustering [26,27]. The second one, C-H, allows indicating the optimal number of clusters by the highest value [28]. The set of t-invariants was divided into 20 t-clusters using Average Linkage method and Pearson distance measure. Each of these clusters consists of certain biological subprocesses, which are described in Table 3. Moreover, occurrences of the particular subprocesses in a given t-cluster are presented in Figure 4 (the more intense the color, the given subprocess appears in the greater number of t-clusters). Subprocesses, which may play a significant role in the modeled system, can be found through their frequency of occurrence in t-clusters.
The purpose of the t-cluster analysis is to find subprocesses that may play a significant role in the functioning of the modeled system. It is possible to identify subprocesses that may be more crucial than others. It seems that if a certain subprocess occurs in a larger number of clusters, it may be stimulated by various, independent processes. Due to this fact, it may be concluded whether this subprocess is stimulated more or less frequently. However, the analysis of t-clusters is quite general, and such conclusions need to be confirmed by applying more precise analyzes. For example, an essential mechanism of the investigated system is the regulation of hepcidin expression. The presented model includes several independent pathways involved in hepcidin expression, including HJV-BMP6-SMAD4 pathway (subprocess (n)), free HFE (subprocess (f)), and JAK-STAT3 pathway modulated by cytokines IL-6 and IL-1β (subprocess (h)). Therefore, the frequency of these subprocesses in all t-clusters was considered. As shown in Figure 4, subprocess (n) occurs in 6 t-clusters, subprocess (f) occurs in 6 t-clusters, and subprocess (h) occurs in 13 t-clusters. Therefore, it seems that hepcidin expression is more frequently induced by cytokines (subprocess (h)) than other pathways. However, this type of conclusion requires additional confirmation by conducting more detailed analyzes. The results of the t-cluster analysis are described in more detail in Section 3.  Table 2, where information about transitions contained in each MCT set is also included. Table 3. List of subprocesses included in t-clusters.

ID Biological Meaning
(a) Iron is exported by Fpn to blood circulation and then it can be absorbed by mammalian cells via endocytosis. To be precise, iron ion Fe 3+ binds to transferrin (Tf) and next binds to transferrin receptor protein 1 (TfR1), which results in iron absorption.
(b) Iron can be engaged in Fenton reaction, but it occurs rarer than iron absorption in enterocytes.
(c) IRP/IRE mechanisms play important roles: in case of increase of iron concentration holo-IRP1 is formed and it leads to increase of Fpn and to decrease of DMT1 and TfR1.
(d) Iron is exported to enterocytes by DMT1, where is absorbed by Ft.
(e) IRP/IRE mechanisms play important roles: in case of decrease of iron concentration apo-IRP1 is formed and it leads to increase of DMT1 and TfR1 and to decrease of Fpn.
(f) Expression of hepcidin by free HFE (in case of iron overload or iron normal status Tf bind to TfR1, which prevents HFE binding).
(g) In case of iron deficiency HFE binds to TfR1, which lead to decrease of hepcidin expression.
h) Expression of hepcidin by cytokines IL-6 (JAK-STAT3 pathway) and IL-1β (i) Vitamin A deficiency additionally stimulates gene Hamp expression by increase of IL-6 and IL-1β.
(j) Hepcidin binds to Fpn and leads to degradation and internalization of Fpn. This mechanism results in iron accumulation in enterocytes and macrophages and in consequence it leads to decrease of iron concentration in serum.
(l) Inflammation leads to hepcidin expression by cytokines.
(m) TfR1 is necessary to iron absorption by mammalian cells. (t) Vitamin A deficiency and iron deficiency impair of erythropoiesis, which results in increased phagocytosis of malformed and undifferential erythrocytes. This mechanism in consequence leads to accumulation of iron in spleen and to decrease of Hamp mRNA level.
As mentioned before, the analysis of t-clusters may turn out to be quite general and thus may cause some phenomena to be ignored and not noticed. Therefore, it is also worth using more detailed analyzes, such as a significance analysis, a knockout analysis, or an analysis of the average number of transition firings (simulation). Conducting additional analyzes will allow to confirm the results obtained from the cluster analysis.
The significance analysis was performed for elementary processes related to hepcidin expression, both stimulation and inhibition. The results of the significance analysis for selected transitions are presented in Table 4. This table contains subprocesses that stimulate hepcidin expression and subprocesses that inhibit hepcidin expression (all subprocesses are listed in the "Subprocesses" column). Each subprocess consists of at least one elementary process; the transitions' IDs and the names of elementary processes are included in the columns "ID" and "Name of elementary process", respectively. The column, "Significance" is divided into "Frequency trans./t-inv.", which includes information about the occurrence frequency of a given transition in supports of t-invariants, and "Percentage ratio" which includes occurrence frequency expressed as percentage, hereinafter called also significance of given transition. In addition, Table 4 includes also the results for the analysis of an average number of the firing of transitions, which are included in the last column "AvgT". This analysis allows obtaining information about the frequency of transition activation. For example, transition t 59 (corresponding to hepcidin expression via the JAK-STAT3 pathway) is part of the larger subprocess of stimulating hepcidin expression by the cytokines IL-6 and IL-1β. This transition occurs in 252 supports of t-invariants which consists 63% of all t-invariants (402 t-invariants). The significance analysis is related to knockout analysis, which means that a knockout of transition t 59 excludes 63% of all modeled subprocesses. Additionally, the average number of firings of transition t 59 is equal to 46.94. This value is one of two the highest values among all elementary processes positively regulating hepcidin expression.  The significance analysis allows determining which elementary process is more significant for the functioning of the whole model. For example, based on Table 4, among the elementary processes that stimulate hepcidin expression, transitions t 59 and t 60 are the most important (these transitions correspond to elementary processes which positively regulate hepcidin expression by cytokines). In this case, these results confirm the conclusions of the t-cluster analysis, which has shown that hepcidin expression induced by cytokines (subprocess (h)) may be more significant than induction by other subprocesses. Moreover, another more detailed analysis is the analysis of the average number of firings of transitions. Contrary to the significance analysis, this is not a structure-based but a simulation-based analysis. Such simulation analysis allows obtaining information about the frequency of firings of a particular transition. For example, based on Table 4, transition t 60 is the most frequently fired among all elementary processes positively regulating hepcidin expression. Therefore, the elementary subprocess which is a positive regulation of hepcidin expression by cytokines occurs more frequently than other elementary processes positively regulating hepcidin. A detailed comparison of the results obtained from the mentioned analyses for the mechanisms stimulating and inhibiting hepcidin expression is described in Section 3.

Discussion
The proposed model is focused on the selected aspects of the interaction between vitamin A and iron, where different levels of iron concentration was taken into account (iron deficiency, normal iron status, iron overload). The performed analyses made it possible to study the interactions between particular subprocesses, and in consequence, a better understanding of the modeled phenomenon.
The analysis of t-clusters focused on the role of hepcidin in the regulation of iron homeostasis. There are several independent mechanisms and biological components influencing positively the expression of hepcidin, which are as follows: 1.

3.
Free HFE; in the case of iron overload or normal status Tf binds to TfR1 and prevents HFE binding (subprocess (f)) included in 6 t-clusters, see Table 4).
Moreover, the cluster analysis showed that hepcidin expression induced by IL-6 and IL-β may be more significant than induction by other subprocesses. This fact can be associated with vitamin A deficiency, which additionally stimulates inflammation through an increase of the mentioned cytokines. However, other subprocesses, like subprocess (n) corresponding to HJV-BMP6-SMAD4 pathway and subprocess (f)) corresponding to the situation in which there exists free HFE, occur in the same number of clusters. Therefore, the role of these subprocesses should be equally important for the functioning of the system. These results were compared to the results obtained from the more detailed analysis, which is the significance analysis (see Table 4). The results of the significance analysis showed that the most significant elementary processes are transitions associated with induction of hepcidin expression by cytokines (t 59 , t 60 with a significance of 63%). On the other hand, the elementary processes related to the stimulation of hepcidin expression by HJV-BMP6-SMAD4 pathway and by free HFE have a similar significance level (t 73 -29%, t 64 -26%, respectively). Thus, the results of the significance analysis confirm the results of the tcluster analysis for subprocesses related to the positive stimulation of hepcidin expression. The t-cluster analysis confirm the above results, because it has shown that subprocess (h) corresponding to induction of hepcidin expression by cytokines may be the most significant compared to the others, and subprocesses (n) and (f) although less crucial than subprocess (h), may be equally important for the functioning of the system. In addition to estimating the significance of selected transitions, an average number of firings of transitions was also considered (see the column "AvgT " in Table 4). As can be seen, the most often fired are transitions t 59 , t 60 and t 82 corresponding to stimulation of hepcidin expression through cytokines, next in order is transition t 73 corresponding to stimulation by HJV-BMP6-SMAD4 pathway, and finally transitions t 55 and t 64 corresponding to the presence of free HFE. Despite the fact that the last two elementary processes are equally important in the sense of the structure of the entire model, they may occur with different frequencies.
The next stage of the analysis of the results is focused on the mechanisms leading to inhibition of hepcidin expression, which are as follows: 1.
The t-clusters analysis revealed that among the subprocesses inhibiting hepcidin expression, inhibition of HJV-BMP6-SMAD4 pathway and impairing of erythropoiesis are more relevant than the other mechanisms. HJV-BMP6-SMAD4 pathway can be inhibited by two different mechanisms. The first one is related to an activation of furin in case of iron deficiency. Furin releases HJV from a membrane as a result of a proteolytic reaction. Soluble sHJV blocks BMPRs receptors resulting in inhibition of HJV-BMP6-SMAD4 pathway. The second one is a reduction of BMP6 (in case of iron deficiency and VAD), which is an essential element in HJV-BMP6-SMAD4 pathway.
Similarly, as for the subprocesses related to the stimulation of hepcidin expression, a more detailed analysis was performed for the inhibitory mechanisms. The significance analysis has shown that the most significant elementary process is the transition associated with inhibition of HJV-BMP6-SMAD4 pathway (t 85 with 19% significance), next in order are elementary processes related to HFE binding to TfR1 (t 62 and t 63 , both with 17% significance), followed by elementary processes related to SMAD7 induction by BMP6 (t 74 and t 75 with 10% significance), and finally elementary process corresponding to impairing of erythropoiesis (t 85 with 9% significance). These results of the significance analysis are not consistent with the t-cluster analysis. The t-cluster analysis has shown that subprocess (o) corresponding to inhibition of HJV-BMP6-SMAD4 pathway and subprocess (t) corresponding to impairing of erythropoiesis seem to be the most significant, while only the first one is really crucial. Moreover, subprocess (t) has the slightest importance in the significance analysis of the selected elementary processes. As mentioned earlier, the t-cluster analysis may be too general. It turns out that the most significant mechanisms inhibiting hepcidin expression are the inhibition of HJV-BMP6-SMAD4 pathway, as well as the binding of HFE and TfR1. Apart from the significance analysis, the simulation of the average number of firings was also performed. The fact that a given elementary process is significant from the point of view of the system structure does not mean that it is fired most often. In Table 4 in the column "AvgT " it can be seen that induction of SMAD7 by BMP6 is most often fired (this mechanism indirectly inhibits HJV-BMP6-SMAD4 pathway, which requires the presence of BMP6). However, transitions t 74 and t 75 corresponding to the mentioned subprocess are not so much significant for the structure of the modeled system. The other mechanisms of inhibition of hepcidin expression (i.e., inhibition of HJV-BMP6-SMAD4 pathway or binding HFE to TfR1), which are structurally important, occur with much lower frequency than induction of SMAD7 by BMP6.
An important aspect of the model analysis is the effect of vitamin A deficiency and its relationship with iron. The cluster analysis identified a few subprocesses directly related to vitamin A deficiency based on the frequency occurrence of subprocesses in t-clusters:

3.
VAD and low iron level leads to decrease of BMP6, which results in inhibition of hepcidin expression (subprocess (s) included in 12 t-clusters).

4.
VAD and iron deficiency impair of erythropoiesis leading to negative regulation of hepcidin expression (subprocess (t) included in 12 t-clusters).
The results of the t-cluster analysis (the frequency occurrence of the mentioned above subprocesses in t-clusters) revealed unquestionable influence of vitamin A deficiency and iron deficiency on inhibition of hepcidin expression. There exist different literature reports concerning the influence of vitamin A deficiency on hepcidin expression. The results presented in [1] clearly describe that vitamin A maintains iron homeostasis by positively modulating hepcidin expression in the liver. On the other hand, the results reported in [2] indicate that vitamin A deficiency also promotes reduction of hepatic Hamp mRNA levels, rather than the expected up-regulation of Hamp gene expression.
The experimental procedures of these two studies ( [1,2]) leading to the various results are briefly summarized below: • Studies described in [1]: Given these discrepancies (between the results presented in [1,2]), two different subprocesses were considered in the next stage of the analysis. The first one, in which vitamin A deficiency has an influence on the increase of BMP6, leading to hepcidin expression via HJV-BMP6-SMAD4 pathway (see subprocess (p)) in Figure 4). The second one, in which BMP6 induces SMAD7, which negatively regulates hepcidin (see subprocess (k)) in Figure 4). These two opposite mechanisms are included in the same number of t-clusters. Therefore, it seems they may be equally important for the modeled system, which does not solve the problem of existing different literature reports.
The results of the t-cluster analysis may be too general, thus performing clustering as the only analysis may lead to an omission of some meaningful relationships. Nevertheless, clustering may narrow the area of search and suggest a further direction of research. Therefore, the analysis of the average number of transition firings was performed for the mentioned before subprocesses (k) and (p).
The results of the analysis of the average number of firings of selected transitions are presented below: This analysis showed that the elementary process associated with activating HJV-BMP6-SMAD4 pathway (which stimulates hepcidin expression) is more often fired than the subprocess leading to a negative regulation of hepcidin expression by involvement of BMP6 in SMAD7 induction.

Methods
Petri nets are mathematical objects that in a natural way can be used for modeling and analysis of systems composed of concurrent processes. Their origins are in the area of theoretical computer science, and for decades they were used mainly for modeling technical systems. However, during the last two decades, it appeared that they are very well suited for modeling and analysis of biological systems [10][11][12][13].
One of the reasons of this suitability is a structure of Petri nets, being a directed weighted bipartite graph. It means that a Petri net is composed of two disjoint sets of vertices. Vertices from one of these sets are called transitions and they usually correspond to some elementary active components of a modeled system (e.g., chemical reactions), while vertices from the other set, called places, are counterparts of elementary passive components of the system (e.g., substrates or products of chemical reactions). The vertices can be joined by an arc but only vertices of various types can be connected in this way. The arcs correspond to some causal relations between passive and active components of the net.
However, a Petri net is not a graph-there is one more, very important kind of components, i.e., tokens. They reside in places and can flow from one place to another through transitions. This flow is a crucial feature of Petri nets since it corresponds to flow of information, substances, signals, etc. through the modeled system. It is governed by a simple principle, called transition firing rule. According to it a transition becomes active if in every of places directly preceding it (such places are called pre-places of this transition) the numbers of tokens are equal to at least the weight of an arc joining a given place with the transition. An active transition can be fired, what means that tokens flow from its preplaces to its postplaces (i.e., places directly succeeding the transition), wherein the numbers of flowing tokens are equal to the weights of respective arcs.
Petri nets, as mathematical objects, can be described and analyzed using mathematical methods but they also have an intuitive and useful graphical representation. In this representation, transitions are depicted as rectangles or bars, places as circles, arcs as arrows, tokens as dots, or positive integer numbers located in places and weights as numbers labeling arcs (when a weight is equal to one, it is usually omitted in this representation for simplicity). The graphical representation is very helpful in understanding a structure of the modeled system and especially useful at the stage of developing the model and simulating it [10,13,22,30].
However, despite that the graphical representation of Petri nets is very useful, it is not well suited for a formal analysis of them. Hence, another representation is also often used, i.e., an incidence matrix. Such a matrix A is composed of n rows corresponding to places and m columns corresponding to transitions. Entry a ij of this matrix is a number equal to a difference between the numbers of tokens present in place p i before and after firing transition t j .
An analysis of a Petri net-based model of a biological system can be based on tinvariants. Such an invariant is vector x being a solution of the following equation: The length of vector x is equal to m, i.e., the number of transitions, and each entry in this vector correspond to a transition. With t-invariant x there is associated a set of transition called its support. The elements of the support are transitions which correspond to positive entries in x. More formally, a support of t-invariant x is set s(x) = {t j : x j > 0, j = 1, 2, . . . , m}.
t-invariants correspond to subprocesses, which can be especially important for the functioning of the modeled system. If every transition t j ∈ s(x) is fired x j times, then a distribution of tokens over the set of places (such a distribution is called a marking of a Petri net) does not change. The state of the modeled system becomes unchanged. Hence, t-invariants correspond to subprocesses which do not change the state of the modeled biological system [22]. Moreover, the Petri net model should be covered by t-invariants, i.e., each transition should belong to a support of at least one t-invariant. In a biological context, every elementary biological process modeled by a transition contributes to a behavior of the biological system. Therefore, the analysis of t-invariants may lead to a better understanding of interactions between the subprocesses of the modeled biological system. However, there is an assumption that the analyzed biological system is in a steady state.
There is also another type of sets of transitions important in the analysis of Petri net-based models of biological systems, i.e., Maximal Common Transition sets (MCT sets). Elements of such a set are transitions, which exclusively belong to supports of exactly the same t-invariants. From this follows that MCT sets divide the set of all transitions into disjoint subsets (i.e., every transition belongs to exactly one MCT set). It should be noted that some of MCT sets may contain only one transition-they are called trivial MCT sets. In general, MCT sets (except the trivial ones) correspond to some functional modules of the modeled biological system [13,22,31].
Apart from that transitions can be grouped into MCT-sets, also t-invariants can be grouped into sets called t-clusters. They contain those t-invariants which are similar to each other according to some similarity measure. Such clusters can be determined using standard clustering algorithms. However, there are many algorithms of this type and there are also many similarity measures which can be used. So, despite that to determine a collection of t-clusters standard algorithms and measures can be used, the task is not a trivial one, since the algorithm as well as the similarity measure should be properly chosen. Moreover, the number of clusters should also be determined in a proper way [21,32].
In addition to the mentioned analyzes, more detailed analyzes can be also carried out, i.e., a significance analysis with the completion by a knockout analysis and an analysis of an average number of transition firings [33]. The significance analysis allows determining a significance of a given elementary process based on the occurrence frequency of a transition (corresponding to such an elementary process) in all supports of t-invariants. The significance of a given transition is a ratio of the number of t-invariant supports containing this transition to the number of all t-invariants and it is expressed in percents.
The significance analysis is directly associated with the knockout analysis, which evaluates the significance of a given transition (or a set of transitions) based on the quantity of excluded subprocesses as a result of a knockout. A current number of t-invariants can be determined once again after the knockout of a selected transition. On this basis, it may be calculated how many subprocesses were excluded in consequence of the knockout. The more subprocesses were excluded, the more important the given elementary process is.

Conclusions
Using a systems approach based on Petri nets allows the analysis of biological phenomena as systems characterized by a complex net of interactions. To study the presented system in detail, the t-cluster analysis, the significance analysis, and the analysis of the average number of transition firings were performed. The above analyzes allowed to distinguish the most important mechanisms positively stimulating the expression of hepcidin and the most important mechanisms negatively regulating the expression of hepcidin. In addition, apart from distinguishing the most significant subprocesses/elementary processes for the functioning of the model, the average number of firings of particular transitions were also determined (it means that the frequency of occurrence of individual elementary processes).
In the case of the studied system focusing on the effect of vitamin A deficiency on the regulation of iron homeostasis, a thorough understanding of the biological problem is difficult due to contradictory literature information. The study of the proposed model aims to understand the analyzed system's nature better and systematize the knowledge. The final analysis also allows resolving some doubts as to the discrepancy in literature reports, where VAD leads to negative regulation of hepcidin expression [2] (in this paper, the authors emphasize the fact that the achieved results are different than expected) or to positive regulation of hepcidin expression [1]. A graphical representation of this problem is shown in Figure 5, where the mechanism in which VAD negatively regulates hecidin expression is associated with induction of SMAD7 by BMP6, while the mechanism in which VAD positively regulates hecidin expression is associated with stimulation of HJV-BMP6-SMAD4 pathway. The t-cluster analysis showed that these two opposite subprocesses may be equally important for the functioning of the modeled system. However, the more detailed analyses showed differences. Both the significance analysis as well as the analysis of the average number of transition firings indicated that VAD more frequently positively stimulates hepcidin expression, and this mechanism is more significant than mechanism inhibiting hepcidin expression indirectly by VAD.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The