Exploring Wound-Healing Genomic Machinery with a Network-Based Approach

The molecular mechanisms underlying tissue regeneration and wound healing are still poorly understood despite their importance. In this paper we develop a bioinformatics approach, combining biology and network theory to drive experiments for better understanding the genetic underpinnings of wound healing mechanisms and for selecting potential drug targets. We start by selecting literature-relevant genes in murine wound healing, and inferring from them a Protein-Protein Interaction (PPI) network. Then, we analyze the network to rank wound healing-related genes according to their topological properties. Lastly, we perform a procedure for in-silico simulation of a treatment action in a biological pathway. The findings obtained by applying the developed pipeline, including gene expression analysis, confirms how a network-based bioinformatics method is able to prioritize candidate genes for in vitro analysis, thus speeding up the understanding of molecular mechanisms and supporting the discovery of potential drug targets.


Introduction
In recent years, the process of tissue-regeneration has been widely studied to discover treatments for controlling the progression of lesions and for restoring diseased or injured tissue functions. Understanding how damaged, defective (e.g., cancer [1]), or lost tissues can be regenerated is one of the major challenges in biomedical research and in regenerative medicine. The wound healing process plays a key role in tissue regeneration; therefore, a better understanding of the molecules involved in this process, and how they interact, provides new insights into the development of drug treatments.
The healing process is dependent on interactions between many cell types, extracellular matrix, growth factors, and mediators in a specific temporal phase [2]. Schematically, the healing process

Results
Literature analysis and knowledge extracted from public repositories (Section 4.1) led to the selection of a small list of genes related to the wound healing process or to the Rigenera ® action. We will refer to this list as the input genes. The literature search (performed according to Section 4.1) allowed us to extract a list of seven genes: Tnf, Cxcl2, Ccl12, Fgf5, Wnt5a, Col3a1, Fosb, and Pgk1; the latter has been included to expand the scope of the network. In fact, besides being a housekeeping Pharmaceuticals 2017, 10, 55 3 of 18 gene, Pgk1 is involved in epidermal regeneration after certain stimuli [15,16]. The input genes and the related literature findings are reported in Table 1. Table 1. Input genes set.

Gene Symbol
Uniprot ID Protein Names Ref.

Assessing the Starting Node Set
We verified how the input genes are involved in wound healing and enhanced by Rigenera ® via measuring gene expression after a skin scratch test [29] (methods in Section 4.1.1). Figure 1 shows how all the nodes, with the exception of the untested housekeeping Pgk1, are stimulated by the Rigenera ® action. Moreover, all the genes, except for Tnf, are activated after the skin scratch test. These seven genes are important in the three phases of wound healing; inflammation, proliferation, and remodeling. Measured gene expressions are shown in Figure 1. At five hours, the input gene expressions were significantly different between no Scratch-no Regenera vs. yes Scratch-yes Regenera when taken as a system of seven genes (Z Stouffer = 0.04-Equation (1)), as well as for three individual gene expressions (Fgf5, p-value = 0.04; Wnt5a, p-value = 0.03; Col3a1, p-value = 0.04; t-test unequal variance-Section 4.1.1.). As expected at 1 h, the input genes were also significantly dysregulated in three genes (no Scratch-yes Rigenera, Tnf, p-value = 0.03; no Scratch-yes Rigenera, Cxcl2, p-value = 0.02; no Scratch-yes Rigenera, Fosb, p-value = 0.015; t-test unequal variance-Section 4.1.1.).

Assessing the Starting Node Set
We verified how the input genes are involved in wound healing and enhanced by Rigenera ® via measuring gene expression after a skin scratch test [29] (methods in Section 4.1.1). Figure 1 shows how all the nodes, with the exception of the untested housekeeping Pgk1, are stimulated by the Rigenera ® action. Moreover, all the genes, except for Tnf, are activated after the skin scratch test. These seven genes are important in the three phases of wound healing; inflammation, proliferation, and remodeling. Measured gene expressions are shown in Figure 1. At five hours, the input gene expressions were significantly different between no Scratch-no Regenera vs. yes Scratch-yes Regenera when taken as a system of seven genes (ZStouffer = 0.04-Equation (1)), as well as for three individual gene expressions (Fgf5, p-value = 0.04; Wnt5a, p-value = 0.03; Col3a1, p-value = 0.04; t-test unequal variance-Section 4.1.1.). As expected at 1 h, the input genes were also significantly dysregulated in three genes (no Scratch-yes Rigenera, Tnf, p-value = 0.03; no Scratch-yes Rigenera, Cxcl2, p-value = 0.02; no Scratch-yes Rigenera, Fosb, p-value = 0.015; t-test unequal variance-Section 4.1.1.). We set up a bioinformatics pipeline to infer other genes potentially involved in the wound healing process, a so called candidate gene set, based on (a) the starting genes and (b) on the topological knowledge extracted from a protein network. An outline of the presented method is depicted in Figure 2. We set up a bioinformatics pipeline to infer other genes potentially involved in the wound healing process, a so called candidate gene set, based on (a) the starting genes and (b) on the topological knowledge extracted from a protein network. An outline of the presented method is depicted in Figure 2.
We set up a bioinformatics pipeline to infer other genes potentially involved in the wound healing process, a so called candidate gene set, based on (a) the starting genes and (b) on the topological knowledge extracted from a protein network. An outline of the presented method is depicted in Figure 2.  Overview of the proposed computational systems biology method. The network construction is performed through the selection of input genes (shown with red diamonds) involved in the wound healing process. The knowledge about Protein-Protein Interaction (PPI) (STRING repository [30]) is then used to build the network. Next, the topological network analysis selects candidate nodes (i.e., cluster hubs and bridge nodes) and identifies other key genes in the wound healing process. Finally, the simulation of the Rigenera ® autologous micro-graft action is performed in the TNF signaling pathway from the Kyoto Encyclopedia of Genes and Genomes (KEGG) [31] through Boolean networks.

A Protein Network for Wound Healing
We built the wound healing network using the input genes and the STRING repository [30] (Section 4.2). The obtained network, depicted in Figure 3, consists of 446 nodes and 24,757 edges.
To prove the network's meaningfulness, we performed a KEGG and Gene Ontology (GO) [32] enrichment analysis (Section 4.2.1), finding 62 pathways and 40 GO terms significantly associated with the network (threshold 10 −3 , hypergeometric test), of which 13 pathways and 19 GO terms specifically pertain to wound healing. The complete list is reported in Supplementary Tables S1 and S2.
repository [30]) is then used to build the network. Next, the topological network analysis selects candidate nodes (i.e., cluster hubs and bridge nodes) and identifies other key genes in the wound healing process. Finally, the simulation of the Rigenera ® autologous micro-graft action is performed in the TNF signaling pathway from the Kyoto Encyclopedia of Genes and Genomes (KEGG) [31] through Boolean networks.

A Protein Network for Wound Healing
We built the wound healing network using the input genes and the STRING repository [30] (Section 4.2). The obtained network, depicted in Figure 3, consists of 446 nodes and 24,757 edges.
To prove the network's meaningfulness, we performed a KEGG and Gene Ontology (GO) [32] enrichment analysis (Section 4.2.1), finding 62 pathways and 40 GO terms significantly associated with the network (threshold 10 −3 , hypergeometric test), of which 13 pathways and 19 GO terms specifically pertain to wound healing. The complete list is reported in Supplementary Tables S1 and S2.
(a) (b) legend) refer to the nodes that, according to ClusterONE, did not end up in any of the clusters due to their topological properties [33]. The hub nodes of each cluster are shaped with triangles.
In the following, we will focus on (a) cluster hubs to study the wound healing machinery with bioinformatics, and (b) bridge nodes to unveil their potential as drug targets with in vitro experiments.

Pivotal Genes in Wound Healing: Hub Nodes of Clusters as the Key Elements in Regulation
Cluster hubs are pivotal nodes in highly connected modules (i.e., clusters) and they have a highly significant number of neighbors. Therefore, these nodes are intrinsically important for understanding the underlying genetic mechanisms.
To identify the hub nodes within a cluster (i.e., cluster hubs), we first partitioned the network with ClusterONE [33]. Of the eleven obtained clusters (Table 2), we retained seven showing a p-value < 0.05 (one-sided Mann-Whitney U test). Next, for each cluster, hub nodes were ordered according to their degree (i.e., the number of direct node neighbors), and only the ones ranking beyond the 90th percentile (k > 203.5) were retained (Section 4.3.1) as a candidate gene set (cluster hubs). legend) refer to the nodes that, according to ClusterONE, did not end up in any of the clusters due to their topological properties [33]. The hub nodes of each cluster are shaped with triangles.
In the following, we will focus on (a) cluster hubs to study the wound healing machinery with bioinformatics, and (b) bridge nodes to unveil their potential as drug targets with in vitro experiments.

Pivotal Genes in Wound Healing: Hub Nodes of Clusters as the Key Elements in Regulation
Cluster hubs are pivotal nodes in highly connected modules (i.e., clusters) and they have a highly significant number of neighbors. Therefore, these nodes are intrinsically important for understanding the underlying genetic mechanisms.
To identify the hub nodes within a cluster (i.e., cluster hubs), we first partitioned the network with ClusterONE [33]. Of the eleven obtained clusters (Table 2), we retained seven showing a p-value < 0.05 (one-sided Mann-Whitney U test). Next, for each cluster, hub nodes were ordered according to their degree (i.e., the number of direct node neighbors), and only the ones ranking beyond the 90th percentile (k > 203.5) were retained (Section 4.3.1) as a candidate gene set (cluster hubs). Internal weight refers to the total weight of the edges contained in the cluster; External weight denotes the total weight of the edges that connect the cluster nodes with the rest of the network [33].
Cluster nodes were enriched with KEGG pathways and GO terms. In particular, clusters were characterized by dozens of statistically significant GO terms and KEGG pathways (p-value < 0.001, hypergeometric test), as summarized in Table 3 (detailed results are reported in Supplementary Tables S3-S16).

Pivotal Genes in Wound Healing: Bridge Nodes
Bridge nodes are measures of the betweenness of clusters. They are sparsely connected and located among clusters. Identifying bridge proteins is generally done to unveil viable targets for pharmaceuticals [34,35]. By applying the bridging procedure (Section 4.3.2), we identified 45 bridge nodes. Next, we calculated the S tot score (Equation (4)) of each bridge node by considering the wound healing-related KEGG pathways and/or GO terms (Supplementary Tables S17 and S18). The final ranking is reported in Supplementary Table S19, and it represents the estimated likeliness of Pharmaceuticals 2017, 10, 55 7 of 18 a bridge node being pivotal in wound healing. Based on the resulting prioritization, we conducted gene expression analysis of the top three scored genes, namely Nfkb1, Rela, and Tnfrf1a (applying the procedure in Section 4.1.1). The experimental results reported in Figure 4 support the findings predicted by the topology analyses. In fact, at five hours, a trend was observed between no Scratch-no Regenera vs. yes Scratch-yes Regenera for the three output genes (bridge nodes) (Z Stouffer = 0.105-Equation (1); Tnfrf1a, p-value = 0.086; Rela, p-value = 0.22; Nfkb1, p-value = 0.5; t-test unequal variance-Section 4.1.1.).
Bridge nodes are measures of the betweenness of clusters. They are sparsely connected and located among clusters. Identifying bridge proteins is generally done to unveil viable targets for pharmaceuticals [34,35]. By applying the bridging procedure (Section 4.3.2), we identified 45 bridge nodes. Next, we calculated the Stot score (Equation (4)) of each bridge node by considering the wound healing-related KEGG pathways and/or GO terms (Supplementary Tables S17 and S18). The final ranking is reported in Supplementary Table S19, and it represents the estimated likeliness of a bridge node being pivotal in wound healing. Based on the resulting prioritization, we conducted gene expression analysis of the top three scored genes, namely Nfkb1, Rela, and Tnfrf1a (applying the procedure in Section 4.1.1). The experimental results reported in Figure 4 support the findings predicted by the topology analyses. In fact, at five hours, a trend was observed between no Scratchno Regenera vs. yes Scratch-yes Regenera for the three output genes (bridge nodes) (ZStouffer = 0.105-

Simulation of KEGG TNF Signaling Pathway Dynamics
Following the suggested pipeline, we modeled the KEGG TNF signaling pathway dynamics with a Boolean network (BN) approach (Section 4.4). The resulting network is illustrated in Figure 5 and consists of 44 nodes and 49 edges.

Simulation of KEGG TNF Signaling Pathway Dynamics
Following the suggested pipeline, we modeled the KEGG TNF signaling pathway dynamics with a Boolean network (BN) approach (Section 4.4). The resulting network is illustrated in Figure 5 and consists of 44 nodes and 49 edges.
Next, we simulated the Rigenera ® autologous micro-graft action in the BN by initializing the state (upregulated) of its target genes whose expression is known, i.e., the input nodes ( Figure 1) and the candidate gene set (bridges) (Figure 3) overlapping with the TNF pathway. These nodes, Tnf, Rela and Nfkb1, are all set to 1 and are highlighted in red in Figure 5. On the other hand, initial values of the other network nodes are randomly assigned to 0 or 1 (Section 4.4). We next applied Odefy, which performs the conversion of BNs into continuous ODEs, with continuous time descriptions, and gives as an output the node behaviors after a certain stimulus (here Rigenera ® action). In order to be independent from the random assignment of non-target nodes, we applied Odefy with 1000 Monte Carlo simulations, obtaining 1000 behaviors of all network nodes (one example is shown in Figure 6) after the Rigenera ® stimulus. The analysis of the results over the simulations allowed us to rank the genes according to the frequency of their state changes. In detail, we looked for the genes that consistently switch their state over the simulations, indicating they have been stimulated by Rigenera ® . The ranked list is showed in Supplementary Table S20 and  Next, we simulated the Rigenera ® autologous micro-graft action in the BN by initializing the state (upregulated) of its target genes whose expression is known, i.e., the input nodes ( Figure 1) and the candidate gene set (bridges) ( Figure 3) overlapping with the TNF pathway. These nodes, Tnf, Rela and Nfkb1, are all set to 1 and are highlighted in red in Figure 5. On the other hand, initial values of the other network nodes are randomly assigned to 0 or 1 (Section 4.4). We next applied Odefy, which performs the conversion of BNs into continuous ODEs, with continuous time descriptions, and gives as an output the node behaviors after a certain stimulus (here Rigenera ® action). In order to be independent from the random assignment of non-target nodes, we applied Odefy with 1000 Monte Carlo simulations, obtaining 1000 behaviors of all network nodes (one example is shown in Figure 6) after the Rigenera ® stimulus. The analysis of the results over the simulations allowed us to rank the genes according to the frequency of their state changes. In detail, we looked for the genes that consistently switch their state over the simulations, indicating they have been stimulated by Rigenera ® . The ranked list is showed in Supplementary Table S20    Next, we simulated the Rigenera ® autologous micro-graft action in the BN by initializing the state (upregulated) of its target genes whose expression is known, i.e., the input nodes ( Figure 1) and the candidate gene set (bridges) (Figure 3) overlapping with the TNF pathway. These nodes, Tnf, Rela and Nfkb1, are all set to 1 and are highlighted in red in Figure 5. On the other hand, initial values of the other network nodes are randomly assigned to 0 or 1 (Section 4.4). We next applied Odefy, which performs the conversion of BNs into continuous ODEs, with continuous time descriptions, and gives as an output the node behaviors after a certain stimulus (here Rigenera ® action). In order to be independent from the random assignment of non-target nodes, we applied Odefy with 1000 Monte Carlo simulations, obtaining 1000 behaviors of all network nodes (one example is shown in Figure 6) after the Rigenera ® stimulus. The analysis of the results over the simulations allowed us to rank the genes according to the frequency of their state changes. In detail, we looked for the genes that consistently switch their state over the simulations, indicating they have been stimulated by Rigenera ® . The ranked list is showed in Supplementary  The top 10 gene list was analyzed by MouseMine [36]. In Table 4, the GO terms related to wound healing and tissue regeneration are reported.  The top 10 gene list was analyzed by MouseMine [36]. In Table 4, the GO terms related to wound healing and tissue regeneration are reported. Table 4. GO biological processes associated with the top 10 genes stimulated by Rigenera ® .

Automatically Assembling a Gene Network
Our pipeline constructs a PPI network by expanding a list of wound healing-related genes ( Table 1). An initial assessment of gene expression levels confirms the input gene list as pertaining to wound healing or to the action of Rigenera ® (Figure 1). In fact, all the genes respond positively to the skin scratch test. Tnf, Cxcl2, Ccl12, and Fosb are upregulated after the administration of Rigenera ® even in the absence of the scratch. When Rigenera ® is administered along with the skin scratch test, all the genes show an upregulation, mostly higher than all the other experimental settings, especially after five hours (statistical significance reported in Section 2.1).
After network enrichment, p-values of wound healing GO terms and KEGG pathways show how the network well represents the wound recovery process (Supplementary Tables S1 and S2). In the significantly enriched lists, different wound healing-related concepts emerge. Within them there are inflammation (GO:0006954) and angiogenesis (GO:0001525), along with their related terms. The inflammatory response, a well-known phase of the wound healing process, besides being present as a GO term, is underlined by the presence of cytokine-related [37] terms (GO:0016493, GO:0008009, GO:0006935) and pathways (mmu04062, mmu04060, mmu04750). The network is also enriched by the term extracellular matrix (GO:0031012) reflecting the pertinence of its nodes to the tissue regeneration.
The Wnt pathway is activated for cell fate specification, cell proliferation, and cell migration. Wnt is also an important signaling pathway during skin wound healing [38]. Studies show how the expression of Wnt-5a in the wound promotes differences mimicking regeneration, including the development of epithelia-lined cysts in the dermis, hair follicles, and sebaceous glands, without the development of tumors [39]. Notably, an excess of Wnt signaling pathway concepts stands out among the significantly associated terms (GO:0060070, GO:0035567, GO:0090263, GO:0016055, GO:0060071, GO:0042813, GO:0017147) and pathways (mmu04310).
Once the network is assembled, we identified pivotal cluster hubs and bridge nodes. The first are important as actors to understand the wound healing genetic orchestration, and the latter as candidate drug targets.

Unveiling Cluster Hub Nodes
Functional regions, or clusters, allow the segregation of the network in functional modules. These modules are consistently characterized by specific processes (Table 3), reflecting how topologically-differentiated areas of the network are specialized in distinct aspects of wound healing.

•
Cluster 1 gathers chemokine/chemotaxis processes, related to cell migration towards chemical gradients. In particular, chemotaxis plays an important role during the inflammatory phase of healing processes [40]. • Cluster 2 is labeled by extracellular matrix/skin development terms. • Cluster 3 gathers genes involved mainly in DNA regulation of transcription. • Cluster 4 gathers regenerative processes and cell growth, being labelled by cancer, pluripotency, Wnt, and mTOR signaling pathways. • Cluster 5, like cluster 2, is characterized by processes of cell adhesion and cell-cell junction formation. This is confirmed by the presence of the Rap1 and Ras signaling pathways, both involved in cell proliferation, survival, growth, migration, differentiation, or cytoskeletal dynamism. • Cluster 6 is characterized by inflammation and immune response processes. In fact, it is enriched by terms such as MAPKs, TNF signaling pathway, inflammation regulation, and leukocyte migration. The presence of osteoclast proliferation calcium ion terms is consistent with the involvement of this cluster in new bone formation. Although not directly relevant in wound healing, the fact that bone formation terms were grouped in the same cluster indicates how the clustering technique successfully gathered similar genes in a consistent, meaningful fashion. • Cluster 8 is labeled by terms mostly related with Glycine and sugar metabolism. Pgk1 appears in this cluster, the housekeeping gene involved in glycolysis [15,16].
Within the clusters, cluster hub genes are identified by network theory topological rules. These hubs are genes reflecting the general functions of their respective clusters, and are associated to consistent specific GO terms and KEGG pathways (see Table 3).

Unveiling Bridge Nodes
Bridge nodes have been found to be associated to wound healing-related processes (Supplementary Tables S17 and S18) other than the ones identified for the whole network (Supplementary Tables S1 and S2). This may be due to their special property as individual actors in bridging the different functions related to the various gene clusters. For example, bridge nodes are related to the vascular endothelial growth factor (VEGF) pathway (mmu04370). VEGF plays an important role during the wound healing process, since it is known to be an inducer of cell migration through chemotaxis. It is considered one of the most important factors for its effects on the healing process, including angiogenesis and collagen formation. Different bridge nodes are also related to the apoptotic process (mmu04210). Apoptotic genes play an important role in different physiological processes by removing damaged and potentially dangerous cells. This event is known as apoptosis-induced proliferation or compensatory proliferation and it is important in tissue regeneration [41]. A recent study also shows the involvement of the JAK-STAT signaling pathway (present in the list as mmu04630) in the promotion of compensatory proliferation. This pathway is the principal signaling mechanism for a wide array of cytokines and growth factors [42].
Three genes, Rela, Nfkb1, and Tnfrsf1a, are identified as main bridges. Experiments on gene expression show how their response, with or without the Rigenera ® stimulus, is dependent on a skin scratch test. Nfkb1 (Figure 4, panel (a)) responds to the scratch test by increasing its expression level, peaking after five hours. Its expression is also stimulated by Rigenera ® in the absence of a skin scratch test. The highest fold change is however reached only when Rigenera ® is administered after the scratch test. Rela (Figure 4, panel (b)) shows a limited increase in its expression after scratch, or after the sole administration of Rigenera ® . It shows its highest levels of expression when Rigenera ® is administered after a skin scratch test. Finally, the expression of Tnfrsf1a (Figure 4, panel (c)) does not seem to be activated by the scratch test; however, it seems to react to the administration of Rigenera ® after five hours, without or with the scratch tests. The analysis of the biological results, confirmed that these genes tend to be related to the action of Rigenera ® on a skin scratch (Section 2.4). The biological experiments with n = 3 in four groups (no Scratch-no Rigenera ® , no Scratch-yes Rigenera ® , yes Scratch-no Rigenera ® , yes Scratch-yes Rigenera ® ) were underpowered to evaluate all combinations and control for statistical multiplicity. Therefore, our method suggests three targets, among all possible candidates, for drug discovery or repositioning that may be potentially involved in the wound healing process.

A Boolean Network to Study the TNF Signaling Pathway
Finally, with the modeling of the TNF signaling pathway into a Boolean network we demonstrate the possibility of simulating the action of a treatment (Rigenera ® ) in order to provide qualitative behavior of the genes in a biological pathway. The model pathway TNF has been chosen because of its involvement in chemokine and cytokine pathways and its consequent role in the inflammation phase of wound healing [17]. The analysis of the results obtained by applying Odefy enabled the identification of a ranked list of the pathway genes based on the effect of the simulated Rigenera ® action. The top 10 genes significantly altered by the Rigenera ® action simulation are Bag4, Pik3r1, Pik3cb, Map2k6, Map3k7, Mapk10, Mapk11, Pik3ca, Map3k14, and Atf2. Their key role in wound healing was then confirmed by the GO terms related to these genes ( Table 4). The p38 MAPKs play a crucial role during the inflammatory phase and they are already known to be involved in tissue regeneration [43]. Moreover, the gene Map2k6 is the major MAPK11 activator in response to the presence of cytokines, such as TNF-α and IL-1. The Activating transcription factor 2 (Atf2) is a member of the leucine zipper family of DNA-binding proteins. Atf2 is able to respond to stress and may influence cell proliferation, inflammation, apoptosis, oncogenesis, neurological development and function, and skeletal remodeling [44]. The activation of Atf2 complexes increases the transcription of genes involved in inflammation such as cell adhesion molecules and cytokines, which are important for the recruitment of leukocytes to the site of injury. Finally, for the genes Pik3r1, Pik3cb, and Pik3ca, we found evidence in the literature (i.e., [45]) in addition to the GO terms about their relation to the wound healing process and tissue regeneration. Future studies should include conventional autologous micro-graft controls in order to identify Rigenera ® -specific networks, which are currently confounded with micro-graft action.

Summary
Our method correctly identifies Rela, Nfkb1, and Tnfrsf1a as drug targets in the wound healing process. It shows how network theories can successfully harvest biological knowledge applied to physiological processes, such as wound healing, and can unveil: (a) candidate targets for drugs/treatments; (b) pivotal genes for the process machinery; and (c) recommendations for in vitro experiments; finally, it simulates the action of a treatment in biological pathways to predict gene responses.

Materials and Methods
We developed a bioinformatics pipeline (Figure 2) consisting of four different steps: (1) Selection of the input genes; (2) Network construction; (3) Identification of candidate genes (cluster hubs and bridges); (3) Gene expression in vitro experiments; and (4) Simulation of a treatment action. The developed pipeline has been implemented in Python 2.7 and Matlab and it is available upon request from the authors. All network figures were created with Cytoscape [46].

Selection of the Input Gene Set
The first step of the proposed approach includes the selection of a set of genes involved in the process under study; in this case, wound healing and, consequently, to Rigenera ® . The gene list is generated from the literature and public repositories. To this end we looked into PubMed [47], GenBank [48], Uniprot [49], and KEGG [31], and we manually retrieved a gene list on which the network will be constructed, reported in Table 1. We obtained peer-reviewed studies by performing an automated search in PubMed by using combinations of keywords such as "wound-healing", "gene expression", "mice", and "mouse". This search yielded hundreds of publications; however, we retained only the studies whose findings were experimentally validated. Furthermore, we restricted the search to the experiments where normal vs. scratched tissue or similar approaches were applied (e.g., not considering diabetic mice or experiments involving other treatments). Then relevant studies were selected by first screening titles and abstracts, and then by analyzing full texts and experiment results. In addition, we selected as candidates the genes present on the Mouse Wound Healing RT 2 Profiler™ PCR Array (Qiagen, Valencia, CA, USA [20]). The genes identified with this procedure were then used to query public repositories to verify their relation to wound-healing. This allowed the selection of a restricted pool of genes on which to perform in vitro experiments to verify their association with Rigenera ® and wound healing. In this way, we can assess the methodology as a proof of concept to investigate even poorly understood/studied biological processes and treatments (i.e., processes where only a handful of genes is known to be influential). A detailed description of the method to select the input genes is provided in Supplementary Method S1.

Gene Expression in vitro Experiments
In order to assess if the gene set identified is related to the wound healing process, we performed gene expression assays for the selected genes.
All in vitro experiments were performed using murine cells, which were cultured from the murine model C57BL/6. By performing a scratch assay, an in vitro damage was created according to the skin scratch test standard practice [50]. Micro-grafts obtained by using the Rigenera ® protocol have been added autologously to this model. Cells were cultured in DMEM supplemented with 10% FBS and treated for 1 and 5 h at 37 • C in 5% of CO 2 .
Gene expression analysis was performed with qRT-PCR on the input genes. For this analysis, total RNA was extracted using the PureLink RNA Mini Kit (Invitrogen) and treated with the DNA-free Kit (Invitrogen). 250 ng RNA was reverse transcribed into cDNA with the Superscript III Reverse Transcriptase First-Strand Synthesis SuperMix (Invitrogen). Quantitative real-time PCR was performed with the Platinum SYBRGreen qPCR SuperMix-UDG (Invitrogen). GraphPad Prism 6 was used to graph the results of the biological replicates ± SEM. The ethical approval code is P056/2017. To assess the statistical significance of the results, we performed the t-test controlling for distinct variances to take into account the heteroscedasticity of the results [51]. Once it was assessed that the housekeeping expression values (log2 transformed) of the samples have a similar distribution across the different conditions (i.e., no Scratch-no Rigenera vs. yes Scratch-yes Rigenera), we could directly perform the t-test assuming unequal variance on log2 transformed expressions of the tested genes. Next, on the obtained p-values we performed an analysis of the biological validation results as a coordinated gene set/systems/network-level signal. In order to detect a coordinated systems/network-level among input genes, we first summarized the p-values of the input gene expression using the Stouffer Z-transform, a meta-analysis [52][53][54] that produces a joint p-value for the inputs (Equation (1)). The Z-transform test converts the t-test p-values, Px, for each of n transcripts of the gene set into a standard normal curve (Zx). The sum of these Zx's, divided by the square root of the number of tests, n, has a standard normal [54].

Protein-Protein Interaction Network Construction
The PPI network related to the wound healing is derived by integrating the knowledge on PPI from the STRING repository [30]. STRING allows the extraction of PPIs derived from different evidences, i.e., Neighborhoods, Gene Fusion events, Co-occurrence events, Co-expression data, Experimental data, Database information, and Text-mining associations. All these sources of evidence are benchmarked and calibrated against previous knowledge. Through this procedure (detailed in [55]), STRING assigns a confidence score to each known or predicted association. This score varies between 0 to 1 and it assesses the quality of the retrieved association based on the related type of evidence (low confidence: scores <0.4, medium: 0.4 to 0.7, and high: >0.7) [55]. In our approach, network nodes correspond to proteins and two proteins are connected by an edge if an interaction between them is found in STRING. Therefore, the entire set of network nodes is determined by considering the input gene set and all its neighbors in STRING (depicted with light green in Figures 2 and 3). The connectivity pattern among nodes is defined according to the STRING interactions and takes into account the confidence score. In order to restrict the edge set to the most reliable STRING associations, we introduced two constraints: (i) PPIs need to be derived from Experimental and Database evidence (considered more reliable with respect to, for example, text-mining evidence); and (ii) the STRING confidence score has to be higher than 0.7 (high-confidence association). Confidence scores have been also integrated into the network model as weights on the network edges.

Enrichment analyses
To assess the meaningful representation of the network and biologically interpret the network genes, we conducted KEGG [56] and Gene Ontology (GO) [32] enrichment analysis of the network nodes. We detected the enrichment on the basis of the hypergeometric test and ranked KEGG pathways and GO terms according to the obtained p-values corrected with a false discovery rate (FDR). Statistically significant processes and terms are finally selected by considering only p-values < 0.001.
The resulting significant pathways are then compared to reference lists of KEGG pathways and GO terms known to be related to wound healing. To this end, we manually selected 15 KEGG pathways and 47 GO terms based on their relation with wound healing (Supplementary Tables S1 and S2).

Identification of Candidate Genes
The third step of the proposed procedure aims to identify candidate genes, i.e., network nodes that can be suggested as targets for treatments or other genes that can have a key role in the process. In our approach, such genes are selected based on in-silico analysis of the constructed PPI network. Candidate genes are first selected based on their topological role in the network, and then enrichment selects the best candidates for further investigation.
Since the structure of the PPI networks reflects the biological properties of the system represented [12,13,57,58], to select candidate genes we introduced two constraints based on the topological property of the network. Candidate genes can correspond to hubs of clusters and bridge nodes.

Cluster Hubs
In biological networks, and in particular in PPI networks, the topological network structure is characterized by modules. Therefore, network nodes can be clustered into groups (clusters) and these clusters reflect specific biological functions and/or biological processes [13,58]. In graph theory, different clustering algorithms can be applied to detect network clusters [59][60][61]. Our approach uses ClusterONE [33] as it takes into account network weights (i.e., STRING confidence scores) and it allows cluster overlapping. This last property is very important in PPI networks since proteins may have multiple functions and therefore the corresponding nodes may belong to more than one cluster. ClusterONE also provides a p-value to select significant clusters, computed through a one-sided Mann-Whitney U test performed on the in-weights and out-weights of cluster vertices. A low p-value means that the in-weights are significantly larger than out-weights, so it is more likely that the cluster is a valid finding and not the result of random fluctuations. To distinguish significant clusters from insignificant ones, the authors suggest a p-value threshold of 0.05 [33].
Once the significant clusters are identified, in order to assess the meaningfulness of each cluster, KEGG and GO enrichment analyses are performed as described for the network as a whole (Section 4.2.1.). Next, we located cluster hub nodes (i.e., nodes having a number of neighbors higher than average) with consideration of the top 10% of the highest degree nodes. This threshold was suggested by previous works [34,35].

Bridge Nodes
Bridging centrality (BR) is a measure used in graph theory that can discriminate bridge nodes, i.e., the nodes that are crucial to dispatch information to the network topological structures. Bridge nodes usually have fewer neighbors than hubs, and are typically located between highly connected regions (i.e., network modules) [62][63][64].
The procedure of identifying bridge nodes begins with computing the BR measure for all network nodes. The BR of a node is defined as the product of the betweenness centrality (BC) and the bridging coefficient (BCoeff). BC measures the local properties of the network, while BCoeff measures global ones. The BR of a node n was calculated as proposed by Vitali et al. [34,35] following the equation: BR(n) = RWBC(n) × BCoeff(n) (2) where RWBC is the random walk betweenness centrality [65] and the BCoeff is computed as: where D(n) is the degree of node n (i.e., the number of neighbors of n), and N(n) is the set of neighbors of node n. Bridge nodes in the network are identified as the nodes whose BR values are in the highest percentiles (i.e., >90%) of the bridging centrality, similar to cluster hub nodes [34,35]. Because of their topological properties, the higher the BR of a node, the more information that flows through it. This makes bridging nodes good drug target candidates.
Once bridge nodes have been selected, a ranking score function is then applied to identify the most significant candidates for further in vitro experiments. For each bridge node b, we automatically retrieved its related KEGG biological pathways and GO terms. Next, by counting the number of pathways and terms related to wound healing, we obtained two first scores S KEGG (b) and S GO (b). These scores are then rescaled between 0 to 1 and a global score S tot is assigned to each bridge node b according to the formula: Based on the resulting ranked list, it is possible to plan in vitro experiments for the top scored genes.

How Do Genes React to Rigenera ® Stimulus? A Simulation of the Kegg Tnf Signaling Pathway
The fourth step of the proposed approach aims at investigating the effect of Rigenera ® on wound healing and suggesting other drug target candidates. To this end, we shift from a PPI representation of the process to a more detailed model. PPI networks can, in fact, provide a static view of the process, because the direction of the interaction (e.g., inhibition or activation) is often unknown. Instead, the molecular network underlying a specific process is highly dynamic: molecules, proteins; and genes associate, dissociate, and interact. In a previous work, we showed how signaling pathways from KEGG can be automatically modeled as BNs [34]. In BNs, nodes are genes and each regulatory reaction in a pathway is translated into a logic formula by parsing the KEGG pathway interactions, according to the conversion rules shown in Supplementary Table S21. BNs can be easily converted into Ordinary Differential Equations (ODEs) for simulation purposes.
The simulation of the Rigenera ® effect in the TNF signaling KEGG pathway (mmu04668) modeled as a BN is performed by initializing the input gene set and the bridge node states according to the treatment action. This means that their state is set to 0 or 1 if the drug, respectively, decreases or increases expression. Since only the initial states of these nodes in the pathway are known, following the procedure exposed in [34], we performed 1000 Monte Carlo simulations randomly assigning 0/1 states to the other nodes. Simulations were performed using the MATLAB-toolbox Odefy [66]. Finally, the states of all network nodes after the treatment stimulation are analyzed in order to identify other candidate key genes in the wound healing process. To rank gene candidates, we count the number of times that a gene changes its state from its initial value across all the Monte Carlo simulations, i.e., if the initial gene value passes a 0.5 threshold in any direction then the gene is considered switched.
Supplementary Materials: The following are available online at www.mdpi.com/1424-8247/10/2/55/s1. Table S1: KEGG pathways significantly associated to the whole network, Table S2: GO terms significantly associated to the whole network ; Table S3 to Table S16: GO/KEGG enrichment of the clusters; Table S17: KEGG pathways related to the wound healing and to at least one bridge node; Table S18: GO terms related to the wound healing and to at least one bridge node; Table S19: Bridge node score S; Table S20: Ordered list of genes changed over the simulations; Table S21: Table to convert KEGG associations into Boolean rules; and Method S1: Input gene panel setup.