Pathway crosstalk perturbation network modeling for identification of connectivity changes induced by diabetic neuropathy and pioglitazone

The pathway Crosstalk perturbation network model

In this work, the pathway crosstalk perturbation network (PXPN) is proposed as a model for integrating high-throughput perturbation biological data to gain insights into the changes in communication between functional biological processes. Figure 1 illustrates a schematic representation of the elements in the model, while the formal definitions of the elements in the PXPN model are provided in Additional file 1. The model is agnostic to the type of high-throughput data, the pathway description, and the statistical measure or algorithm used to define enrichment.

Basically, the PXPN model consists of four steps: 1) identification of perturbed pathways between two physiological states, 2) identification of crosstalk among the perturbed pathways, 3) identification of perturbations in the crosstalk regions between perturbed pathways, and 4) Network integration. A pseudocode representation of this model is available in Additional file 2. The scripts used in the current study are available in our Github repository (https://www.github.com/hurlab/pxpn_neuropathy).

Step 1 involves taking i) an expression dataset with information on two physiological conditions and ii) a list of pathways, defined by a single set of inclusion criteria and curatorial rules, such as those obtained from the same database. These are used as input for an enrichment function to obtain a list of pathways that are considered to be “enriched,” which in this context indicates a perturbation in pathway activity between the physiological states.

Step 2 involves looking for crosstalk between the pathways that were identified as perturbed in the previous step. Crosstalk between pathways is found if the pathways share genes, or, in other words, if the lists of genes for two pathways overlap (see definition in Additional file 1). The intersection between these two lists represents the genes that belong to the crosstalk region (or regions) between these pathways. Importantly, our research focused only on the crosstalk regions of pathways that were identified as perturbed in step 1.

Step 3 involves taking the list of crosstalk regions previously identified and using the expression data and the same enrichment function on it. Doing this allows identifying which crosstalk regions are perturbed themselves. Perturbation in the crosstalk region between two pathways indicates a change of activity in the molecules that are shared between pathways observed between the two physiological conditions, which in turn indicates a perturbation in the communication between pathways, as the expression of genes that connect the two pathways is being collectively perturbed.

Step 4 involves integrating the outputs of steps 1 and 3 into a network model. This is done by representing the perturbed pathways from step 1 as nodes in a graph, then establishing undirected links between pairs of pathways if the crosstalk region between them was identified as perturbed in step 3. The resulting undirected graph is referred to as a Pathway Crosstalk Perturbation Network, which represents the pathways that are perturbed between two biological states, along with the perturbations found in their crosstalk regions. This network model can be further analyzed from a graph-theoretic perspective.

The present research focused on the changes in pathway perturbation and communication through the analysis of changes in topology, modular structure, and connectivity, in PXPNs associated to physiological transitions. Given two phenotypes, such that one can give way to the other sequentially, the transition from the first to the second phenotype may involve the perturbation of a set of biological functions, which can be modeled using a PXPN. The progression from a physiologically healthy state to a pathological state of disease is an important biomedical case. This pathological state may, in turn, by using pharmacological agents, advance to a state of partially restored functionality. Hypothetically, a “perfect” drug could induce a final transition from the pharmacological state of partial functionality back to a healthy state indiscernible from the original physiological state. Each of these transitions can be modeled as three different PXPNs that represent the perturbations associated with each transition. As a case study, this model was implemented with data from a study of the effects of pioglitazone on murine T2DM-associated neuropathy.

RNA-Seq data

RNA-Seq raw data were obtained from a previous study on pioglitazone’s effects on diabetic complications [12] using leptin receptor-deficient db/db mice, a model of T2DM. In brief, male C57BLKS (BKS) db/+ and db/db mice (BKS.Cg-m+/+Leprdb/J) were purchased from the Jackson Laboratory (Bar Harbor, ME). Mice were fed with or without 15 mg/kg pioglitazone (112.5 mg pioglitazone/kg chow for a dose of 15 mg/kg to the mouse) for 11 weeks starting from 5 weeks of age. Pioglitazone treatment normalized renal function and improved small nerve functions but did not improve large fiber functions. Four complication-prone tissues—sciatic nerve (SCN), dorsal root ganglia (DRG), glomeruli, and kidney cortex—were collected at 16 weeks of age and examined for their genome-wide gene expression profiles using RNA-Seq (HiSeq 2000 paired-end read length of 100 bases). The current study focused on the three groups of SCN data, including db/+ (nondiabetic denoted as “healthy” group), db/db (diabetic denoted as “disease” group), and db/db + Pio (diabetic with pioglitazone treatment denoted as “treatment” group). There were n = 6 samples in each group.

The raw sequencing reads were first cleaned by removing reads containing the adapter or poly-N and removing low-quality (quality score < 30) reads using Trimmomatic [14]. FastQC version 0.10.1 [15] was used to assess the quality of raw reads. Clean reads were mapped to the mouse reference genome GRCm38 (mm10) using Hisat2 [16]. The mapping summaries—such as the percentage of reads that were uniquely mapped, multiple mapped, or unmapped—were then collected from the log files of Hisat2 runs. FeatureCounts [17] was used to count reads mapped to individual genes. Only uniquely mapped reads were used in the counting step. Then, the counting metrics were collected from the summary file of each FeatureCounts run. Genes were omitted with zero expression across all samples for the correlation and differential expression analysis. Fragments per kilobase of exon per million mapped reads (FPKM) as a measurement of transcript expression were calculated using an in-house script.

Pathway enrichment algorithm

The Reactome [18] collection of pathways was used in this study. We used the complete set of Reactome murine pathways available through the Graphite R/Bioconductor package [19]. The generally applicable gene set enrichment (GAGE) [20], a cutoff-free enrichment algorithm, was used to identify significantly enriched pathways that were perturbed by diabetes or treatment. The algorithm was run considering undirected perturbations, with an enrichment significance threshold set to q-value < 0.05.

Network analysis

Calculations for topological parameters—degree, clustering coefficient (CC), network density, average path length, and number of connected components (islands in the network)—were carried out using Igraph [21] for R, NetworkX [22] for Python, and Cytoscape 3.3.0 [23]. Additionally, communities (subsets of nodes with high intraconnectivity and low outbound connections) were detected using the Infomap algorithm [24], as implemented in the Igraph package.

Implementation for the study of physiological transitions in the murine diabetic neuropathy model

For this study of murine diabetic neuropathy, the previously described RNA-Seq expression dataset and pathways from the Reactome database were the inputs of the model. The differences between these groups represent the transitions observed in a patient. First, the patient transitions from a physiologically functional state to a pathological state (health to disease, denoted as HTD). Given therapy, the patient transitions from the pathological state to a pharmacologically modulated state (disease to treatment, denoted as DTT). Finally, with successful therapy the patient transitions back to the physiological state (treatment to health, denoted as TTH). Three networks, each representing one of these physiological transitions, were constructed. It is proposed that changes in pathway connectivity in different transitions indicate changes in the overall impact of a particular pathway activity in the observable phenotype.

Null model

To evaluate the significance of the topological parameters of these three PXPNs, an ensemble of 5000 networks was generated for each transition using a null model by randomly rewiring the edges, with a rewiring probability proportional to the size of the intersection between two pathways (measured as the Jaccard index). For each network, each topological parameter was compared against the null model using a Z-test. This model allowed assessing whether particular topological properties of the obtained network differed from a randomly generated network containing the same number of nodes and edges (not all edges are possible, as not all pathways crosstalk with each other).

A second null model, for evaluating the overall capacity of the method of obtaining non-trivial network structures from gene expression measurements, was employed. For each of the three comparisons (HTD, DTT, and TTH), an ensemble of 1500 random expression datasets was generated by shuffling the gene labels of the original RNA-seq data. These data were used to generate PXPNs using the established pipeline and compared against those of the comparisons.

Network overview

Using the proposed approach, each transition between physiological states was represented as a network with characteristic structural features. The generated networks may be found in Additional files 3, 4, and 5. Figure 2 illustrates the HTD network, which represents pathway alterations associated with the transition from the healthy state to the pathological neuropathic condition. Between these two states, 104 pathways were altered, and 222 significant alterations to the activity of crosstalk regions were observed. The second network, visualized in Fig. 3, represents the transition from the pathological state to a pharmacologically modulated state (the DTT network). This transition was associated with 78 altered pathways and 149 crosstalk perturbations among them. Finally, the TTH network, as illustrated in Fig. 4, describes the alterations found between the pharmacologically treated state and the healthy state, which would represent the perturbations observed in the transition back to the healthy state. This TTH network included 110 altered pathways, with 213 edges representing perturbed crosstalks between them. Additional file 6 illustrates the overlap among the pathways in these three networks; these pathways and their q-values for each transition are listed in Additional file 7.

Each physiological transition involves a precise combination of perturbed pathways, with specific communication patterns between them. This is evident in terms of similarity among the networks as each network has a portion of exclusive and shared nodes. In all three transitions, 53 pathways were perturbed and therefore represented as nodes in the network. The most similar networks in terms of both nodes and edges were the DTT and TTH networks (node Jaccard index = 0.55; edge Jaccard index = 0.45); the most different were the HTD and TTH networks (node Jaccard index = 0.49; edge Jaccard index = 0.43). Additional files 8 and 9 contain similar values among the three networks in terms of nodes and edges, respectively.

As each transition was associated with a specific network, each network was associated with specific structural characteristics. The communication between pathways associated with each transition between physiological states was different, resulting in a unique network topology. This in turn was reflected in network properties, such as the average path length, the clustering coefficient, and the distribution of nodes into connected components and communities inside connected components.

Global network parameters

The three experimentally derived networks had clustering coefficient values that were higher than the ones expected from the null model (for instance, the average clustering coefficient value for the HTD null model networks was 0.149). Nevertheless, the clustering coefficient values of all networks were comparable (ranging from 0.513 to 0.620). Average path lengths were slightly higher in the experimentally derived networks than in the null model networks in the cases of the HTD and TTH transitions. Interestingly, in the case of the DTT transition, the average path length was considerably lower than the one predicted by the null model (1.99, compared with the predicted 3.26). This result suggests that the transition induced by pharmacological treatment on the pathological state involves perturbations with short-range pathway communication, whereas the perturbations from and to the physiological state require longer-range changes in communication.

All three networks had low density for edges. Since not all edges were biologically possible, as not all pathways are able to crosstalk, comparing the number of edges against the total possible crosstalks between the pathways in each network is important; we refer to this as the specific density of a PXPN. In the HTD network, 19% of possible crosstalks were perturbed, where 26 and 20% of possible crosstalks were perturbed in the DTT and TTH networks, respectively. Derived from this is the observation of pathways that while they have a large crosstalk potential with other perturbed pathways in a transition, they appear disconnected nonetheless. For instance, in the HTD network, the “regulation of insulin secretion” could potentially connect to 32 pathways, but it was disconnected, indicating that the crosstalk of this pathway was not altered during the transition from health to disease. This lack of observed connectivity indicates that, in this transition, at least at the level of gene expression perturbation, this pathway has little system-wide influence.

Connectivity and modular structure of networks

Our proposed model allowed the representation of the alterations between physiological states in pathway activity and communication as a graph. The structures of these networks were nontrivial and different from those of the random networks because the connections reflect the perturbation of crosstalk regions associated to each physiological transition. Therefore, there are differences with respect to the connected components (subgraphs in which any pair of nodes has a path between them) and communities (modules inside a connected component in which the nodes belonging to the same module have a higher number of edges between them than the nodes outside of the module). Differences in the organization of these networks are indicative of specific communications between biological processes that were altered in each physiological transition.

The HTD network (Fig. 2) was composed of 22 connected components, 14 of which were single nodes. The largest connected component contained 66 pathways (~ 63% of all 104 pathways), which were related to the “neuronal system” processes. This component also included other smaller communities related to “solute carrier (SLC)-mediated transmembrane transport,” “extracellular matrix (ECM) organization,” “G-protein coupled receptor (GPCR) ligand binding,” and a community containing diverse pathways, such as “hemostasis” and “GPCR downstream signaling.” The second largest connected component contained 10 pathways (~ 10% of all nodes), which were mainly related to “lipoprotein metabolism.”

The DTT network (Fig. 3) was composed of 16 connected components, 7 of which were single nodes. In this network, the largest connected component contained only 21 pathways (~ 27% of all 78 pathways), which were similar to those found in the “neuronal system” community in the HTD network. The second largest connected component (19 pathways corresponding to ~ 24% of the network nodes) was composed of three communities: one comparable to the “lipoprotein metabolism” community in the HTD network, another related to the “metabolism of lipids,” and the other containing three pathways related to “retinoids.” Three other connected components were comparable to the communities found in the HTD network, such as “ECM organization” (8.97%), “GPCR ligand binding” (8.97%), and “SLC-mediated transmembrane transport” (6.41%). Another connected component with four pathways related to the TCA cycle was also found.

The TTH network (Fig. 4) was composed of 19 connected components, 11 of which were single nodes. This network was dominated by the largest connected component, which contained 59 pathways (~ 54% of all 110 pathways). The communities in this component were similar to those in the largest component of the HTD network, including the “neuronal system,” “SLC-mediated transmembrane transport,” “ECM organization,” “GPCR-ligand binding,” and “GPCR-downstream signaling.” Interestingly, a new community emerged, containing the “glycolysis,” “gluconeogenesis,” and “glucose metabolism” pathways, which were connected to the “ECM organization” community. A notable difference between this network and the HTD network was the changes in the communication of the GPCR community, which became disconnected from the “ECM organization” community and connected to the “neuronal” community.