A structural analysis of the hypoxia response network

Datasets

This study focused on the structure of the hypoxia-inducible factor-1 (HIF-1) (Semenza, 1999; Pugh & Ratcliffe, 2003) signaling pathway of the HRN and retrieved related information from the entries in the Kyoto Encyclopedia of Genes and Genomes (KEGG) database (Kanehisa & Goto, 2000). KEGG provides data resources related to the high-level functions and operations of cells, organs, ecosystems, and other life systems based on molecular-level information, especially information used for parsing molecular datasets generated by gene sequencing and other high-throughput experimental techniques (Kanehisa et al., 2017). There are two essential elements of KEGG discussed in this paper: the KEGG PATHWAY and KEGG ORTHOLOGY databases. The KEGG PATHWAY database is a collection of knowledge-based hand-drawn maps of pathways based on molecular interaction networks covering metabolic, genetic information processing, environmental information processing, cellular processes, biological systems, human diseases, and drug discovery. KEGG PATHWAY provides a concrete explanation of the meaning of the markers of the HIF-1 signaling pathway. The KEGG ORTHOLOGY (KO) database is a set of ortholog assemblages manually defined to denote the nodes (boxes) in the KEGG PATHWAY maps. The distinctive identifier that determines each KO entry is referred to as the K number (‘K’ followed by a five-digit number).

We analyzed the data in the KEGG Markup Language (KGML) file. These data contain information that corresponds to the HIF-1 signaling pathway in KEGG PATHWAY. The KGML, which enables the automated mapping of KEGG pathways, is an extensible markup language (XML) representation of the KEGG pathway. It is conducive to the computer-aided analysis and model building of gene/protein and chemical networks. There are two types of graphic elements in the KGML-based metabolic pathways. One is a rectangle, which represents an enzyme, connected by a “relationship.” The other is a circle, which represents a compound, connected by a “reaction.” In non-metabolic pathways, there are only rectangular elements, which indicates that there are only proteins joined by “relationships” in these pathways. The HIF-1 signaling pathway addressed in this study belongs to the metabolic pathways.

Overall analytical methods

The methods of this study are based on network analysis. The network of the HIF-1 signaling pathway is a metabolic network. It is a directed unprivileged network; i.e., the direction of a line between two nodes needs to be considered rather than the weight in the network. The corresponding analytical methods applied in this study are as follows: Process the text in the KGML file of the HIF-1 signaling pathway. Refine the network topology information and transform it into a network file that is identifiable for network analysis. Analyze the network indicators and draw the network topology maps. Analyze and test the structural complexity of the network. The generation of the network file is based on regular expression using the text-editing software Notepad++ and functions using the spreadsheet software Excel. The process of network analysis is based on network analysis models using Pajek. These methods and tools are comprehensively utilized to realize the qualitative and quantitative integration analysis of the network structural properties of the HIF-1 signaling pathway.

Biological network analysis tools are used to identify, analyze, visualize, or simulate nodes (reactants) and edges (reactions) from various input data types, including mathematical models of biological networks. There are four commonly used tools, namely Gephi, Networkx, IGraph, and Pajek. All of them are free for use, and they can handle large graph size. Gephi and Pajek are GUI based network analysis tools, whereas Networkx and IGraph are used in a programming language. We have chosen Pajek as the analysis tool based on the following reasons: First, Pajek supports a multi-relational network graph. There are different reactants in the HIF-1 signaling pathway, such as orthologous enzymes, compounds, and groups. And the kind of their relationships is diverse, including activation reactions, expression reactions, and inhibition reactions. Therefore, the network of the HIF-1 signaling pathway is a multi-relational network which can be analyzed by Pajek. Second, Pajek support partition feature enables the identification and extraction of the HIF-1 signaling pathway core network. Third, Pajek supports more network graph layout algorithms, and its algorithms for analyzing the corresponding network features are better. Forth, Pajek is a standalone software that is easy to learn and easy to use (Akhtar, 2014).

Topology

The network topology indicators include the node number, the line number, the line value, the network density, and the network diameter. The node is an entity that forms the network. The lines between the nodes represent the connections between network entities. The line value is the weight of a line, which indicates the strength of the relationship between network entities. The diameter of the network is the longest path of the shortest distance between pairs of nodes in the network, which is the number of maximum steps required to connect any pair of nodes in the network. The network density is equal to the numeric ratio of the actual connections to the possible connections. It indicates the degree of closeness of the relationships between the nodes in the network.

Centrality

Centrality is an essential concept of network analysis. A highly centralized network supports the convenient transmission of information. The central node has a critical influence on the transmission of information in the network. The metrics for centrality used in this paper include the centrality of the degree, the closeness, and the betweenness.

Cohesive subgroups

Some entities in the network relate to each other so tightly that they form a small local group known as a cohesive subgroup. The number of cohesive subgroups in the network as well as the local and global associations between subgroups are analyzed via cohesive subgroup analysis. This analysis includes component analysis, K-core analysis, and island analysis based on different analytical indicators and perspectives (De Nooy, Mrvar & Batagelj, 2018). The component analysis is mainly applied in this study. The component is an essential indicator of the cohesive subgroup analysis, which refers to the largest connected subnetwork in the network; i.e., there is a way to reach other nodes between any nodes in the subnetwork. Biometabolism networks often contain components that are not connected, and the most significant component is often the one that needs attention.

Bow-tie structure

The bow-tie property is an essential component of these properties. Broder et al. (2000) revealed that there is a bow-tie structure in the topology structure of the Internet at the macro level.

Figure 1 shows the bow-tie structure. The nodes of the input section can reach the nodes of the core section of a strong connection, while the reverse is not true. The nodes of the core section of the strong connection can reach each other nodes. The nodes of the core section of the strong connection can reach the nodes of the output section, while the reverse is not true. The nodes of the pipe section can reach the nodes of the output section via the input section. The tendril section connects with the input section and the output section. The nodes of the input section can reach the nodes in the same section or the nodes of the output section. There are no connections between the nodes of the non-connected component and other sections.

Ma & Zeng (2003) divided the macrostructure of the biological metabolic network into four parts: the giant strong component (GSC), substrate subset (S), product subset (P), and isolated subset (IS). These components correspond to the strongly connected core, input, output, and disconnected components of the Internet bow-tie structure, respectively, while the exact opposite situation is not observed. Any two metabolites within GSC can be generated from each other by a series of reactions. Any metabolite in S can be converted into the corresponding metabolite in GSC, but not vice versa. Any metabolite in P can be transformed via a series of the corresponding metabolites in GSC, while the exact opposite situation is not observed. The metabolites in IS cannot be converted from those in GSC, nor can they be converted into the corresponding metabolites in GSC.

Structural complexity

Qian et al. provided a strict definition of complex networks (Xuesen, Jingyuan & Ruwei, 1993). A network with some or all of the following characteristics is known as a complex network. These characteristics are self-organization, self-similarity, attractor, small world, and scale-free. We start with a scale-free perspective to verify the sophisticated attributes of the network. The scale-free property refers to the invariance of the network’s scale, meaning that the degree of the network nodes obeys the law of a power exponential distribution (i.e., the power-law distribution). It is also known as the Pareto distribution rate or the Zipf rate. This distribution law describes the distribution characteristics of the degree of network nodes. A small number of nodes in the network tend to present a large number of connections, but the number of connections of most nodes is very limited. This paper uses the Pareto distribution rate: (4)${\displaystyle P_r\left[X\ge x\right]\sim x^{-a}}$ where x is the probability corresponding to the actual measurement scale of the network, and a is the scale parameter of the power rate distribution.

The law of the degree distribution of the network nodes needs to be fitted to the power exponential distribution law to verify this characteristic of the core network of the HIF-1 signaling pathway. That is, the fitting value, $\hat{a}$, of the scale parameter, a, is obtained. Then, the characteristic is judged by the fitting effect.

The power-law fitting of the degree distribution

The general method of power-law fitting is to use the least-squares linear fitting based on double-logarithmic coordinates and to use the R² test to measure the fitting effect. However, Goldstein et al. believe that the scale result obtained by this method shows a significant error in relation to the corresponding value of the actual measurement network scale (Goldstein, Morris & Yen, 2004). Therefore, Clauset and Barabási et al. proposed a maximum likelihood estimation method for power rate fitting and tested the fitting results using KS statistics (Barabási, Albert & Jeong, 1999; Clauset, Shalizi & Newman, 2009). At present, no research shows that all the data corresponding to the actual measurement scale of the studied network obey the power-law distribution, but there is a critical value. The data corresponding to these scales obey the power-law distribution only when the value corresponding to the actual measurement scale, x, is higher than x_min. This study used Clause’s universal method, applicable to both discrete and continuous data, to estimate x_min. For actual network-scale data, the following formula is used to estimate the scale parameter, $\widehat{a}$, of the power-law distribution: (5)${\displaystyle \widehat{a}\simeq 1+n{\left[\sum _{i=1}^{n}ln\frac{x_i}{x_{min}-0.5}\right]}^{-1}}$ where x_min can be accurately obtained by calculating the maximum difference, L, between the actual measurement-scale data and the corresponding data of the fitted model: (6)${\displaystyle L=\underset{x\ge min}{max}\left|S\left(x\right)-P\left(x\right)\right|}$ where S(x) is the corresponding data of the actual measurement, P(x) is the corresponding data of the fitted power-law distribution model, and x_min, which minimizes L, is the optimal value.

Similarity test

This study used the K-S (Kolmogorov-Smirnov) method to test the “distance” L between the actual measurement data and the fitted power-law distribution model. The model constructed from actual measurement data is denoted as M. This model produces n sets of data. There is a set of data whose “distance”, L, from the corresponding fitted model, M, is greater than the “distance”, L, between the actual measured data and the fitted model, M. The number of such datasets is m. $\frac{m}{n}$ is represented as p and is known as the p value. If the p value is large (close to 1), it can be considered that the statistical fluctuations alone cause the difference between the actual measurement data and the fitted model. If the p value is small, there is room for adjustment in the rationality of the fitted model. If p ≤ 0.1, it can be considered that the actual measured data do not obey the power-law distribution (Ping, 2014).

Topology analysis

Figure 2 shows the overall topology of the HIF-1 signaling pathway network based on the KGML file from KEGG. Quantitative analysis of the network indicators is performed, and the results are as follows: The numbers of nodes and lines are 85 and 61, respectively. The network has a ring and ringless density of 0.00844291 and 0.00854342, respectively. The average node degree (i.e., the average number of connections) is 1.43529412. The network diameter is 7, and the average distance between the nodes is 2.26761. There are 59 orthologous enzymes, 16 compounds, eight pathways, and two unknown groups in this metabolic network that have no references in KEGG. The numbers of activation reactions, expression reactions, and inhibition reactions are 29, 26, and 6, respectively.

The core network of the HIF-1 signaling pathway is identified based on cohesive subgroup analysis. The network is extracted for the convenience of further study, as shown in Fig. 3. The quantitative indicators of the network are as follows: the number of nodes is 34, and the number of lines is 34. The network has a ring density of 0.02941176 and a ringless density of 0.03030303. The average node degree is 2. The network diameter is 7, and the average distance between the nodes is 2.57018. There are 24 orthologous enzymes and 10 compounds in this metabolic network. The number of unknown groups that have no references in KEGG is 1; the number of activation reactions is 27; the number of expression reactions is 2; and the number of inhibition reactions is 5.

Centrality analysis

Cohesive subgroup analysis

First, the distribution of the subset of the overall network is determined. Then, the core network based on the distribution is determined. The minimum component size is 1; i.e., the smallest subnet can be an isolated node without connections. Table 4 shows the analysis results.

The data in the table show that the overall network of the HIF-1 signaling pathway includes 25 subnets. The minimum size of the subnet is 1; i.e., there is only one node in the subnet. The maximum size is 34; i.e., 34 nodes are in the subnet.

We colored different subnets differently for the convenience of analysis in this paper. The number of the subnet to which the node belongs is identified. The subnet distribution map of the overall network of the HIF-1 signaling pathway is obtained based on Pajek, as shown in Fig. 2. Figure 2 depicts the subnet distribution of the overall network of the HIF-1 signaling pathway graphically. For example, “[6.00] 20” means that node 20 belongs to the number 6 subnet.

Table 4 and Fig. 2 show that the number 6 subnet contains most nodes of all the subnets and has meaningful reaction relationships. Therefore, subnet 6 is the core network of the HIF-1 signaling pathway.

Bow-tie structure analysis

The Bow-tie analysis module of Pajek software was used to analyze the macrostructure of the core network of the HIF-1 signaling pathway. Table 5 shows the distribution of the structure.

We colored the different structures in the core network topology of the HIF-1 signaling pathway, as shown in Fig. 3.

The specification of the vertices of the core component of the HIF-1 signaling pathway is shown in Table 6.

Structural complexity analysis

We analyzed the distribution of network node degrees. That is, the scale of the core network of the HIF-1 signaling pathway is measured. The degree distribution of the core network of the HIF-1 signaling pathway is shown in Table 7.

Based on the above methods, the power-law fitting results of the core network scale of the HIF-1 signaling pathway in this study are shown in Fig. 4.

Based on Clauset’s method, the power-law fitting model for the core network of the HIF-1 signaling pathway can be calculated as follows: (7)${\displaystyle P_r\left[X\ge x\right]\sim x^{-3.28}}$

The scale parameter of the power-law distribution of the core network of the HIF-1 signaling pathway is 3.28, and x_min = 2. That is, the nodes with a degree of 2 or more in the network obey the Pareto distribution rate. The similarity detection index of this fitted model is much larger than 0.1 and closer to 1. This indicates that the fitting effect is good.