Systematic analysis of transcription-level effects of neurodegenerative diseases on human brain metabolism by a newly reconstructed brain-specific metabolic network.

Network-oriented analysis is essential to identify those parts of a cell affected by a given perturbation. The effect of neurodegenerative perturbations in the form of diseases of brain metabolism was investigated by using a newly reconstructed brain-specific metabolic network. The developed stoichiometric model correctly represents healthy brain metabolism, and includes 630 metabolic reactions in and between astrocytes and neurons, which are controlled by 570 genes. The integration of transcriptome data of six neurodegenerative diseases (Alzheimer's disease, Parkinson's disease, Huntington's disease, amyotrophic lateral sclerosis, multiple sclerosis, schizophrenia) with the model was performed to identify reporter features specific and common for these diseases, which revealed metabolites and pathways around which the most significant changes occur. The identified metabolites are potential biomarkers for the pathology of the related diseases. Our model indicated perturbations in oxidative stress, energy metabolism including TCA cycle and lipid metabolism as well as several amino acid related pathways, in agreement with the role of these pathways in the studied diseases. The computational prediction of transcription factors that commonly regulate the reporter metabolites was achieved through binding-site analysis. Literature support for the identified transcription factors such as USF1, SP1 and those from FOX families are known from the literature to have regulatory roles in the identified reporter metabolic pathways as well as in the neurodegenerative diseases. In essence, the reconstructed brain model enables the elucidation of effects of a perturbation on brain metabolism and the illumination of possible machineries in which a specific metabolite or pathway acts as a regulatory spot for cellular reorganization.

Introduction

Neurodegenerative diseases are disorders associated with functional loss of brain cells, especially with advancing age. Reported increase in the average life expectancy throughout the world means a large increase in the number of people with these diseases in the coming decades [1-3], making it essential to understand how these diseases develop. A full understanding of the mechanisms of neurodegenerative diseases will lead to improvement of successive treatment techniques.

Metabolism has an essential role in the cell, therefore, it is of substantial importance to perform a systematic investigation of the effect of neurodegenerative diseases on brain metabolism. Such an analysis can also reveal the common pathways perturbed by neurodegenerative diseases as well as disease-specific metabolic pathways. Limited research has been performed on this, however, a systematic and comparative analysis is still missing [4,5]. Such an analysis necessitates a detailed and curated representation of metabolic reactions occurring in the brain. Consequently, brain-specific metabolic stoichiometric models have recently been reconstructed [6-8]. These models can be used as scaffolds for computational analyses to draw hypotheses, as done for other metabolic systems in literature [9,10].

In this study, we have developed a novel brain-specific stoichiometric metabolic model which covers a large number of metabolic pathways in and between astrocytes and neurons, the two major cell types in the brain. We then used this model as a scaffold to analyze transcriptional changes associated with six common neurodegenerative diseases: Alzheimer’s disease (AD), Parkinson’s disease (PD), Huntington’s disease (HD), multiple sclerosis (MS), amyotrophic lateral sclerosis (ALS), and schizophrenia (SCH). The changes induced by these diseases on brain metabolism at the transcriptional level were systematically and comparatively analyzed by integrating corresponding transcriptome data of the diseases with the here-developed metabolic network of brain. This enabled a comprehensive analysis of the effect of neurodegeneration on metabolism from the perspective of different related diseases.

Results

Validation of the reconstructed model

We have reconstructed a novel brain metabolic model based on a previous study [6] by (i) transforming the lumped reactions in the original model into elementary reactions and (ii) adding new literature-supported pathways, as described in detail in Methods section. The new brain model comprises 630 reactions (571 internal, 59 exchange) and 524 metabolites (465 internal, 59 external). This means a 3-fold increase in the number of reactions compared to the original model [6]. 253 of the reactions occur in neurons while 299 of them are astrocytic. 59 reactions are for the exchange of metabolites, and there are 19 reactions depicting the transfer of metabolites between the two cell types. The model is available both in excel format and in SBML format as a supplement (Supplementary Files 1 and 2). 460 out of 571 internal reactions have gene information obtained from HumanCyc [11]. Due to the high coverage (81%) of gene-associated reactions in the model, it is possible to use this brain model for the integrated computational analysis of brain metabolism and transcriptome data. Since the total number of genes in the model is 570, it was named iMS570. Some of the reactions are controlled by one or more genes whereas some of the reactions do not have gene information. The number of the matched genes is 496 (87%) covering 447 of 571 internal reactions for GEO: GSE26927 dataset used in this study.

Note that the present model, iMS570, was developed purely in a brain-specific manner. Its gene coverage is much higher than the other gene-associated brain metabolic model (iNL403) reported in literature so far (570 vs. 403) [7], which was derived from a generic human metabolic reconstruction. The curation of the present model was achieved in extensive iterative steps such that it can be used in flux-balance based simulations. Flux balance analysis [12] of resting healthy brain metabolism performed with the newly reconstructed model iMS570 revealed a meaningful flux distribution consistent with experimental major flux splits in central metabolism (Table 1). Similar to the original model [6], the resting-state fluxes were calculated using two consecutive simulations: (i) linear programming with the objective function of maximization of the sum of these fluxes through glutamate/glutamine/GABA cycles (primary objective), (ii) after fixing the value of the sum of these three cycle fluxes to the optimum value, quadratic programming with the objective function of minimization of the Euclidean norm of fluxes (secondary objective). The primary objective provides a tight coupling between the two cell types, and it was shown to outperform several other alternative objective functions [6]. On the other hand, linear objective functions usually result in non-unique flux distributions, known as the problem of alternate optima [13]. Our secondary objective function eliminates the alternate optima, resulting in a unique flux distribution [6]. Both of the objective functions have physical meanings. The first one ensures a tight coupling of astrocytes and neurons through the well-known intercellular cycles [14,15]. The second one ensures channeling of fluxes through all pathways with minimal use of enzymes [16,17]. The FBA-based prediction of metabolic fluxes at resting state when combined with a proper objective function is one of the superiority of the iMS570 model over the other comprehensive brain model available in literature [7]. Besides, the FBA results render a verification of the reconstructed model.

Table 1

Predicted flux results in this study, in comparison to the experimental results in resting (healthy) state. Predictions by the original model [6] are also given for comparison.

% Flux ratio	This study	Original model	Experimental
Lactate release flux (r11) with respect to CMRglc	7.2	4.5	3–9
			[132–135]
Glutamate/glutamine cycle flux (r95) with respect to CMRglc	73.1	68.0	40–80
			[14,136,137]
Relative oxidative metabolism of astrocytes (rTCA,A/rTCA,total, r25/(r25 + r69))	33.9	35.0	30
			[136,138,139]
Total lipid synthesis with respect to CMRglc	3.1	2.8	2
			[140]
Total PPP flux with respect to CMRglc	5.5	5.6	3–6
			[141,142]
Pyruvate carboxylase flux (r12) with respect to CMRglc	11.1	11.7	10
			[136,138,143]

The model was reconstructed in such a way that pathways which are of major importance to brain metabolism are included. Although there are still reactions/pathways not covered by iMS570, we have performed an additional analysis to show that missing reactions will not affect current results and conclusions. To this aim, we removed 72 reactions from iMS570 belonging to the purine and pyrimidine nucleoside metabolisms, and repeated the reporter metabolite and reporter metabolism calculations for PD. Among about 40 associated reporter metabolites (Table 2), only three (astrocytic ADP, ATP and GTP) were not captured in this highly reduced model. Additionally, only methionine metabolism (Table 3) was not any more a reporter metabolism in addition to the excluded pathways. This clearly shows that pathways which may not be included in the current reconstruction do not have any practical effect on the conclusions made.

Table 2

Reporter metabolite analysis results in different neurodegenerative diseases. “+” sign states significant changes with p-values less than 0.05 (reporter metabolites), and “−” sign represents no significant change (A: astrocytes, N: neurons).

Metabolite	AD	ALS	HD	MS	PD	SCH
ADP/ATP (A)	+	−	+	+	+	−
ADP/ATP (N)	+	+	+	+	+	+
AMP (A)	−	−	+	−	−	+
AMP (N)	−	−	−	−	−	+
GTP (A,N)	−	−	−	−	+	+
NADmit/NADHmit (A)	−	−	+	−	−	−
NADmit/NADHmit (N)	−	−	−	−	+	−
NADP/NADPH (A)	−	−	−	+	−	−
NADPmit/NADPHmit (N)	−	−	−	−	−	+
CytCox (A,N)	+	+	+	+	+	−
CytCred (A,N)	+	+	+	+	+	−
Ubiquinol (A,N)	+	+	−	−	+	−
Ubiquinone (A,N)	+	+	−	−	+	−
Hc (A,N)	+	+	+	+	+	+
Pyruvate (A)	−	−	−	−	+	−
Pyruvate (N)	−	−	+	−	+	+
Fumarate (N)	−	−	−	−	−	+
Malate (A,N)	−	−	−	−	+	−
Succinate (A,N)	+	−	−	−	−	−
Phosphoenol–pyruvate (A,N)	−	−	−	−	−	+
2-Phosphoglycerate (A,N)	+	−	−	−	−	+
Bicarbonate (A)	−	+	−	−	−	−
Bicarbonate (N)	+	+	−	−	−	−
Serine (A)	−	+	−	−	+	−
Serine (N)	−	−	−	−	+	−
KIV (A)	−	−	−	−	+	−
KMV (A)	−	−	−	−	+	−
3-Phosphohydroxypyruvate (A)	−	−	+	−	−	−
Acetamidobutanal (A)	−	−	−	+	−	−
S-3-hydroxy-isobutyrate (A)	−	+	−	−	−	−
Propionyl-CoA (A)	−	+	−	−	−	−
Cholesterol (A)	−	−	−	+	−	−
Desmosterol (A)	−	−	−	+	−	−
Dimethylallyl_diphosphate (A)	+	−	−	−	−	−
Lanosterol (A)	−	−	−	+	−	−
Lathosterol (A)	−	−	−	+	−	−
24-25-Dihydrolanosterol (A)	−	−	−	+	−	−
5Alpha-cholesta-7-24-dien-3beta-ol (A)	−	−	−	+	−	−
Fatty acid (A,N)	+	−	−	−	−	+
Palmitoyl-CoA (A,N)	−	−	−	−	+	−
1-Acyl-sn-glycerol-3-phosphate (A,N)	+	−	−	−	−	+
Myo-inositol-(1-3-4)-trisphosphate (A,N)	−	−	−	−	+	−
Myo-inositol-(1-3-4-5)-tetrakisphosphate (A,N)	−	−	+	+	+	−
Myo-inositol-(1-4)-bisphosphate (A,N)	−	−	−	−	+	−
Myo-inositol-(1-4-5)-trisphosphate (A,N)	−	−	+	+	+	−
Myo-inositol-(4)-monophosphate (A,N)	−	−	−	−	+	−
Phosphatidyl-1D-myo-inositol-4-5-bisphosphate (A,N)	+	−	−	−	+	−
Oxidized glutathione (A,N)	−	−	−	−	−	+
Reduced glutathione (A,N)	−	−	−	−	−	+
5-Amino-levulinate (A,N)	+	−	−	−	−	−

Table 3

Reporter pathway analysis results in different neurodegenerative diseases. “+” sign states significant changes with p-values less than 0.05 (reporter pathways), and “−” sign represents no significant change.

Pathway	PD	AD	ALS	HD	MS	SCH
Glycolysis	+	+	+	+	−	+
Pentose phosphate pathway	−	−	+	−	−	−
TCA cycle	+	+	−	+	+	+
Oxidative phosphorylation and ATPase	+	+	+	+	+	+
Glutamate–glutamine cycle	+	−	−	+	+	+
Ketone body metabolism	+	−	+	+	−	+
Creatine metabolism	−	+	+	−	+	+
Purine nucleoside metabolism	+	+	+	+	+	+
Pyrimidine nucleoside metabolism	+	+	−	−	+	−
Gaba cycle	+	+	+	+	−	−
Aspartate metabolism	−	−	−	−	+	−
Asparagine metabolism	−	+	−	+	+	+
Alanine metabolism	+	−	+	+	−	+
Glycine–serine metabolism	+	−	+	−	−	−
Leucine metabolism	+	−	+	+	−	+
Valine metabolism	+	−	+	+	+	−
Isoleucine metabolism	+	−	+	+	+	−
Methionine metabolism	+	−	−	−	−	−
Fatty acid synthesis	−	+	−	+	−	−
CDP-diacylglycerol biosynthesis	−	+	−	−	−	−
Cholesterol synthesis	−	−	−	−	+	−
Phosphatidylethanolamine metabolism	−	−	−	+	−	−
Phosphatidylcholine metabolism	−	+	−	+	+	−
Sphingomyelin metabolism	+	−	−	−	−	−
Inositol metabolism	+	+	+	+	+	−
Reactive oxygen species pathway	−	+	−	+	+	+

Reporter metabolite analysis (RMA)

The developed brain model was integrated with the transcriptome data (GEO: GSE26927) [18] for a comprehensive investigation of neurodegenerative diseases of AD, PD, HD, ALS, MS and schizophrenia. RMA was applied to the data as described in Methods section, and the reporter metabolites specific to each disease were tabulated comparatively (Table 2). ATP and ADP in neurons are found to be reporter metabolites in all investigated neurological diseases together with Hc, the proton in oxidative phosphorylation. When the shared reporter metabolites for Parkinson’s disease, Alzheimer’s disease, and amyotrophic lateral sclerosis are closely compared (Table 2), most of them are found in oxidative phosphorylation metabolism (ubiquinol, ubiquinone, cytochrome C). Several inositol derivatives are reporter metabolites in Parkinson’s disease.

As for disease-specific reporter metabolites, ketoisovalerate (KIV), alpha-keto-beta-methylvalerate (KMV), palmitoyl-CoA, and malate are found to be reporter metabolites only in PD. Significant transcriptional changes around succinate, 5-amino-levulinate and dimethylallyl-diphosphate are specific only to AD. Succinate plays a role in the citric acid cycle, an energy-yielding process and the latter two metabolites are known to take part in the biosynthesis of secondary metabolites. S-3-hydroxy-isobutyrate, bicarbonate and propionyl-CoA satisfy the reporter metabolites criteria only in ALS. The reporter metabolite specific to HD is 3-phosphohydroxypyruvate, which is involved in serine biosynthesis. The transcriptional changes significant only in multiple sclerosis are around lanosterol, lathosterol, 5alpha-cholesta-7-24-dien-3beta-ol, desmosterol, cholesterol, 24-25-dihydrolanosterol, and acetamidobutanal. Fumarate, oxidized and reduced glutathione, and phosphoenol-pyruvate are reporter metabolites solely in schizophrenia. The disease-specific reporter metabolites are potential candidates as biomarkers for the identification of the corresponding diseases.

The number of common reporter metabolites for these 6 diseases is shown in Fig. 1, where the diagonal cells represent a total number of reporter metabolites for the corresponding neurological diseases. Parkinson’s disease has the highest number of reporter metabolites (40) compared to the other neurodegenerative diseases. The lowest number of the reporter metabolites is found in amyotrophic lateral sclerosis, indicating a relatively unaffected metabolism in this disease.

Fig. 1

Distributions of the common reporter metabolites for neurodegenerative diseases.

The highest number of common reporter metabolites is 16 in the investigated diseases and found to be between Parkinson’s disease and Alzheimer’s disease. In general, schizophrenia has the lowest number of shared reporter metabolites with the other diseases. This indicates that transcriptional changes taking place in schizophrenia are quite different from the others. In fact, that result is expected as schizophrenia is also known to be a psychiatric disorder. The comparison of the common reporter metabolites in different neurodegenerative diseases (Table 2, Fig. 1) elucidates the effect of the transcriptional changes on brain metabolism due to different kind of perturbations (diseases).

Fig. 2 compares the diseases (observations) and the metabolites (variables) in a PCA bi-plot. PCA bi-plots are useful in visualizing the relationship between observations and variables. The figure demonstrates that ALS, SCHZ and AD exhibit relatively similar phenotypes in terms of reporter metabolites whereas MS and PD exhibit distinct profiles. Additionally, metabolites which are scattered around a disease indicate a specific role for these metabolites for the related disease. For example, lathosterol, cholesterol, and desmosterol populate around MS. These metabolites are specifically significantly affected from the perturbation caused by MS, as also revealed in Table 2. On the other hand, these metabolites do not have any reporter role in diseases such as AD, ALS and SCHZ since they are located on the opposite side of the graph.

Fig. 2

PCA bi-plot comparing all investigated diseases in this study. The plot is based on the Z-scores of metabolites. The relationship between the diseases and the effect of metabolites on the diseases are visualized.

Reporter pathway analysis (RPA)

Reporter pathway analysis was performed based on the results of RMA to explicate the pathways significantly affected in the perturbed (disease) state. The pathways having a p-value of less than 0.05 are defined as reporter pathways and they involve significant change due to disease condition. A comparison of the reporter pathways in Table 3 indicates that fatty acid synthesis is regulated in both Alzheimer’s disease and Huntington’s disease significantly. Phosphatidylcholine metabolism is affected in Alzheimer’s disease, Huntington’s disease, and multiple sclerosis. Creatine metabolism is significantly perturbed in Alzheimer’s disease, amyotrophic lateral sclerosis, multiple sclerosis, and schizophrenia (Table 3). Disease-specific reporter pathways are as follows: aspartate and cholesterol metabolisms for MS, phosphatidylethanolamine metabolism for HD, sphingomyelin metabolism for PD, CDP-diacylglycerol biosynthesis for AD and pentose phosphate pathway for ALS. There is no unique reporter pathway for schizophrenia.

Fig. 3 compares common reporter pathways for the investigated neurological diseases, offering a more systematic demonstration of their effect on cell metabolism. Huntington’s disease has the highest number of reporter pathways (17) whereas schizophrenia has the lowest number of the reporter pathways (11). The highest number of common reporter pathways (12) is found between Parkinson’s disease and Huntington’s disease. Although the number of the reporter metabolites in Parkinson’s disease is considerably higher than that in Huntington’s disease, total number of reporter pathways in Huntington’s disease is higher than that in Parkinson’s disease. Similarly, despite having the highest number of common reporter metabolites, Parkinson’s disease and Alzheimer’s disease have the lowest number of common reporter pathways. These cases show that the averaged effect of metabolites over pathways is much more significant than the individual effects.

Fig. 3

Distributions of the common reporter pathways for neurodegenerative diseases.

Computational binding site analysis

Potential transcription factors for the studied diseases were determined based on binding sites of genes neighboring the reporter metabolites, as described in Methods. Identified over-represented motifs for each reporter metabolite for the up- and down-regulated genes are given in detail in Supplementary File 3, together with a lumped version of all identified motifs irrespective of corresponding reporter metabolite for each neurodegenerative disease. According to the results, a number of transcription factors common in more than one disease was pinpointed, among which are paired box (Pax), forkhead box (FOX) and Sox protein families as well as proteins such as Krüppel-like factor (KLF) 4, Upstream stimulatory factor (USF1) 1, and SP1.

Discussion

Energy metabolism

According to the RMA and RPA results, energy-related metabolisms including glycolysis, TCA cycle, oxidative phosphorylation and ATPase, ketone body, and creatine metabolisms are significantly regulated in the neurodegenerative diseases.

Experimental studies on Parkinson’s disease showed that oxidative phosphorylation does not function properly due to the reduced activity of NADH-ubiquinone reductase (Complex I) and NADH cytochrome c reductase [19-21]. Indeed NADH was detected to be a reporter metabolite for PD in our analysis (Table 2). Mitochondrial studies of multiple sclerosis similarly revealed decreased activities of respiratory chain complexes I and III [22]. Succinate, the AD-specific reporter metabolite (Table 2), was reported as a potential biomarker for this disease [23]. The alterations in TCA cycle observed by our computational analysis for PD and AD (Tables 2 and 3) are in agreement with the studies showing a decrease in the α-ketoglutarate dehydrogenase complex enzyme due to Parkinson’s disease [24,25], as well as with the findings reporting changes in a number of mitochondrial enzymes that function in TCA cycle on Alzheimer’s disease [26-28]. Significant perturbations pinpointed by RPA in schizophrenia are supported by studies demonstrating dysregulation of several TCA cycle enzymes and energy metabolism [4,29,30]. Our results on HD are in agreement with the fact that energy metabolism in Huntington’s disease was altered because of dysfunction of oxidative phosphorylation pathway and TCA cycle [31,32]. The reported change in lactate-to-pyruvate ratio [33] and reduced pyruvate dehydrogenase activity [31] in HD are also compatible with the fact that neuronal pyruvate is a reporter metabolite for this disease. Regulation of energy metabolism through altered glycolysis, as predicted for all the diseases except MS (Table 3), can be explained by the role reported for Paired box (Pax) 6 protein, since its mutation leads to abnormal glucose metabolism [34]. Indeed, this transcription factor was predicted to have a role in HD, PD and schizophrenia in our analysis (Supplementary File 3). The potential role of pax family in neurodegenerative diseases is also stated in literature [35].

Another pathway related to neurodegenerative diseases is reactive oxygen species, which is an indicator of oxidative stress in the cell. This pathway was identified to be reporter in AD, HD, MS, and schizophrenia. It is widely known that neurodegenerative diseases are associated with oxidative stress [36,37]. Interestingly, PD and ALS seem to be not regulated at the transcriptional level in terms of this pathway according to our results although the opposite was reported in literature. Our binding-site analysis, however, revealed FOXO3 as a regulating factor for PD and AD (Supplementary File 3), which is in accordance with the role of this protein in oxidative stress [38,39]. Therefore, we can state that the interrelation of PD with oxidative stress was not on the pathway level but through the controlling transcription factors. Indeed, a mechanism was suggested in PD for the link between FOXO3 and FOXF2, another identified motif in our analysis, and neurogenesis [40]. A number of sox family proteins were identified as potential regulators in PD (Supplementary File 3). Among those, the dysregulation of Sox2, which was also predicted as a regulator in schizophrenia, led to oxidative stress and, hence, neurodegeneration [41,42]. Krüppel-like factor (KLF) 4, which is identified in our analysis in all the investigated diseases except Schizophrenia (Supplementary File 3), have recently been assigned a role in Parkinson’s disease through oxidative stress [43].

Creatine helps to supply energy to all cells in the body, and deficiencies in the creatine biosynthetic pathway lead to the malfunction of brain and various severe neurological defects [44]. Our reporter pathway analysis points out alterations in creatine metabolism in different neurodegenerative diseases (AD, ALS, schizophrenia), in agreement with several literature reports [45-48]. There are some studies pointing to creatine supplementation as a potential therapeutic in neurological diseases for energy homeostasis [49-51]. Similar to creatine, ketone bodies, metabolites from another reporter pathway for some of the investigated diseases, are also used in therapeutic strategies to provide energy [52].

Amino acid metabolism

Reporter pathway analysis indicates that amino acid metabolism is affected in the case of neurodegenerative diseases. Methionine metabolism is regulated only in Parkinson disease. In line with this finding, plasma analysis of L-Dopa treated Parkinson’s patients demonstrates a decreased level of methionine and an increased level of homocysteine in methionine metabolism [53]. Significant perturbation in glutamate–glutamine cycle was found by RPA for PD, HD, MS and schizophrenia. The experimental study performed in cerebrospinal fluid stated that glutamate is declined and glutamine level is raised notably in PD [54]. RPA identified alterations in alanine, valine, leucine, and isoleucine metabolisms due to Parkinson’s disease in agreement with cerebrospinal fluid studies which also revealed changes in these amino acids, among others [55]. Similarly, in the case of Huntington’s disease, results of our reporter pathway and reporter metabolite analyses are in accord with the literature findings indicating increased levels of alanine, isoleucine, and glutamate [56,57]. Alteration in valine metabolism in multiple sclerosis observed by RPA is also supported by another experimental study [58]. This study also reports a change in tryptophan and some others, but these are not captured by our analysis. In schizophrenia, significant perturbations observed in leucine and glutamate metabolisms by RPA can be linked to the literature studies reporting alterations in the levels of these amino acids, among others [57,59].

All the above mentioned findings indicate that amino acids, the precursors of neurotransmitters, play key roles in brain. There are abundant evidences in literature such that complex interactions among various amino acids and metabolites in brain support proper functioning of metabolism, leading to a healthy state. Disturbances in the levels and balance among them cause cognitive impairments associated with many psychiatric disorders, such as schizophrenia, and with neurological disorders, such as Parkinson’s, Alzheimer’s disease and many others.

Lipid and related metabolisms

The studied neurodegenerative disorders were found to cause significant perturbations in lipid metabolism (Tables 2 and 3). Having an important role in nourishing brain cells, inositol metabolism is predicted to be affected in PD, AD, HD, ALS and MS, except schizophrenia, whereas sphingomyelin metabolism is only affected in Parkinson disease. In PD, sphingolipid structures are altered due to interactions with alpha-synuclein protein [60]. In Alzheimer’s disease, fatty acid synthesis, phosphatidylcholine metabolism, CDP-diacylglycerol biosynthesis, and inositol metabolism are determined as reporter pathways by RPA. Inositol metabolites are known to play a major role in transducing signals in the nervous system [61,62]. Identification of inositol metabolism as a reporter pathway in AD is in line with literature finding stating a changed level of inositol-related compounds [63,64]. In amyotrophic lateral sclerosis, inositol metabolism is the only reporter pathway among the other lipid-related metabolisms. This result is in accordance with experimental study showing perturbed composition of inositol due to this disease [65]. Reported link between cholesterol metabolism and AD and PD in literature could not be captured by our analysis. Cholesterol synthesis is influenced only in multiple sclerosis according to the RPA results. This finding is in agreement with a study, which demonstrates a decline in cholesterol, lathosterol, desmosterol and lanosterol in multiple sclerosis [66]. The identification of cholesterol precursors as MS-specific reporter metabolites points to their potential as biomarkers for MS pathology, as suggested by a recent experimental study [67]. Although there are similarities, the effects of the different neurodegenerative diseases are different in pathway level in the brain.

Many of the reactions in lipid metabolism are associated with ATP, and cofactors, NADH and NADPH. ATP is a reporter metabolite in all investigated diseases, and the cofactor NADH is a reporter in HD and PD while NADPH is a reporter in MS and schizophrenia (Table 3). The identification of hubs (high connection nodes) in the metabolic network as reporters was also stated before [10]. The binding site analysis of neighbor genes of these hub reporters predicted many common putative regulators, some of which were discussed in previous sections. SP1 and USF1, which are known to regulate the genes of lipid metabolism, are also among them (Supplementary File 3). USF1 is a ubiquitously expressed upstream transcription factor with multiple roles including the regulation of apolipoprotein E [68], an AD-risk gene and known to transport cholesterol from astrocytes to neurons [69]. SP1 belongs to zinc-finger protein family with widespread functions, among which is a reported essential role in the regulation of cholesterol metabolism [70]. The transcriptional activity of SP1 is known to change in Huntington’s disease with neuroprotective role [71,72], in Alzheimer’s disease, [73,74] and in multiple sclerosis [75]. Similarly, literature findings associate the upstream transcription factor 1 (USF1) with AD [68], and MS [76].

Concluding remarks

A comprehensive stoichiometric model of brain metabolism covering interactions in and between astrocytes and neurons was reconstructed via an extensive literature survey. The model, termed iMS570, has the associated gene information for each included reaction whenever possible, with much higher gene coverage compared to the only other similar brain-specific model. After it was shown that the model correctly represents healthy resting brain metabolism in terms of the distribution of major fluxes, we used it as a scaffold for a detailed and systematic analysis of the metabolic effects of six neurodegenerative diseases comparatively. The identified reporter metabolites revealed parts of brain metabolic network perturbed by these diseases at transcriptional level. Such metabolites respond to the perturbations by a change in the expression levels of the surrounding genes, either to adjust their level into a new level or to keep their level constant. Therefore, they are potential biomarkers for the pathology of the related diseases, and they can be used as potential drug targets. Major effects around energy metabolism, reactive oxygen species and lipid metabolism were explained by common transcription factors regulating the genes neighboring the reporter metabolites. Furthermore, the identified transcription factors had roles on neurodegeneration, strengthening the potential functioning of the reporter metabolites as regulatory hot spots. The results, which are presented in the form of reporter metabolites, reporter pathways, and putative transcription factors, present a detailed catalogue of similarities and distinctive properties of the studied diseases at the level of metabolism, and can therefore be used as basis to develop and test several hypotheses. A close clinical inspection is needed to further elaborate on the differences in the affected metabolic pathways observed for these neurodegenerative diseases.

Methods

Transcriptome data

Transcriptome data belonging to six neurodegenerative diseases (AD, PD, HD, ALS, MS and Schizophrenia) were downloaded from Gene Expression Omnibus (Dataset series number: GEO: GSE26927). The dataset was originally generated by Prof. Richard Reynolds and co-workers [18] by using Illumina BeadChip platform for the analysis of common neuroinflammatory pathways in these diseases. Disease samples are from different regions of the brain since the most affected part of the brain is different for each disease. The downloaded dataset including 118 samples was subjected to outlier analysis by using Sammon mapping [77,78] within MATLAB 2012b. Sammon mapping is a dimension reduction method used for high dimensional data and it helps the identification of outliers. In Sammon mapping, distance matrix between the arrays is used in place of original data matrix. A new representation of the arrays in a lower-dimensional space is aimed by minimizing the total stress with reference to the original one. This method is performed on MATLAB with the commands of ‘pdist’ giving Euclidean distances between each pair of observations in the data matrix, and ‘mdscale’ performing non-metric multidimensional scaling on the Euclidean distances. The results of Sammon mapping are plotted to visualize samples among replicates which show a distinct behavior from the rest. In this way, outlier samples are determined. The identified outlier samples were discarded from the dataset before further analysis. One such outlier sample for HD and MS disease sets were identified (GSE: GSM663064 for HD, GSE: GSM663079, respectively), and two outlier samples for schizophrenia healthy set were identified (GSE: GSM663106 and GSE: GSM663108).

Model development

A novel brain metabolic model is developed based on a previous study [6] by (i) expanding the lumped reactions in the original model and (ii) adding new reactions that existed in the literature. Transforming lumped reactions into individual reactions led to a better representation of metabolic pathways. For example, the number of reactions in lipid metabolism is increased about fivefold after unlumping. It is aimed to integrate transcriptome data with the brain reaction model, and gene information controlling each reaction is obtained from HumanCyc, (www.humancyc.org) [11]. New reactions were appended to the model in four steps: (i) Addition of reactions missing in the original model following a detailed investigation of central carbon metabolism, (ii) Addition of pathways for essential and non-essential amino-acids that were not covered in the original model, (iii) Addition of pathways which are reported to have highly expressed controlling genes based on the transcriptomic data of healthy brain, (iv) Addition of other pathways which are reported to be important for brain metabolism. The ultimate goal is to construct a metabolic model of brain with a correct representation of healthy metabolism. The steps followed from model development stage to the final analyses are summarized in a flowchart in Fig. 4. Supplementary File 1, which lists all reactions in iMS570, also provides a complete overview of the model expansion by giving the reference for each reaction in the model.

Fig. 4

Flowchart depicting the model development steps for iMS570 and further analyses.

Central carbon metabolism

A literature survey yielded few reactions related to the central carbon metabolism which were not accounted for in the original model [6], but included in the current model. In glycolysis, 3-phosphoglycerate is also produced from 1-3-biphosphoglycerate via 2-3-Disphospho-D-glycerate without yielding ATP in both astrocytes and neurons (r15, r16, r59, r60) [79]. In addition to the already considered NADPH dependence, NADH is also used in the production of glutamate from alpha–ketoglutarate in astrocytes and neurons (r90, r93) [80]. Phosphoenolpyruvate carboxykinase, transforming oxaloacetate to phosphoenol pyruvate, is expressed in both astrocytes and neurons (r14, r58) [81,82]. In addition to isocitrate biosynthesis, citrate is also converted to oxaloacetate and acetyl-CoA in an ATP-dependent manner in neurons as well as astrocytes enzymatically (r28, r72) [83]. The conversion of citrate to oxaloacetate and acetate is taken into account in astrocytes part of the present brain model (r29) because citrate lyase subunit beta-like protein, responsible for this reaction, was reported to be an astrocyte-specific enzyme [84].

Amino-acid pathways

Pathways for asparagine, histidine, methionine, threonine, arginine, and proline amino acids were added into the present model following an extensive literature survey. All other amino acids were already included in the original model [6].

Asparagine is a non-essential amino acid, and its biosynthesis is taken into account only in neurons, since asparagine synthetase, responsible for enzymatic biosynthesis of glutamate and asparagine from glutamine and aspartate is found to be neuron-specific (r106) [84]. Histidine, an essential amino acid, is taken up by the brain [85] and converted into neurotransmitter histamine in neurons by histidine decarboxylase (r107) [86,87]. Neurons are also capable of synthesizing histamine N-methyltransferase, producing S-adenosyl-L-homocysteine and methylhistamine from histamine and S-adenosyl-L-methionine (r108) [88]. Methionine, an essential amino acid, is taken up by the brain [85]. It is converted into cysteine via the transsulfuration pathway in astrocytes and neurons (r181–r186, r189–r194) [89,90]. In this pathway, homocysteine is formed by three reactions in series and reacts with serine to produce cystathionine which is transformed into cysteine, alpha-ketobutyrate (2-oxobutanoate), and ammonia. Homocysteine is also recycled to methionine due to the activity of 5-methyltetrahydrofolate-homocysteine methyltransferase enzyme in the brain (r187, r195) [91]. The conversion of alpha-ketobutyrate to propionyl-CoA is thought to be active in astrocytes due to the control of this reaction by astrocyte-specific branched-chain α-keto acid dehydrogenase complex (r188) [84]. Threonine is an essential amino acid and can cross the blood-brain barrier with a slow rate compared to other amino acids [85]. Its metabolism contains only two enzymatic reactions in the brain cells (r196–r201) [92]. First, threonine is converted to 2-amino-3-ketobutyrate by threonine 3-dehydrogenase. Then, 2-amino-3-ketobutyrate is catabolized to glycine and acetyl-CoA via 2-amino-3-ketobutyrate ligase. Arginine is taken up by the astrocytes (r607) [93]. It is used to produce agmatine (r446) [93,94] and ornithine [95] (r447) in astrocytes. Arginine is also transferred to neurons by astrocytes to synthesize agmatine (r453) [96], citrulline (r455) [97], and ornithine (r458) [95]. In addition to its uptake, arginine is also produced in neurons from aspartate and citrulline (r456–r457) [95,98,99]. Ornithine is taken up by both cell types, besides its production in astrocytes and neurons [100]. In the cytosol of astrocytes, pyrroline-5-carboxylate, an intermediate metabolite in glutamate synthesis from ornithine, and NADPH are degraded into proline and NADP+ (r179). In astrocytic mitochondria, proline and FAD are converted to pyrroline-5-carboxylate and FADH2 (r180). These two proline producing and consuming reactions are considered in only astrocytes in the brain model, because Pyrroline-5-carboxylate reductase-like protein and Proline dehydrogenase were stated to be astrocyte-specific [84].

Active pathways based on transcriptome data

Pathways controlled by the highly expressed genes were identified based on the transcriptome data of a healthy brain. To this aim, gene expression data of different brain regions for healthy individuals were used (control samples of GEO: GSE26927 dataset [18], storing the transcriptome data for 6 common neurodegenerative diseases). Since the aim of the present study is to develop a general model for the whole brain, the brain model constructed here is not region-specific. The identification of metabolic genes with high expression was achieved in 5 steps: (i) Genes with expression values higher than 3-fold of the average expression values of all genes in the same transcriptome data were extracted. (ii) Among these genes, those which are expressed in only one of the investigated regions were discarded to minimize false positives. (iii) Of the remaining genes, the ones associated with cellular metabolism were identified by matching with reaction-controlling genes from HumanCyc. (iv) Metabolic pathways corresponding to these active genes were extracted using the information given by HumanCyc when available. If a pathway includes many reactions, but genes controlling only few of its reactions are found to be active in this pathway according to transcriptome data analysis, the related pathway was disregarded since such genes are active probably because of different functions other than metabolic control. (v) The candidate pathways were added into the present brain model if further support from literature was found as to the activity of the pathway in the brain. This procedure led to the involvement of following pathways into the brain model: inositol, creatine, polyamine, heme, pyrimidine and purine nucleosides.

Inositol is taken up by the brain [101]. It is used to produce diacylglycerol, complex inositol phosphate, and phosphatidylinositol (r388–r415) [102,103]. Another important metabolism consuming arginine and functioning in both neurons and astrocytes is creatine metabolism (r474–r479) consisting of ornithine and guanidinoacetate production from glycine and arginine [104]. Guanidinoacetate is then transformed to creatine to be used in creatine phosphate synthesis (r475–r1476, r478–r479) [105]. Agmatine synthesized in neurons from arginine is converted into putrescine (r454) [106,107] to be utilized both in polyamine metabolism including interconversion of precursor putrescine, spermidine, and spermine in neurons (r467–r473) [108] and in GABA synthesis from putrescine in astrocytes (r461–r464) [109,110]. The production of putrescine from ornithine is neuron-specific (r459) [111] and glutamate production from ornithine and alpha-ketoglutarate is astrocyte-specific (r448–r450) [112]. The metabolism of heme, an iron containing chemical compound, is included in the present study in both neurons and astrocytes since it has active genes satisfying our criteria and supported by literature [113,114]. Heme is synthesized from glycine and succinyl-CoA by eight reactions in series and degraded to bilirubin via biliverdin intermediate enzymatically (r480–r499) [115]. Adenosine, guanosine, uridine, and cytidine are nucleosides encountered in brain metabolism and they are mainly synthesized in liver and transported into the brain through the blood-brain barrier [116]. These nucleosides, in addition to a five-carbon sugar, contain the bases adenine, guanine, uracil and cytosine, respectively. These nucleosides are transformed into related nucleotides in the brain, which are utilized in metabolic processes or recycled to their nucleosides (r500–r571) [116].

Addition of other pathways

Pathways for ketone body metabolism, cardiolipin and sphingomyelin were included in the present model due to their importance for brain metabolism.

Acetoacetate and 3-hydroxybutyrate are the most important ketone bodies in human brain energy metabolism. They are produced by liver and transported by blood to brain where they can be used in ATP production as an alternative to glucose [117]. The uptake rates of ketone bodies from blood-brain barrier into brain are very small in comparison with glucose in resting condition; however, they significantly change in case of starvation in which uptake rates of ketone bodies increase with a decline in glucose uptake rate [118]. Acetoacetate and 3-hydroxybutyrate are converted into acetyl-CoA (r438–r444, r127–r128, r164) which participates in TCA cycle reactions or other metabolic reactions both in astrocytes and neurons [119]. Cardiolipin is a phospholipid and found in the mitochondrial membrane [120]. It is synthesized from CDP-diacylglycerol and Glycerol-3-phosphate by three enzymatic reactions in series in the brain (r366–r371) [121]. Sphingomyelin, a type of sphingolipid, is localized in the brain [120,122]. Precursors of sphingomyelin synthesis are serine and palmitoyl-CoA. Astrocytes use its own serine; however, neurons utilize serine transported from astrocytes in sphingomyelin synthesis since there is no serine synthesis in neurons (r372–r376, r378–r382) [123-125]. Sphingomyelin biosynthesis includes five metabolic reactions with 3-dehydrosphinganine, sphinganine, dihydroceramide, and ceramide as the intermediate metabolites. Sphingomyelin is also converted to ceramide and phosphoryl-choline due to the availability of sphingomyelinase in the brain (r377, r383) [126].

Flux balance analysis

The model was validated by applying flux balance analysis [12] for the prediction of major flux splits in the resting (healthy) state. Simulations were performed in MATLAB 2012b environment. Optimization (linear programming and quadratic programming) was performed using MATLAB bindings of CPLEX ILOG optimization algorithms provided by IBM’s Academic Initiative Program. The healthy-case FBA simulation was performed using the exact constraints as reported for the original model [6], as well as few other following constraints. The uptake rates for newly added amino acids were constrained as 0.0025 [85], 0.0017 [85], 0.0008 [85], 0.0020 [127], and 0.0031 [85] in units of μmole/g tissue/min for histidine, methionine, threonine, arginine and ornithine, respectively. Uptake rates of ketone bodies acetoacetate and 3-hydroxybutryrate were constrained to be zero since this pathway is known to be inactive in resting state [118].

Reporter metabolite analysis

Reporter metabolites are defined as metabolites around which the most significant transcriptional changes occur in response to a genetic or environmental perturbation, including diseases [10,128]. In this way, alterations in a metabolic network are elucidated at the transcriptional level. For this, genes are assigned p-values based on the differential transcriptome data using t-test. p-value of each gene (p), representing the significance of differential gene data, is transformed into Z score (Z) by using the inverse normal cumulative distribution. Then, each metabolite in the metabolic network in question is scored based on the aggregated Z-scores of the genes controlling the reactions in which the metabolite is involved, i.e., each metabolite is assessed by the Z-scores of these k neighbor genes.

In the case where a reaction is associated with more than one gene, the one with most significant change is taken into account [10]. Then, Z scores are corrected by using mean (μ) and standard deviation (σ) of the aggregated Z scores of many sets of randomly selected k genes from the metabolic network.

In the last step, corrected Z scores are converted to p-values, and metabolites having p-values less than the selected cut-off are named as reporter metabolites. For reporter metabolite analysis (RMA), p-value of 0.05 was selected as cut-off. Those metabolites with only a single neighbor gene were not considered as reporter metabolites in evaluating the results. RMA calculation was performed using the online BioMet Toolbox [129] ttest2 command of Statistics Toolbox of MATLAB was used to calculate p-values based on two-sample t-test.

Reporter pathway analysis

A novel analysis, termed reporter pathway analysis (RPA), was performed to predict significantly affected metabolic pathways due to the investigated diseases with the primary goal of making it easier to analyze RMA results. Here, each pathway is scored based on the p-values of its metabolites, calculated through RMA. The scoring is based on Eqs. (1) and (2), leading to pathway p-values. The pathways having p-values less than the selected cut-off (0.05) are defined as reporter pathways, i.e., pathways significantly affected due to the perturbation. Calculations were performed using BioMet Toolbox [129]. Note that we hereby opted for a global network-centered effect of the diseases on pathways by using metabolite p-values to calculate pathway scores since metabolite p-values include information on the interactions with other pathways.

To investigate whether genes associated with a reporter metabolite are controlled by common transcription factors, computational binding site analysis was performed. For each reporter metabolite, its neighboring genes were grouped into two sets with respect to up- and down-regulation [10]. For the sets with at least 4 gene members, DNA sequence regions with −950 to +50 base pairs were scanned from the JASPAR database [130] for the significantly over-represented motifs of the stored transcription factors (p-value: 0.05). The analysis was performed using the internet-based software tool, Pscan [131].

Principal component analysis (PCA)

Metabolites which did not show a significant change in any of the studied diseases (p-value cut-off: 0.05) were discarded from the set. Z-scores of the remaining metabolites were used in the PCA analysis. The analysis was performed using pca function of Statistical Toolbox of MATLAB.