Multiple sclerosis (MS) is a multifocal inflammatory autoimmune disease (1). Even though MS is usually considered a white matter disease, but several studies have demonstrated the involvement of gray matter impairment in conjunction with cortical and deep (2–6) leading to progressive neuronal damage in genetically sensitive hosts (1). MS is a complex multicomponent demyelinating disease and its pathophysiology consist of redox, autoimmune, vascular, and neurodegenerative systems, to name a few. The clear-cut mechanisms of MS triggering, its development, and progression are still obscure. In MS, impairing of the myelin sheath of neural axons in the Central Nervous System (CNS) is observed (7, 8). MS shows a long range of symptoms e.g., from pathological processes to severe physical disabling. Gender preferences, genetical factors, and geographical differences have been reported for people suffering from MS. Over outburst of MS, monocytes, which are a preserved subset of white blood cells, are activated by interferon-β (IFN-β) (7). MS study is quite vivid, using Omics data, many authors have used gene networks to get some insight into molecular mechanisms of MS (9–19). The integration of information gleaned from a variety of resources encompassing transcriptomics, genomics, proteomics and patient clinical data could boost our understanding of the mechanism(s) underpinning the reason for this disease (20). In this regard, we can explore the signaling pathways involved in MS (21), apply logical networks to model signaling pathways in MS (22) and use networks to combine information on transcriptome-interactome data from MS studies (17). We can also apply theory of biochemical systems for improving therapeutic drugs in re-myelination (15), create molecular networks based on transcription factors and genes expressed in mononuclear cells in MS patients (23), and design reactive networks between distinct miRNA and target genes in T cells (23). This approach will help explain the molecular mechanisms of the MS disease (12). Supplementary Table 1 shows some examples of network-based studies used with different MS biological data. As one can see, Bayesian Network (BN) modeling paradigms have rarely been applied in this setting. BN uses probability theory to reason under uncertainty. BN as a graphical scheme (directed acyclic graph) consists of a qualitative part (structural model) and a quantitative part (local probability distributions), which allow for a different kind of probabilistic inference, and quantitatively measures even the smallest impact of a variable or set of variables on others. This sort of modeling is of great importance in transcriptomic studies, since it can reveal both qualitative and quantitative elements of learned gene networks. BN has previously been used in several transcriptomic studies (10, 24, 25).

Many existing categories of gene networks identify groups of related genes as gene sets, making experimental follow-up a formidable task. With BN, it is possible to determine whether a gene is a driving source of changes in its gene network or not, since both in-degree and out-degrees of connectivity of each gene can be readily verified. The more out-degree gene has, the higher likelihood of being a possible regulator one. This would be a crucial characteristic for example when looking for potential drug targets. However, it is likely that a specific transcription factor defining a particular cell type that drives pathology, may not have a large number of out-degrees while still being crucial. To this end, if a particular gene is expressed across different cell types, for example like S100A6, then it may be correlated with different genes, but this may be a spurious correlation. Therefore, to leave off possible artifacts, we should use extra source of information when interpreting the results. Today, MS research is increasingly data-driven—a trend that arguably shall continue at a much higher rate in times to come. To tackle these large amounts of heterogeneous data, and to derive insight into MS disease, many interdisciplinary scientists have started using a variety of computational tools. In this study, we aim to gain much insight into the regulatory transcriptional gene network underlying MS using systems biology approaches in the context of BN, that may yield mechanistically interpretable results.

The fundamental idea behind this analysis was to shed some light into gene-gene interactions underpinning MS disease with regard to cause and effect (36). In this study, we reused GSE17048 experiment data which contained the profiled mRNA expression for all known genes in whole blood from 144 health individuals, 99 with MS (43 PPMS, 36 RRMS, and 20 SPMS). As meta-data of GSE17048 shows in the Gene Expression Omnibus–NCBI, in the conducting the experimental design, whole blood mRNA expression was compared between different types of MS and age-matched healthy control. The nature of probability distribution induced by a gene regulatory BN will allow diverse probabilistic gene queries to be answered in linear time. This makes BN to be practically appealing. The results of comparison of the network structures determined from various algorithms, including Hill Climbing, Tabu Search, Max-Min Hill Climbing, and Restricted Maximize algorithms with different scoring functions, are shown in Table 2. Some key properties of BN are fundamental in judging estimated results.

Figure 2 shows some properties of trained BN gene networks.

Topological parameters can characterize the location of genes in a gene network (37). Using NetworkAnalyzer, the following network topological parameters were calculated in our data. This was based on clustering coefficient (0.003), number of nodes (1,707), connected components (857), network diameter (26), network radius (1), shortest paths (192,463), characteristics path (8.313), the average number of neighbors (1.992), network density (0.0), isolated nodes (854), number of self-loops (0), multi-edge node pairs (0), and analysis times (1.467). The nature of probability distribution induced by a gene regulatory BN allowed diverse probabilistic gene queries to be answered in linear time. However, many structural BN parameters may be important. One of the key parameters (shown in Table 2) is the branching factor. This parameter plays a significant role in development of the gene network. Each Node (gene) will have its own branching measure, which will determine the out degree of that gene. If the branching factor value is not uniform in the network, an average branching factor can generally be calculated. This value turned out to be different depending on the type of algorithm used. Max-Min Hill Climbing returned a higher average than Restricted Maximize. In terms of system level understanding of research, the higher the branching factor, the more frequently gene regulators can be identified in the network. Biological networks have a modular architecture (38). MCODE can find connected and dense areas of the gene network based on network topology measures. In our analysis, 12 different modules were detected using MCODE, among which, 7 modules had 3 nodes; 3 edges with different interaction modes; 3 modules had 6 nodes and 7 edges; 1 module 15 nodes and 17 edges and finally 1 module had 6 nodes and 6 edges (Supplementary Figure 1). The active subnetworks were obtained using jActiveModules. The jActiveModules comprised 5 modules, where ILLMN_1742167 (TUBA1C), ILLMN_1665909 (LASP1), and ILLMN_1713636 (S100A6) were seen to be enriched modules (Figure 3). The number of modules detected by this method was different than those identified with the MCODE based method.

Figure 3 shows predicted modules in different modes of interaction. Module-level analysis explores the organization of biological systems and reconstructs module networks [A network module is a group of nodes (hub genes) that work together. Node or Vertex is a connection point or a branch point or an end point in a connection. And the path that connects the nodes to each other is called edge]. Figure 3 shows a module-level view of our gene regulatory BN network that denotes a high-level representation of the regulatory machinery of the MS gene network topology. Dense module searching of two MS Genome-Wide Association Study (GWAS) datasets identified several genes (GRB2, HDAC1, IL2RA, JAK2, KEAP1, MAPK1, RELA, and STAT3). These genes were enriched for glial cell differentiation (14). CytoHubba provides a user-friendly interface for discovering important nodes in biological networks (32). CytoHubba considers the shortest path between groups of nodes. Among the 11 proposed algorithms, MCC fitted better than the others. In Figure 4 and Table 3, we present the top 10 identified probes. Many of the genes, such as TUBA1C, LASP1, and S100A6 shown in Figure 4, are close to the hub genes and were actually identified as hub genes by other algorithms such as CyTargetLinker. The iRegulon software then allowed us to identify regulons using motif discovery in a set of regulated genes. Identified transcription factors affecting the hub genes are listed in Supplementary Table 2 and Supplementary Figure 2 and their common factors identified are given in Table 4. The most significant, the STAT5A protein, mediates the responses of many cell ligands, such as IL2, IL3 and different growth hormones. In this study, the gene identifiers were uploaded to Metascape and used in conjunction with KEGG pathways, GO biological processes, Reactome gene complexes, canonical and CORUM pathways (39). The results of the enrichment analysis, including descriptions, function, ontology, expression, etc. are shown in Table 5, Supplementary Table 3, and Supplementary Figure 3.

Genes were ranked from top to bottom based on degree, closeness and betweenness [higher degree (hub), higher betweenness (throat) and higher closeness centrality (shortest distance with other genes in the network)]. In terms of these parameters, three genes (LASP1, TUBA1C, and S100A6) showed a significant correlation with MS disease. These were thus identified as hub genes (Supplementary Table 4). In this study, probes ILLMN_1665909, ILLMN_1742167, and ILLMN_1713636 had high degrees of 15, 13, and 11, respectively and were identified as hub probes. In total of 850 probes had zero input edges, 200 probes had 1 in-degree. ILLMN_1665909, with the highest out-degree (13 out-degree) and 2 in-degree (mapped to human LASP1) plays an important role in regulating activity. Its encoded cytoplasmic protein binds focal adhesion proteins and plays a role in cell signaling, migration, and proliferation. ILLMN_1742167, with 12 out-degree and 1 in-degree mapped to the human tubulin gene (TUBA1C), and ILLMN_1713636 with 9 out-degree and 2 in-degree mapped to the S100A6 gene (Figure 2).

At present, the cause of MS is not fully understood, but knowledge of the genetic factors involved is essential for effective diagnosis and identification of the most appropriate MS therapeutic interventions. In this study, three genes (LASP1, TUBA1c, and S100A6) with high degree, high closeness centrality and high betweenness measures were highlighted as potential MS candidate regulator markers. These three genes (LASP1, TUBA1c, and S100A6) seem to be the most significant in the MS disease process. S100A6 functions in a wide range of cell types as a member of the S100 family and this family expression in MS patients could be considered as a diagnostic biomarker for MS. Its inhibition of demyelinating nerve cells suggests that S100 proteins could act as a candidate therapeutic target in MS (40). Komatsu et al. reported increased expression of S100A6 (Calcyclin), a calcium-bound protein of the S100 family, in human colorectal adenocarcinoma (41). Peterova et al. reported an overexpression of S100 protein-encoding mRNA in both colorectal cancer cell lines and surgically resected specimens of colorectal cancer (42). A study by Bartkowska et al. (43) showed that in response to different stress conditions, the level of S100A6 decreased in several brain structures, indicating that S100A6 may modulate stress responses. The genome-wide methylation array has identified a few hypomethylated immune-related genes, amongst them S100A6 which shows up-regulation in autoimmune encephalitis patients (44). Even though S100A6 is involved in many biological phenomena, its biological activity is still unknown (45). At the transcriptional level, upstream stimulatory factor and Nuclear factor-kappa B (NF-κB) activates the S100A6 gene promoter, although p53 might act indirectly to suppress transcription of the S100A6 gene (46). TUBA1C is a member of Microtubules which are vulnerable to degradation and disorganization in a variety of neurodegenerative diseases (47–49). Malfunction of microtubules (e.g., TUBA1C) is also considered as the central physiopathological mechanism of neurodegenerative diseases. The abnormalities in the regulatory pathways of microtubules disrupt the properties and functions of microtubules, leading to nerve damage (50). A decreased expression of the TUBA1C gene in Parkinson's disease has already been demonstrated by quantitative analysis of gene expression (51). LASP1)The LIM and SH3 protein 1), a focal adhesion adaptor protein, is an actin-binding, signaling pathway-regulated phosphoprotein which localizes within multiple sites of dynamic actin assembly. It has the potential to interact with various molecules, and is highly expressed in the adult CNS. Microarray data has revealed that alterations in LASP1 proteins affect cell migration, adhesion, and cytoskeletal organization (52). LASP1, significantly expressed by CNS neurons, is localized at synaptic sites (53).

A couple of significant transcription factors (TFs) that interact with these hub genes were identified in this study. The YY1 TF (Yi and Yang 1) is a multifactorial protein that, depending on the cell tissue, can activate or suppress gene expression (54). It is expressed in the nervous system. The YY1 promoter lacks the usual TATA box but has a rich GC sequence and therefore resembles a large subset of housekeeping and growth regulator genes. These features suggest that it may play an important role in development. In the CNS, myelination is performed by oligodendrocytes. YY1 function in oligodendrocytes was first reported by Berndt et al. (55). YY1 activates the promoter of myelin lipids and has been identified as an important player in myelination of the central nervous system during growth. In multiple neurodegenerative diseases, YY1 function is degraded through distinct mechanisms, including protein utilization, protein degradation, and ectopic nuclear/cytoplasmic shuttle (N/C). These disorders inhibit YY1 transcriptional activity and lead to gene transcriptional abnormalities that contribute to disease pathogenesis. A future goal in YY1 research is to discover other potential mechanisms that lead to YY1 dysfunction in neurodegenerative diseases, such as ectopic changes after translation (56). The other TF identified in study was Nuclear Factor of Activated T Cells 3 (NFATc3), a member of NFAT family. NFATc3 acts as signal integrators because their function is to bind STAT3, c-Jun, CREB, and ATF3 factors at specific DNA binding sites. NFATc3 cannot be regulated alone and act as calcium-dependent transcription factors. The antigen-mediated T cell receptor (TCR) mediates multiple signaling cascades, including phospholipase C (PLC) -dependent pathways that are secondary messengers of inositol-1,4,5-triphosphate (IP3) and diacylglycerol (DAG). IP3 binds to the IP3 receptor in the endoplasmic reticulum (ER) and releases Ca²⁺ ions into the cytoplasm (57). In this way, NFATc1-4, activates intracellular calcium via dephosphorylation (35). The findings show that NFATc3 is defined as a marker of a specific subset of astrocytes that are activated in response to lesions, as well as some degree of heterogeneity among astrocytes that may have consequences for cells in the nervous system (58). Preliminary findings in neuroblast cells have shown that various treatments that alter tubulin polymerization, such as reducing the mineral zinc, prevent the transfer of NFATc3 to the nucleus. In agreement with a functional relationship between NFAT and microtubules, it has been observed that the degradation of several proteins that control the proper organization of the microtubule network, and the actin-cytoskeletal linker, disrupts the nucleus and transcriptional activity of NFAT. Overall, it indicates the involvement of microtubules in NFAT nuclear stimulation (59). The LASP1 gene enhances NFAT2 nuclear translocation by activating the nuclear factor Akt (60). NFAT can affect processes such as axon growth, synaptogenesis, Schwann cell differentiation, and myelination (58). In general, it can be concluded that increase of the expression of LASP1 and S100A6 genes and decrease the expression of the TUBA1C gene in multiple sclerosis disrupts NFAT transcriptional activity. Although the role of NFAT in regulating the immune system is well established, our knowledge of NFAT in human disease is limited. The function of NFAT in other aspects of human immune or inflammatory diseases is also largely unknown (61).

The involvement of hub genes identified in this study in other disorders have been reported as well. Patients with MS are known to suffer from a number of digestive problems (62) and studies have shown that LASP1 (63) and S100A6 genes have high expression in the digestive system. A link can therefore be established between the expression of these genes, MS, gastrointestinal problems and possibly other types of human cancers (64). Also, LASP1 plays a crucial role in the growth and metastasis of gastric cancer and other cancers (52, 63, 65–68). For example, LASP1 can cause the progression and metastasis of colorectal cancer (CRC), but its mechanism is still unclear (69). A connection between LASP1 and S100A has reported underpinning LASP1 binds to the calcium-binding protein family (S100A) and increases its expression in colon cancer (Kappa = 0.347, P < 0.01) (70). On this basis, the present study confirmed the importance of three gene expression patterns (LASP1, S100A6, and TUBA1C) for understanding the transcriptome complexity of MS. This leads us to conclude that upregulation of LASP1 and S100A6 genes along with down-regulation of TUBA1C is central to MS pathology. To our knowledge, this is the first report to evaluate the level of expression of the above genes for discovery of a transcriptomic signature for MS disease. These findings provide a potential mechanism for some significant biomarkers responsible for the pathogenesis of MS. However, we still have a long way to go to understand the larger transcriptomic profile for this disease. This study provides initial data to further investigate the possible role of these genes in the pathogenesis of MS.

Results of the present study indicate that the analysis of gene expression data based on gene-gene interaction networks can provide opportunities to determine the genes involved in MS. The importance of three candidate marker genes in this disease were highlighted. These candidate marker genes, LASP1, TUBA1C, and S100A6, identified by the biological systems approach, have been further confirmed in the laboratory. The significant difference in the expression of these three genes in patients with MS will help further research on this disease and its treatment. This useful tool can serve as a good starting point for identifying new therapies and understanding the basic mechanisms controlling normal cellular processes and disease pathologies. It is crucial to point out here that for learning Bayesian gene network in this study, we did not separate sets of possible signaling protein molecules and transactional factors beforehand in our data, and consider them to be parents (causatives) in the learned network. By doing so, the learned Bayesian gene network probably would be biologically much more appealing. We aim to do this in a due course in the future.

The data presented in the study are deposited in the https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE17048 repository, accession number GSE17048.

All procedures were approved and carried out in accordance with Medical Research Ethics Committee in Iran under code IR.UI.REC.1399.076. The patients/participants provided their written informed consent to participate in this study.

NK and MG-Z developed the theoretical formalism, designed the model and the computational framework, and analyzed the data. NK performed experimental lab the analytic calculations. MM-B conducted the backbone of the experiment and contributed to the final version of the manuscript. All authors have read and agreed to the published version of the manuscript.