Mathematical analysis of the influence of brain metabolism on the BOLD signal in Alzheimer's disease.

Blood oxygen level-dependent functional magnetic resonance imaging (BOLD-fMRI) is a standard clinical tool for the detection of brain activation. In Alzheimer's disease (AD), task-related and resting state fMRI have been used to detect brain dysfunction. It has been shown that the shape of the BOLD response is affected in early AD. To correctly interpret these changes, the mechanisms responsible for the observed behaviour need to be known. The parameters of the canonical hemodynamic response function (HRF) commonly used in the analysis of fMRI data have no direct biological interpretation and cannot be used to answer this question. We here present a model that allows relating AD-specific changes in the BOLD shape to changes in the underlying energy metabolism. According to our findings, the classic view that differences in the BOLD shape are only attributed to changes in strength and duration of the stimulus does not hold. Instead, peak height, peak timing and full width at half maximum are sensitive to changes in the reaction rate of several metabolic reactions. Our systems-theoretic approach allows the use of patient-specific clinical data to predict dementia-driven changes in the HRF, which can be used to improve the results of fMRI analyses in AD patients.

Introduction

Blood oxygen level-dependent functional magnetic resonance imaging (BOLD-fMRI) is a standard clinical tool for the detection of brain activation. The physiological response to stimulation is an increased consumption of oxygen coupled with an increase in blood flow. The combination of these two effects leads to local changes in the deoxyhemoglobin concentration in the brain, which results in a shift of the local magnetic susceptibility. This dynamic process can be observed using magnetic resonance imaging.

In Alzheimer’s disease (AD), task-related fMRI has been used to detect early brain dysfunction.[1,2] Already in 2005, Rombouts et al.[3] have shown that peak time and peak height of the BOLD response are affected in early AD.

Several studies with AD patients’ brain tissue[4-6] have reported changes in the expression levels of enzymes important for the regulation of glycolysis and the pentose-phosphate pathway. As the dynamics of oxygen consumption depends on the complicated interplay of those transporters and the enzymes that regulate glycolysis, the pentose phosphate pathway and mitochondrial respiration, the question arises whether the observed changes in enzyme expression in AD can affect the results of BOLD fMRI.

To address this question, we present a mathematical model of brain energy metabolism, neuronal stimulation and the hemodynamic response. The model describes the dynamic sequence from neuronal stimulation to metabolic and vascular response and predicts the effect of changes in enzyme activity, flux distribution or metabolite concentration on the BOLD-shape. Our systems-theoretic approach allows the use of patient-specific clinical data to predict dementia-driven changes in the hemodynamic response function (HRF). We use this model to calculate the sensitivities of peak height, peak time and full width at half maximum for all metabolic reactions explicitly considered in the model. We furthermore simulate different scenarios, including increased expression of key glycolytic enzymes, increased expression of enzymes regulating the flux through the pentose-phosphate pathway and astrocyte hypertrophy. All scenarios are based on data reported in AD patients. The resulting BOLD shapes are compared with the BOLD shapes for the healthy state as predicted by the model. Our goal here is not to model a clinical situation but to delineate the individual contributions of chosen factors to the BOLD signal. Thus, the obtained results may contribute to a better understanding of data gathered in the clinical setting.

Kinetic modelling of brain metabolism

Mathematical modelling of the brain energy metabolism and the hemodynamic response has been used to describe the BOLD signal for more than a decade. To the best of our knowledge, the first work which combined a model of energy metabolism and a description of the neuronal stimulation was published by Aubert et al. 2001.[7] Their model, derived from a simplified representation of erythrocyte glycolysis, includes metabolic and transport reactions as well as a description of neuronal stimulation via dynamic changes in the sodium concentration. Aubert and Costalat[8,9] subsequently refined this model, the most notable steps being the application to functional neuroimaging[8] and the consideration of neuron-astrocytic interaction.9 The final model in this series was published in 2007 and used in vivo and in vitro data to delineate the importance of both neurons and glia cells in BOLD-fMRI studies.[10] Based on these publications, a number of models by other groups analysed different aspects of brain energy-metabolism in more detail. Notable works include the model by Tiveci et al.[11] which analysed the influence of calcium dynamics and the model by Cloutier et al.,[12] which included an explicit description of neurotransmitter cycling and the glycogen storage system.

To link these models with functional neuroimaging, the models in Aubert et al.[8-10] and Tiveci et al.[11] employ a description of the BOLD response based on the balloon model developed and popularized in a series of publications by Richard B. Buxton et al.[13-15] All these models can be considered mechanistic descriptions of the interactions between brain energy metabolism and neuronal stimulation. This is in stark contrast to the HRF models used in current clinical practice, described next.

In clinical functional neuroimaging, the HRF is modelled in a purely descriptive form. The standard shape, as implemented in SPM 12,[16] is often referred to as the canonical HRF. It is a linear combination of two gamma functions Γ(·). We shall here denote the time-dependent HRF as hrf(t):

The names of the parameters are identical to those in SPM 12 and are explained in Table 1.

Table 1.

Parameters of the canonical HRF.

Parameter	Meaning	Default value
p(1)	Delay of response relative to onset	6
p(2)	Delay of undershoot relative to onset	16
p(3)	Dispersion of response	1
p(4)	Dispersion of undershoot	1
p(5)	Ratio of response to undershoot	6
dt	Ratio of scan repeat time and time resolution	–

Note: Parameters of the canonical hemodynamic response function used in SPM 12 (using the same notation).

A similar description of the HRF is available for convolution in the software package FSL[17] under the name double-gamma. There exist several other approaches to model the hemodynamic response, including the use of linear combinations of reference waveforms or the incorporation of the temporal derivative of the canonical HRF.[18] While these approaches can be used to account for a shift in the onset of the BOLD response and model region-specific differences, they do not try to be a representation of the biology. Instead they provide a phenomenological description of the observed data. The shape of the canonical HRF described above, for example, can be adjusted with parameters accounting for a delay of the response and a delay of the undershoot relative to the onset, dispersion of response and undershoot, ratio of response to undershoot and time of onset, but there is no connection of these parameters to biological quantities (see Table 1).

Owing to the descriptive nature of these parameters, there is no straightforward way to include information about a patient’s disease state into phenomenological models of the hemodynamic response. Instead, changes in the HRF that result in changes of the observed BOLD signal are solely attributed to changes in the neuronal stimulation patterns as described next.

Interpretation of variations in the shape of the BOLD signal

In the analysis of BOLD fMRI data, the shape of the BOLD signal is commonly characterized by three parameters: (a) peak height, (b) time-of-peak and (c) full width at half maximum (FWHM). Because the BOLD response is modelled as the convolution of a fixed stimulus function with the HRF, changes in the three characteristic parameters are direct consequences of corresponding changes in the HRF. In the classical interpretation of the HRF shape, changes in its peak height, time-of-peak and FWHM are attributed to changes in firing rate, onset latency and duration of neural activity, respectively.[18] Due to this simplified view of the origin of the HRF, there is no attribution of metabolic influences. In a situation where the metabolic state of different patients in a study can be considered identical, this poses no problem. However, if the metabolic reactions underlying the hemodynamic response are changing, e.g. due to different stages of disease progression, these differences need to be considered in the evaluation of BOLD fMRI data. A prominent example for such a situation is AD, where increases in the activity of different enzymes regulating the pentose phosphate pathway have been observed.[5,6] This finding was supported by Orešič et al.,[19] who identified the upregulation of pentose phosphate pathway in patients progressing to AD. Mechanistic kinetic modelling of brain energy metabolism provides a tool to include this valuable information into the analysis of BOLD fMRI data.

Materials and methods

Description of the model

To analyse the influence of alterations in the metabolic network, we build upon the models by Aubert et al.[9] and, Cloutier et al.[12] and Heinrich and Schuster.[20] As previous modelling approaches restricted the metabolic networks to a description of glycolysis, the TCA-cycle and oxidative phosphorylation only, we extended the model and included a description of the pentose-phosphate pathway (PPP) in neurons and astrocytes.

After inclusion of the PPP, the final model consists of 64 reactions (40 metabolic reactions and 24 transport reactions) and 65 species. Parameter values listed in the following section and the tables in this publication are derived from the literature.[9,12,20]

The model is deposited in the BioModels database and has been assigned the identifier MODEL1603240000.[21]

Compartments and transport reactions

The full model comprises six compartments: Neurons, astrocytes, the extracellular space, capillaries, larger blood vessels (arteries) and veins. The size of the individual compartments is given by their relative amount and is fixed for neurons (Vn = 0.45), astrocytes (Vg = 0.25), the extracellular space (Ve = 0.2), capillaries (Vc = 0.0055) and arteries (Va = 0.0055).[8] The size of the venous compartment is described by an ordinary differential equation derived from the difference between the blood flow into the veins (F) and the blood flow out of the veins (F) with an initial volume of V,0 = 0.02396. At steady state, F is set to a constant value of 0.012 and F is given by with F0 = 0.012, α = 0.5 and τ = 2. This equation takes the viscoelastic properties of the veins into account, τ being the viscosity parameter, which results in a post-stimulus undershoot in the BOLD signal.[15]

Transport reactions govern the exchange of oxygen, glucose, carbon dioxide, lactate, glutamate and sodium between the first five compartments (Figure 1). Transport between the arteries and capillaries is driven by diffusion and convection. For transport between capillaries, astrocytes, neurons and the extracellular space facilitated diffusion is assumed for glucose and lactate. Oxygen exchange between capillaries and neurons as well as between capillaries and astrocytes is modelled as simple diffusion according to Fick’s law. Following the description given in Aubert and Costalat[8] and Vafaee and Gjedde[22], the concentration of intracellular oxygen is coupled to the plasma concentration of oxygen via a Hill relation. All reactions and the respective parameter values are listed in supplementary Table 5.

Figure 1.

Overview of the transport reactions considered in the model. All reactions are reversible. For glucose (Glc), lactate (Lac), oxygen (O2) and carbon dioxide (CO2) the arrows indicate the direction of flux at steady state. For glutamate (Glu) and sodium (Na+) the arrows indicate the flux after stimulation. The notation follows the Systems Biology Graphical Notation (SBGN).[23]

Metabolic pathways

Metabolic reactions for glycolysis, the pentose-phosphate pathway and mitochondrial respiration are modelled in neurons and astrocytes (Figure 2). While glycolysis and mitochondrial respiration has already been described in the publications by Aubert et al.[9] and Cloutier et al.,[12] we applied a more accurate description of the reactions catalysed by hexokinase and phosphoglucoisomerase suggested by Heinrich and Schuster.[20] Additionally, our model includes for the first time a detailed description of the PPP. The arterial concentrations of glucose, oxygen, lactate and carbon dioxide and the extracellular sodium concentration are fixed boundary conditions. All other species concentrations are governed by rate laws (Supplementary Table 4).

Figure 2.

Metabolic processes considered in astrocytes. The same metabolic reactions are considered in neurons. Nodes inside the dashed box present the pentose phosphate pathway, not considered in previous models of brain energy metabolism. In addition to the species already depicted in Figure 1, the following species are shown: glucose 6-phosphate (G6P), fructose 6-phosphate (F6P), glyceraldehyde 3-phosphate (GAP), phosphoenolpyruvate (PEP), pyruvic acid (PYR), creatine (Cr), phosphocreatine (PCr), 6-phospho gluconolactone (G6L), 6-phospho lactonate (P6G), ribulose 5-phosphate (Ru5P), ribose 5-phosphate (R5P), xylulose 5-phosphate (X5P), erythrose 4-phosphate (E4P), sedoheptulose 7-phosphate (S7P), adenosin triphosphate (ATP), adenosin diphosphate (ADP), adenosine monophosphate (AMP), NAD, NADH, NADP and NADPH. The notation follows the Systems Biology Graphical Notation (SBGN).[23]

To include the pentose phosphate pathway, we followed an approach developed by Stanford et al.[24] This approach makes use of the common modular rate law[25] and is based on data from public databases and publications to compute initial guesses for the kinetic parameters. To parameterise the equations of the pentose-phosphate pathway in our model, we collected data from SabioRK,[26] eQuilibrator[27] and public available models where applicable. To achieve a thermodynamically meaningful parameterization based on these data, the collected values needed to be harmonized. Parameter balancing is a statistical method to infer model parameters based on experimental data measured under inconsistent experimental conditions such as temperature and pH value.[28] We used the implementation available on the semanticSBML website with the input values listed in supplementary Table 1.

For species which take part in the glycolytic and the pentose-phosphate branch, the initial concentration was set to the steady state concentration of the respective species in the glycolysis only model. As an initial guess for species concentrations exclusive to the PPP, we used the arithmetic mean of the concentration of the metabolic species already contained in the original model. After parameter balancing, we applied enzyme rescaling to achieve a physiologically meaningful flux distribution between the glycolytic branch and the pentose-phosphate branch. The only exception from this procedure was the calibration of the NADPH oxidase, which constitutes the single source of NADP in the model. As this reaction is a housekeeping reaction to guarantee the availability of NADP, it is basically a proxy for all other NADP sources not considered in the model. To account for this special case, we chose simple mass action kinetics for this reaction and set the rate constant k for NADPH oxidase in neurons and astrocytes to 1 before enzyme rescaling. Enzyme rescaling is described in detail in Kholodenko et al.[29]. The amount of flux through the PPP in healthy steady state is set to about 6%of the total flux, a value taken from Holzhütter,[30] which is close to experimental findings for the flux distribution in adult brain tissue.[31] The full equations for the metabolic model are listed in supplementary Tables 2 and 3. In the absence of stimulation, the model shows a steady state behaviour with species concentrations in line with the values described in the data available.

Stimulation

Neuronal stimulation is described as a combination of (a) release of glutamate from neurons into the extracellular space coupled with an influx of sodium from the extracellular space to the neurons, (b) an increase in the cerebral blood flow. The description of the increase in sodium influx follows an approach taken by Aubert and Costalat[8] The stimulation is described as the sum of a constant term and a gamma function where v1 = 0.041, v2 = 1.44, t0 = 200 and τ = 2.[12] This setup is in line with the notion of habituation of neuronal activity and can be mapped to time courses of local field potential published by Logothetis et al.[32] A time course of the sodium influx is depicted in Figure 3, panel (a). The increase in neuronal sodium is coupled to the increase in extracellular glutamate which initiates the neurotransmitter cycling (see supplementary Table 5). where Δ = 0.42, a = − 4.59186, t 0 = 200, t 1 = 2, t end = 40 and δ = 3, which is in line with the description in Cloutier et al.[12] The increased blood flow influences the size of the venous compartment via F in = F0 +  (equation (2)). The resulting venous volume is depicted in Figure 3, panel (b).

Figure 3.

Time course of the two dynamic inputs into the model and the resulting BOLD response. (a) Sodium influx as a surrogate for neuronal stimulation (v stim). (b) Change in venous volume as a result of the increased blood flow (V v). (c) BOLD response evoked from the interplay of metabolic and hemodynamic response.

The increase in blood flow is described as a bistable switch

The increase in neuronal sodium as a result of the stimulation and the increase in astrocytic sodium in response to glutamate uptake drive the sodium exchange system away from their steady state solution. This leads to an activation of the sodium pumps which transport sodium from neurons and astrocytes towards the extracellular space, thereby consuming ATP. The ATP pool is subsequently replenished by an increased activity of the metabolic network, which leads to concentration changes of capillary oxygen, visible in the dynamic profile of dHb (Supplementary Figure 1).

The time evolution of the dHb concentration is governed by the following equation which describes its dependency on arterial and capillary oxygen concentration where F describing the flow into the venous compartment, F the outflow of the venous compartment and V the venous volume.

BOLD response

The BOLD response is calculated in dependence of the venous volume and the change in deoxygenated hemoglobin (dHb) with the following equation which is taken from Friston et al.[33] This equation has been used in a mathematical model of brain energy metabolism before.[7] The resulting shape (Figure 3(c)) compares well to the shape reported in other modelling studies.[7,9,10]

Model analysis

The purpose of our model is the analysis of the effect of changes in enzyme activity and transporter activity on the shape of the BOLD-response. Although we are aware that a number of different metabolic changes has been implicated to AD, we concentrate our analysis on the following published experimental findings:

The activity of the pentose-phosphate pathway is increased in AD

In the publication by Russel et al.[6] an upregulation of neuronal glucose-6-phosphate dehydrogenase was observed. The authors concluded that the observed increase in PPP activity is consistent with an ‘attempted reductive compensation to oxidative stress’. In the paper by Palmer,[5] a significant increase in glucose-6-phosphate dehydrogenase and 6-phosphogluconate dehydrogenase in the inferior temporal cortex of AD patients is described. The authors use this observation to suggest a ‘state of oxidative stress in the AD brain’. The question whether the upregulation of the PPP is actually a result of increased oxidative damage is outside the scope of our analysis.

The activity of pyruvate kinase, lactate dehydrogenase and phosphofructokinase is increased in several brain regions in AD

Bigl et al.[4,33] published a series of articles reporting increased levels of phosphofructokinase (PFK) pyruvate kinase (PK) and lactate dehydrogenase. An increased activity of PFK in the frontal and temporal cortex of AD patients compared to a control group was reported. In the same study, no change in the activity of aldolase was found.[34] In 1999, the same group of authors analysed the activity of other glycolytic enzymes, including hexokinase, PK and lactate dehydrogenase. In this study, significantly increased levels of lactate dehydrogenase in the basal forebrain and the frontal cortex as well as significantly increased levels of PK in the frontal cortex are reported. While the values of enzyme activity are only given as bars inside of figures, the increase can be estimated to be about 20% in both enzymes. For hexokinase, no significant change could be established.[4]

Astrocyte hyperplasia and hypertrophy have been reported in AD patients

Immunochemical studies of activated astrocytes in hippocampus and enthorinal cortex have shown an increased number of GFAP positive astrocytes.[35] Similar results were published by Robinson,[36] who reported a 1.4-fold increase in astrocytic density in AD patients compared to controls, and Vijayan et al.,[37] who found a general increase in the number of astrocytes. Recent evidence suggests, however, that the observed changes in the number of astrocytes are not due to increased proliferation, but are a consequence of phenotypic changes.[38] This would be in line with several reports of astrocytic hypertrophy as a consequence of reactive astrogliosis in AD.[39,40] Regardless of the underlying phenomena, both hyperplasia and hypertrophy result in an increased astrocytic volume fraction in the observed voxel.

As described above, changes in the BOLD-response are typically attributed to changes not directly related to brain energy metabolism. Our goal was to elucidate how the different changes in metabolism and transport described above would affect the parameters of the BOLD response. To answer this question, we undertook a two-step analysis of our kinetic model. In a first step, we use sensitivity analysis to quantify the effect of changes in the models parameters to the three main shape parameters of the BOLD response: time-of-peak, peak height and full width at half maximum. In a second step, we created scenarios that model the situation described above, analysing the impact of several simultaneous changes to the models parameters on the BOLD shape.

Sensitivity analysis of the three shape parameters of the BOLD-response was performed using Copasi V4.18 for t = 0 to t = 215s using 2,150,000 time steps (i.e. Δt = 0.0001).[36] The chosen time frame guarantees that peak time, peak height and full width at half maximum can be calculated from the sample. For the calculation of peak time and FWHM, the default numerical precision of the ODE solver (relTol = 1e-06, absTol = 1e-12) was used, whereas for the calculation of peak time, a precision of the chosen Δt was achieved. Scaled sensitivities were subsequently calculated with Matlab R2011b (The Mathworks, Natick, MA) according to the following formula where x is the cause (e.g. enzyme concentration) and f(x) is the effect (e.g. peak height). For Δx, the default setting of Copasi V4.18 was used (Δx = 0.001 · x). This formula provides a simple method for computing the finite difference approximation of the derivative and has been proven to yield robust results.[37]

Analysing the influence of metabolic alterations on the BOLD shape

To measure the influence of experimentally observed changes in enzyme expression as described above, we set the values for the enzyme concentration as reported in the specific publication. The following scenarios are considered, for which the arguments were provided above:

Healthy state: This is the model with all parameters unchanged.

PPP activation: According to the experimental findings, the enzyme level for ZWF was increased to 162%and the enzyme level of GND to 126%. The rate constant of the NADPH oxidase was also increased by 162% to model the hypothesized increased demand for NADPH in the antioxidant defence system.

Glycolysis activation: The enzyme activity of lactate dehydrogenase (LDH), PK and PFK was increased to 120% of its original value in neurons and astrocytes.

Increased astrocytic volume: To model the increase in the share of brain volume occupied by astrocytes, independently of the underlying mechanism, the relative amount of astrocytes in the model was increased from 0.25 to 0.3, while the relative amount of the extracellular space was reduced from 0.2 to 0.15 of the total volume, so that the total volume remains unchanged.

Each scenario was simulated and the resulting shift in the flux distribution and the value of the three characteristic shape parameters of the BOLD response were calculated.

Results

We describe here the results of the analyses in the same order as above, including sensitivity analysis, the analysis of the different scenarios and of the influence of changes in blood glucose levels.

Sensitivity analysis

Sensitivity analysis of the three main characteristics of the BOLD response was conducted for the model parameters which describe enzyme concentrations. The results, depicted in Figures 4 and 5, indicate that all three characteristics: peak height, full width at half maximum and peak timing are sensitive to changes in the reaction rate of reactions both from the glycolytic and the PPP branch. In the glycolytic branch (Figure 4), the three characteristics are most sensitive to changes in the reactions HK, PGI forward, PGI backward and PFK. The latter three are the glycolytic reactions with direct connections to the pentose-phosphate pathway. Positive sensitivities in the neuronal reactions HK, PGI forward and PFK for peak height and FWHM correspond to negligible sensitivities of peak time in these reactions. In astrocytes, the positive sensitivity of the peak height in HK and PGI forward is coupled to a negative sensitivity in full width at half maximum for these reactions. While for most reactions, the effect of changes in neurons is more pronounced than in astrocytes, the opposite is true for LDH. Full width at half maximum is the parameter with the most pronounced sensitivity to changes in the glycolytic reactions overall. In the pentose-phosphate pathway, the sensitivities of FWHM are again larger than the sensitivities in peak height. For peak time, all non-negative sensitivities are identical, indicating that the difference is identical to the chosen Δt with which the peak timing was calculated. With the exception of RKI in astrocytes and TKL-2 and TAL in neurons, positive sensitivities for peak height correspond to negative sensitivities in FWHM.

Figure 4.

Sensitivities of peak height, full width at half maximum and peak timing of the BOLD signal for glycolytic reactions. All three properties are sensitive to small changes in the reaction rate of HK, PGI and PFK. With the exception of LDH and PFK, changes in the neuronal reactions evoke overall higher responses compared to changes in astrocytic reactions.

Figure 5.

Sensitivities of peak height, full width at half maximum and peak timing of the BOLD signal for reactions of the pentose-phosphate pathway.

Increased activity of ZWF and GND – PPP activation scenario

Altering the enzyme concentration for the reactions ZWF and GND according to the observations in Palmer[5] led only to a slightly decreased BOLD peak height and full width at half maximum: less than 0.05% in peak height and FWHM. This result was expected, as according to prior sensitivity analysis, the sensitivities of the shape parameters to changes in the PPP are minute.

Increased activity of PFK, LDH and PK – Glycolysis activation scenario

Increasing the concentration of PFK, lactate dehydrogenase and PK by 20% lead to an decrease in peak height by 0.84%, an decrease in peak time by 0.07% and a decrease in FWHM by 5.9%. The alteration in the shape of the BOLD response is more pronounced than the alteration in the PPP activation scenario, but the impact is still small, especially for peak height and peak time.

Astrocytic hypertrophy or hyperplasia – Increased astrocytic volume scenario

Increasing the astrocytic volume led to a decrease in the peak height by 2.38%, a reduction in peak time by 0.1% and a decrease in FWHM by 14% (Figure 6).

Figure 6.

Comparison of the BOLD shape arising from the different metabolic scenarios. No visual distinction is possible between the normal (healthy) scenario and the scenarios ‘PPP activation’. The most prominent change can be observed for the ‘astrocytic volume increase’ scenario. The right panel provides a detailed view of the peak.

Discussion

Summarizing the different results described above, we derive the following main conclusions:

Our extended mathematical model allows an investigation of the effect of metabolic alterations reported in AD brains on the BOLD shape. BOLD-fMRI imaging relies on a suitable model of the underlying hemodynamic response. Using a mathematical model of brain energy metabolism, neuronal stimulation and venous volume change, we are able to demonstrate how various reported changes in brain metabolism influence the three important parameters of the BOLD shape. In our model, the BOLD response is not calculated as the convolution of a HRF with a fixed stimulus pattern, but instead directly calculated from a mathematical description of the underlying metabolic network and the hemodynamic response to a simulated stimulation. This allows us to relate the observed changes in the BOLD shape parameters to changes in the dynamics of the metabolic reaction network. An advantage of our approach is that our model allows a simultaneous assessment of the influences on the shape of the BOLD response, such as changes in the flux distribution between the main glycolytic branch and the PPP, changes in expression of glycolytic enzymes and changes resulting from astrocytosis.

The shape of the BOLD response is sensitive to changes in the metabolic network. The sensitivity analysis revealed that the three characteristic properties of the BOLD signal, peak height, peak timing and full width at half maximum are sensitive to changes in the reaction rate of several reactions. While the calculated sensitivities are small, they are not negligible because they appear consistently. The classic view that changes in the BOLD shape are only attributed to changes in strength and duration of the stimulus does not hold. Contrary to what is suggested in Rombouts et al.,[3] the peak timing is least sensitive to changes in the metabolic network. It is instead the FWHM which is most sensitive to these changes.

Reported changes in enzyme activity have only small influence on the BOLD shape. The analysis of the literature-derived scenarios revealed that the reported changes in enzyme activity have only a small influence on the simulated BOLD shape. The PPP activation scenario, which analysed the impact of an increase activity of the pentose-phosphate pathway, resulted in a BOLD shape that can for practical purposes be considered identical to the healthy state. The glycolysis activation scenario, which assumes increased levels of PFK, PK and lactate dehydrogenase, decreased the BOLD shape by less than 1% and reduced the FWHM by about 2%. While these changes seem negligible, one important aspect of these findings is that all changes in enzyme expression, which have been reported in the literature and are analysed here, resulted in lower values for peak height, peak time and FWHM. This is in line with the observation of a decrease in the BOLD response, reported for example in the cortex of a rat model of Alzheimer’s disease.[38] More importantly, is has been suggested that the increased enzyme levels observed in the glycolytic and the pentose-phosphate pathway can be taken as a proxy for changes in the ratio of neurons and glia cells.[4] If this is the case, their impact on the BOLD shape may be much higher than our initial analysis suggests, as discussed below.

Astrocytosis influences the BOLD shape. The scenario with the most pronounced influence on the BOLD shape models the impact of astrocytic hypertrophy or hyperplasia, an important feature of AD neuropathology. The analysis presented above shows that such changes in the voxel composition exercise much more control on the resulting BOLD shape than changes in enzyme concentration alone. This result indicates that the close interaction between neurons and astrocytes in the regulation of brain energy metabolism might be a decisive factor for the shape of the BOLD response. Given that the observed changes were reported for the hippocampus and the inferior temporal cortex, two regions affected early and extensively in AD, this might be the most relevant result of this study.

The results presented above have been obtained by in-silico experiments, which allow to fix most of the observables and thereby focus on the contribution of individual changing parameters. Our model serves the particular purpose to quantify the contribution of individual metabolic reactions in an otherwise unchanged environment. Our analysis therefore complements the clinical analysis of the BOLD response, in which many other aspects, such as stimulation-induced changes in ionic currents and vascular properties might be affected as well. Mathematical modelling, on the other hand, is a suitable tool to delineate the individual contributions and thereby increase the understanding of complex dynamic phenomena.

While our model represents an important step towards a comprehensive model of brain energy metabolism and the vascular response, two important aspects, which have been also absent in previous models, are not included. The neurovascular coupling in our model is currently realized via time-dependent changes in both metabolism and blood flow. While this is current practice, we expect that models that can include this connection might be able to shed more light on the connections between metabolic activity and the BOLD response. The second aspect, which has the potential to increase the accuracy of the results, is the mechanistic coupling of mitochondrial activity and the production of reactive oxygen species (ROS) with the demand of NADPH. While there exist several publications modelling the production of ROS in the mitochondria, the current state of the art does not allow an inclusion of these models into the model presented here. This is due to the complex nature of the electron transport chain, the main source of ROS, and the lack of time-series data of this process. The available models restrict their analysis to the steady state behaviour[41,42] or use a rule-driven modelling approach to circumvent the difficulties associated with the complexity of kinetic models.[43] We restricted ourselves therefore to mimicking the increase in NADPH demand by variations in the NADPH oxidase reaction. Finally, a feature of the model that might need further attention in future studies is the glutamate-glutamine cycle. In the model presented here, we used the approach chosen by Cloutier et al.,[12] i.e. a simplified description of glutamate shuttling coupled to sodium exchange. Using this approach, it is possible to couple energy requirements in neurons and astrocytes and thereby analyse the close interaction of the different cell types. However, there is increasing evidence that the parameters chosen by Cloutier et al.[12] need to be reconsidered. In particular, DiNuzzo et al.[44] very recently reported a higher ratio of Na + /glutamate and argued that potassium instead of glutamate might be the decisive factor of metabolic coupling of neurons and astrocytes.

To use our model in the clinical analysis of BOLD fMRI data, the differences in peak height, peak timing and FWHM need to be expressed in terms of parameter adjustments of the HRF supplied by the image analysis software. Our results indicate that not only the peak height, and to a smaller extent also peak timing, of the hemodynamic response are affected in Alzheimer’s disease. Instead, it is its width that is most sensitive to changes in the metabolic network. Comprehensive mathematical models of brain energy metabolism such as the one presented here can play an important role in the creation of personalised and disease-specific HRFs. The current state of the art allows accounting for the individual variability of the concentration of blood metabolites such as glucose and oxygen. Future models, however, may also offer the possibility to include information about differences in vascular determinants of the HRF such as arterial compliance. While these results need clinical validation, they offer a potential path towards future disease-specific HRFs.