Blow-up of energy solutions for the semilinear generalized Tricomi equation with nonlinear memory term

1.   Introduction

Obesity and T2D are two chronic metabolic and endocrine abnormalities that are on the rise all over the world. High body mass index ratio is known to contribute to β-cell function decline and inadequate insulin secretion, which ultimately leads to T2D ^[1,2]. The development of insulin resistance in obese individuals is due to the elevated levels of non-esterified fatty acids (NEFAs), cytokines, hormones, and inflammatory substances. The NEFAs secreted from adipose tissue in obese individuals are the key molecular etiologic factors connecting impaired β-cell function and insulin resistance to the risk of T2D. Obesity-related T2D may further increase the risk of developing comorbidities like cardiovascular diseases, hypertension, neurological and neuropsychiatric illness, and autoimmune diseases ^[3–6]. No effective therapy for T2D caused by obesity is available, except the adoption of disease management strategies aimed at controlling blood glucose levels and secondary complications.

Molecular studies have reported the aberrant expression patterns of different genes in adipose tissues and blood samples of obese subjects ^[7] and also in blood samples of diabetics ^[8]. Although, several unique markers for these metabolic conditions have been identified, but they share only few disease genes like PPARG ^[9] and UCP3 ^[10], ENPP1 ^[11] and FTO ^[12] between them. The critical genes and mechanisms linked to obesity and T2D are still far from being fully understood. Therefore, identifying key molecular mechanisms between obesity and T2D may not only advance our understanding about the pathophysiology of these common health conditions, but also contribute to developing new treatment strategies.

In recent times, next generation sequencing methods like RNAseq have provided an advantage to perform the unbiased molecular interrogation of several genes in studying complex diseases ^[13]. NCBI hosted Gene Expression Omnibus (GEO) database consists of high-quality transcriptomics data generated both from microarrays and RNAseq. The enormous amount of publicly resource gene expression data provides a valuable opportunity to identify potential key genes contributing to obesity and T2DM. Furthermore, in recent years, gene network-based bioinformatics assessment of differentially expressed gene signatures is gaining attention as a valuable tool in studying molecular drug targets and to uncover many druggable genes ^[14]. Exploring the novel therapeutic targets or molecular pathways for the existing drugs could also help in patient subgrouping and developing personalized treatment ^[15].

Owing to the sparse data available, the objective of present study was to evaluate the shared genes, gene networks, molecular pathways, and drug targets between obese and T2D subjects. In this pursuit, this study has analyzed the RNAseq transcriptomic data of adipocytes of obese, lean, and diabetic subjects with the help of differential gene expression, integrated protein interaction networking, and molecular drug target identification strategies. The potential genetic markers identified among obese individuals at risk of developing T2D may have the potential to develop targeted therapies in the future.

2.   Materials and methods

The overview of methods used in the present study is represented in the Figure 1.

2.1.   Clinical Samples and RNASeq Expression

The gene expression datasets, comprised of adipocyte mRNA expression profiles from lean (6 samples), obese (6 samples) and T2D (6 samples) subjects, were obtained from the gene expression omnibus (GEO) under the accession number GSE133099 ^[16]. The information about samples is given in the Supplementary File (Supplementary File S1). The details of the original experiment, including RNA isolation, library preparation, sequencing, and data analysis, are yet to be published. The RNASeq reactions of the clinical samples were performed using an Illumina HiSeq 2500 machine.

2.2.   RNASeq Data analysis

The GEO expression data had been quality checked using FastQC ^[17] followed by Trimmomatic ^[18] methods for removing adaptor contents and low quality reads with a Phred Score of < 20. High quality reads were aligned to the hg38 reference genome using STAR ^[19], and mapped reads were counted by Feature Counts ^[20]. The genes that had zero counts were removed from the analysis. Later on, the differential expression analysis was performed using DESeq2 ^[21] based on the Wald test scores followed by the Benjamini-Hochberg procedure for removing false positives in the data. The genes that had an absolute fold change (FC) of ≥ 1.5 and an adjusted p value of ≤ 0.05 were screened as differentially expressed genes (DEGs). The heatmap and volcano plots representing DEGs were created using r packages like pheatmap ^[22] and Enhanced Volcano ^[23], respectively.

2.3.   Construction of the Protein Interaction Network

The protein interaction data was retrieved from the Human Integrated Protein-Protein Interaction Reference (HIPPIE) database ^[24,25]. HIPPIE collects the experimentally confirmed interactions from various sources, which include HPRD, BioGrid, IntAct, DIP, MINT, MIPS, and BIND. All the protein interactions of DEGs obtained from the RNAseq analysis were extracted with a default association score of ≥ 0.4 to create the protein interaction network (PIN). Visualization and calculation of topological parameters of PIN were performed using Cytoscape ^[26] (version 3.8.2). To generate the topological parameters, the PIN was analyzed with the help of Cytoscape plugins like NetworkAnalyzer and CytoNCA ^[27].

2.4.   Identification of key genes in the network

The local topological parameter degree ^[28] and global topological parameters like betweenness ^[29], closeness ^[30] and eigenvector ^[31] centralities were used to identify the most influential nodes in the network.

2.4.1.   Selection of hubs

The genes with high degree of connectivity were considered as hubs. The following formula ^[32], was used for identifying hub genes, where SD denotes their standard deviation.

2.4.2.   Selection of bottlenecks

Using the node betweenness, closeness and eigenvector distribution, genes positioned in the top 90^th quantile were scaled as bottlenecks. The formula for calculating each parameter is as follows:

2.4.3.   Betweenness centrality

where '$s$' and '$t$' are nodes in the network other than '$n$'. ${\sigma }_{st}$ represents the number of shortest paths from '$s$' to '$t$', and ${\sigma }_{st}\left(n\right)$ is the number of shortest paths from $s$ to $t$ that '$n$' lies on ^[33].

2.4.4.   Closeness centrality

The measure of how quickly information flows from a given node to other reachable nodes in the network is given by closeness centrality (CC). The CC of a node n is the reciprocal of the average shortest path length and is computed as follows:

where L(n, m) is the length of the shortest path between two nodes n and m.

2.4.5.   Eigenvector centrality

αmax is the eigenvector corresponding to the largest eigenvalue of the adjacency matrix A.

2.5.   Downstream analysis of key gene signatures

2.5.1.   Functional enrichment analysis

To find out gene ontologies, pathways and diseases associated with the key genes, we filtered them based on topological parameters using the ToppGene ^[34] suite with a threshold of q value of ≤ 0.05. ToppGene uses gene expression, protein domains, protein interactions, ontologies, pathways, phenotypes, and text mining to identify functional enrichment of the input gene list.

2.5.2.   Tissue Specific Analysis

Tissue specific functional gene network modules were constructed using HumanBase ^[35]. This tool presents genome-wide functional interactions of gene networks using Bayesian methodology that is driven by the integrated data of diverse experiments spanning tissue and disease states. We have queried the key gene list derived from topological parameters (hubs, bottlenecks, closeness, and Eigenvalue) in HumanBase to find out their changing functional roles across adipose and pancreatic tissues, at a default q-value threshold of ≤ 0.05.

2.5.3.   Druggability analysis of hub and bottleneck genes

The potential of the drug targets among the hub and bottleneck genes has been analyzed in Open Target Platform ^[36]. A ≥ 0.1 cutoff association score was used to detect the expression status of druggable molecular targets in the adipose tissues of individuals with obese and T2D conditions.

3.   Results

3.1.   Analysis of differentially expressed genes

The analysis of raw count data yielded the expression profiles of 34328 molecular biotypes, including protein coding mRNAs, miRNA, lincRNA, etc., from, which only protein coding (17222) genes were retained for further analysis. Differential expression of genes was performed using DEseq2 to identify the DEGs from 'obesity vs lean' and 'T2D Vs lean' sample comparisons. DESeq2 employs a negative binomial distribution to pick statistically significant differentially enriched genes from an over dispersed data, as is often the case with biological count data. The volcano plot distinguishing significant DEGs in each comparison are represented in Figure 2 A and B. The T2D and lean group analysis identified 1381 DEGs (FC ≥ |1.5|, adj p-value of ≤ 0.05), comprising of 680 upregulated and 701 downregulated genes (Figure 2C). On the other hand, for the obese vs lean comparison, 1281 DEGs (FC > |1.5|, adj p-value of ≤ 0.05) including 759 up and 522 down regulated genes were identified (Figure 2C). The list of the top 10 differentially regulated genes in each category is given in Tables 1 and 2. The comparison of DEGs across T2D and obesity conditions, has revealed the shared expression pattern of 730 (430 up- and 300 down- regulated) genes (Figure 2D). Interestingly, these 730 shared genes presented similar expression patterns, i.e. either up or down regulated in both conditions, underlining the potential molecular connections between T2D and obesity.

3.2.   The protein interaction network

The experimentally validated protein interaction network downloaded from the HIPPIE database was used to make the protein interaction network of T2D and obesity. All the 1932 non-redundant DEGs obtained from both T2D, and obesity conditions were queried in the HIPPIE database to generate their molecular interactome. From the protein network outputs generated by the HIPPIE database (Figure 3A), we have chosen ≥ 0.4 as the cut-off score to generate protein interaction networks (PIN) (Figure 3B). The constructed PIN entailed 1203 nodes and 3211 edges, where nodes indicate proteins and edges indicate the interaction among the proteins. The mean degree of the network was 5.34, and the standard deviation was 7.88.

3.3.   Accessing Network Centrality parameters of the protein interaction network

The network visualization and the calculation of topological parameters were performed in Cytoscape. Local parameters like degree of network and global parameters like betweenness, closeness, and eigenvector centrality were considered to filter the DEGs from the PIN. The nodes with a threshold degree of ≥ 21 were considered hubs based on the mean and standard deviation values of the degree of the network. For global parameters, nodes that were in the 90^th percentile were screened. Finally, 36 hubs, 121 high betweenness, 121 high closeness, and 120 high eigenvector nodes from the network were obtained. There were common nodes among the nodes generated from global parameters due to the higher correlation between the global parameters (Table 3). The plot representing the higher correlation among global parameters is represented in figures 3C and 3D. A total of 190 genes remained after filtering them based on topological parameters and removing the duplicates. Among 190 genes, 70 genes (36.84%) were shared between T2D and obesity. The common genes in each topological parameter are represented in figures 4A and 4B. The top 10 hubs in shared genes based on degree are represented in Table 3. As shown in the heatmap (Figure 5), the pattern of expression of filtered genes was similar in both T2D and obesity but had a clear distinction from the lean samples.

Response to lipid, regulation of cell death, regulation of cell differentiation, regulation of insulin, regulation of phosphorylation, fatty acid transport, response to hormone, glutamate receptor binding, response to glucose (p < 0.05) were the enriched terms for T2D and obesity genes. Functional enrichment indicates major deregulation in insulin and glucose-related metabolism. Obesity has a significant impact on insulin-sensitive tissues, such as adipose tissue and the pancreas, at both biomolecular and functional levels.

The tissue based functional modules in the pancreas and adipose tissue have identified 6 and 8 modules, respectively (Figure 7). The modules detected in adipose tissue with a significant p-value of < 0.05 were enriched in different biological processes such as glyceraldehyde-3-phosphate metabolic processes, cholesterol biosynthetic processes, steroid metabolic processes, and lipid biosynthetic processes. Similarly, the modules of pancreatic tissue were enriched in triglyceride metabolic processes, amyloid precursor protein metabolic processes, and regulation of lipid metabolic and biosynthetic processes. Tissue level enrichment has shown major deregulation of glucose, lipid, and cholesterol related metabolic pathways, which equally contribute to obesity and T2D ^[37].

3.4.   Ranking the key genes based on topological analysis

The genes that were filtered based on protein interaction network and topological parameters like degree, betweenness, closeness, and eigen vector have shown their relatedness to metabolic disorders like T2D and obesity. We next categorized the genes into four categories: Q1, Q2, Q3, and Q4. If a gene is filtered by all four topological parameters, it is categorized as Q1. Similarly, Q2, Q3, and Q4 were assigned according to the number of topological parameters used to filter it. Of all the 4 categories, we focused on Q1 (TGM2, HIST1H1B, HIST1H3H, MEOX2, HIST1H2BE, ARRB1, HSPA1A, FUS, CHD3, FYN, ERBB2 and FN1) and Q2 (HBA1, EEF1A2, FASN, ACACA, HSPA6, IRS1, SVIL, PFN2, ITGB1, SIPA1L1, NCOA3, PLAUR, HP) gene sets for further analysis as they had higher confidence since they fall under multiple topological parameters.

3.5.   Mapping the traits of the signature genes

The Q1 (12 genes) and Q2 (13 genes) class genes were queried in Open Targets platform to map the reported traits. In both T2D and obesity conditions, we discovered that 18/25 (72%) genes were functionally correlated. The snapshot of the association of these 18 genes and GWAS traits related to T2D and diabetes is shown in Figure 8. Genes like ACACA, FASN, FUS, HBA1, HSPA1A, NCOA3, TGM2, EEF1A2, HIST1H1B and HIST1, H2BE were known for their association with body mass index (p ≤ 0.05). The additional associations observed are as follows; IRS1 and ERBB2 in triglyceride levels; SIPA1L1 in bone mineral density (p≤0.05); ARRB1, ERBB2, IRS1, and NCOA3 in HDL cholesterol levels (p≤0.05); CHD3, FN1 and HIST1, H2BE in LDL cholesterol (p≤0.05). Moreover, the association of ACACA, EEF1A2, HBA1, HSPA1A, IRS1, MEOX2 and PFN2 genes with Type 2 Diabetes (p≤ 0.05), FUS with obesity, and CHD3, EEF1A2, HIST1, H2BE, HSPA1A and PFN2 with hypertension (p≤ 0.05) was also reported. The genes screened based on the topological parameters have shown close association of the traits related to T2D and obesity.

3.6.   Druggability analysis

Druggability analysis in the Open Target database of the hubs and bottleneck genes (Q1 and Q2 categories) selected from the topological analysis found 15/25 (60%) genes with a genotype-phenotype association score of > 0.1. Six of these 15 genes (ERBB2, FN1, FYN, HSPA1A, HBA1, ITGB1) were targeted by small chemicals, antibodies, and/or enzyme molecules (Table 4). The ERBB2 gene had the highest number of known drugs, followed by FYN (6) and FN1 (5) (Supplementary File S2).

The ERBB2 gene has 36 known drugs, of which 11 receptor protein-tyrosine kinase inhibitor molecules (Afatinib, Afatinib Dimaleate, Dacomitinib, Lapatinib, Tucatinib, Lapatinib Ditosylate, Vandetanib, Neratinib Maleate, Dacomitinib, Neratinib) are currently in phase 4 clinical trials. For the FYN gene, Dasatinib, a small molecule inhibitor, is currently under phase IV trials for different phenotypes. For the FN1 gene, Ocriplasmin, a fibronectin proteolytic enzyme, is currently undergoing phase IV clinical trials. For the ITGB1 gene, there are ongoing clinical trials with Volociximab antibody (an Integrin alpha-5/beta-1 antagonist) in phase 1 and with Firategrast small molecule inhibitor (Integrin alpha-4/beta-1 antagonist) in phase 2.

The Open Target analysis identified six hub genes that are druggable. These genes, in total, have 52 known drugs in the market. Currently, some more compounds targeting them are undergoing clinical trials for an assortment of disease indications. After studying their efficacy in vitro and in vivo experiments, 52 known drugs can be repurposed to manage T2D and obesity. Drug repurposing is a rapid process to develop revolutionary treatments because the new indication is based on already accessible safety, pharmacokinetic, and manufacturing data.

4.   Discussion

Obesity and type 2 diabetes mellitus are the most common progressive metabolic disorders with complex molecular pathology involving defective lipid metabolism, insulin resistance, and beta-cell functioning. Although previous investigations have reported various potential molecular markers linked with the progression of obesity and type 2 diabetes mellitus, the key common genes, the molecular mechanisms underlying its pathogenicity remain elusive ^[38]. In this study, by using RNA-sequencing data, we identified 1932 unique genes expressed across two comparisons (a) diabetes vs lean healthy control (680 up-and 701 down-regulated genes) and (b) obese vs lean healthy control (759 up and 522 down regulated genes) groups. In recent years, some studies have performed gene expression analysis of obesity ^[7,8] and diabetes ^[39]. But variability in study objectives, datasets, experimental platforms, statistical measures, and a broad range of bioinformatics methods has added more complexity to data analysis and interpretation ^[7,40–42].

Network biology provides an excellent platform for investigating the dynamic interactions that exist between different genes ^[43]. Network biology implements graph theory to reveal molecular interactions and enrichment patterns that are otherwise not solvable by univariate analysis. In this study, the PIN of 1932 genes entailed 1203 nodes and 3211 edges, where nodes indicate proteins and edges indicate the interaction among the proteins. It is important to study the topological changes in biological networks as they are context-specific and dynamic. A divergence from normal regulatory network topology could indicate the pathogenic mechanism, and the genes with the highest network topological changes could be used as biomarkers for disease or as targets for drug development or therapeutic intervention ^[44]. A differential network analysis approach was effectively used to identify 15 potential key players or disease markers in hepatocellular carcinomas ^[45]. By applying different local and global centrality parameters (degree, betweenness, closeness, and eigen vector), we identified 10 hub genes (FN1, FYN, ERBB2, HIST1H3H, CHD3, HSPA1A, ARRB1, HIST1H2BE, HIST1H1B, and FUS) shared between obesity and T2D. Functional enrichment of hub genes represented their involvement in insulin, fatty acid and glucose metabolism, supporting the findings from other studies ^[46]. A recent computational study has identified the differential expression of CEBPD, TP73, ESR2, TAB1, MAP3K5, FN1, UBD, RUNX1, PIK3R2 and TNF genes in tissue samples of obese samples by employing different methods like DEGs, GO and reactome pathway enrichment analysis and transcription factor mapping ^[47].

Human diseases are caused mainly by the disordered balance of tissue-specific processes, and the actual function of genes is heavily dependent on their tissue context ^[35]. Recently, tissue-specific functionally interacting network modules have been developed to discover phenotype-relevant network modules enriched in gene expression networks ^[48]. In this study, we identified the enrichment of network modules connected to different biological processes such as the glyceraldehyde-3-phosphate metabolic processes, cholesterol biosynthetic processes, steroid metabolic processes, and lipid biosynthetic processes in adipose tissue. Whereas the modules reported in pancreatic tissue were enriched in triglyceride metabolic processes, amyloid precursor protein metabolic processes, and regulation of lipid metabolic and biosynthetic processes. Our tissue based functional enrichment of the biological network underscores the disturbances of lipid and energy metabolism in both diabetes and obesity. These findings are supported by other studies ^[49].

Computational network biology measures like hubs are proven to aid in discovering novel drug targets as they can link complex molecular interactions ^[50]. The Open Target Platform identified 6 genes (ERBB2, FN1, FYN, HSPA1A, HBA1, ITGB1) as tractable targets by small chemical, antibody, and/or enzyme molecules. The ERBB2 gene encodes the HER2 receptor, which is related to diabetes and insulin resistance through fatty acid synthase (FASN) activity ^[51]. ERRB2 has also been suggested to play a role in preadipocyte differentiation and obesity ^[52]. The second tractable target, fibronectin, encoded by the FN1 gene, is an extracellular matrix protein that binds to integrins and plays a major role in cell adhesion, growth, migration, and differentiation. Previous research has revealed that fibronectin in adipose tissue regulates adipocyte-specific gene expression and may play a pathophysiological role in obesity-related comorbidities in humans ^[53]. Furthermore, endothelial dysfunction causes insulin resistance, which leads to FN1 gene overexpression and, eventually, an increase in fibronectin levels in the blood ^[54]. Fyn is a tyrosine kinase that promotes the interaction of PIKE-A and STAT5a, thereby regulating insulin signaling and adipogenesis ^[55]. Fyn knockout (Fyn/--) mice have a lean phenotype, which links this kinase to obesity development ^[56]. Fyn inhibition has been shown to prevent the maturation of 3T3-L1 preadipocytes into fully mature adipocytes ^[57] and to increase in vivo fat loss by improving insulin sensitivity and modulating AMPK activity via LKB1 function ^[58]. HSPA1A is a stress inducible protein with functions ranging from cell signaling, immune response, and chronic conditions. Extracellular HSPA1A levels (eHSP70) were found to be higher in obese patients with T2DM when compared to non-obese T2DM patients ^[59]. In type 2 diabetes, pharmacological (e.g., the hydroxylamine derivative BGP-15, which is currently undergoing clinical trials) and physiological (hyperthermic) treatments have begun to be viewed as promising therapeutic options to target the HSPA1A gene ^[60]. Therefore, the 6 drug targets identified in this study, are in line with previous reports indicating their potential drug tractability for obesity and/or T2D.

Our approach, however, has certain limitations. The number of samples that we analyzed here is low, and the samples were from the female population (lean, obese, and T2D with obesity). But, given that we used the secondary data downloaded from GEO, and our lack of control over the study design, this is an unrealistic limitation. We agree that analysis of more samples will help in reducing the background noise and reducing the standard error of the effect fraction, but it may not significantly change the conclusion of this study. Many sex hormone responsive pathways and gene networks are known to undergo gene expression changes between sexes, especially in adipocyte function, inflammation, which leads to differential adipose tissue remodeling ^[61]. We recognize the female only study may be biased significantly towards these gene network differences. To overcome these limitations to some extent, we have used very stringent thresholds for the downstream analysis in the study to remove noise in the data. Furthermore, all the statistical analysis of the data hosted in public databases was carried out on a global scale, irrespective of any ethnic population.

5.   Conclusions

This paper lays forth a comprehensive bioinformatics method for identifying the most important key signatures shared between obesity and T2D from the experimentally validated Protein Interaction Network. A rigorous parametric downstream analysis based on biological insights reveals 10 candidate hub genes (FN1, FYN, ERBB2, HIST1H3H, CHD3, HSPA1A, ARRB1, HIST1H2BE, HIST1H1B, FUS) shared between obesity and T2D that could be regarded as new genetic biomarkers for T2D and obesity. All these hub genes are involved in insulin regulation, glucose response, triglyceride transport, and fatty acid transport. Additionally, we identified six genes (ERBB2, FN1, FYN, HSPA1A, HBA1 and ITGB1) as tractable targets by small chemicals, antibodies, and/or enzyme molecules. The exact functional role of these potential genes can be biologically validated by implementing suitable in vitro and in vivo experimental approaches. Overall, by implementing a broad range of computational biology methods, our study highlights the potential of correlating gene expression to topological parameters, functional relationships, druggability, and the identification of candidate genes involved in the metabolic disorders with common etiological basis.

Author contributions

AB, NS and BB: conceptualization; AB and BB: data curation; AE and BB: formal Analysis; AB: funding acquisition; NS and BB; methodology; AB: project administration; SP, BB: Resources; SP, BB: software; AB, NS and BB: supervision; WA, BB: validation; BB: visualization: AB, AE, ZK, TS, PK, NS, BB; writing original draft: AB, PK, SA, RE, NS, BB; writing, review, and editing: AB, AE, WA, ZK, SP, TS, PK, SA, BB and NS. All authors read and approved the final version of the manuscript.

Acknowledgments

This Project was funded by Deanship of Scientific Research (DSR), King Abdulaziz University, Jeddah, Saudi Arabia under grant no. (KEP-1-140-41). The authors, therefore, acknowledge with thanks DSR technical and financial support.

Conflict of interest

Authors declare no conflict of interests for this article.

Data availability statement

All datasets analyzed for this study are included in the article and Supplementary Material.