https://orcid.org/0000-0002-1906-7198
Figures
Abstract
The neurovascular and neurometabolic couplings (NVC and NMC) connect cerebral activity, blood flow, and metabolism. This interconnection is used in for instance functional imaging, which analyses the blood-oxygen-dependent (BOLD) signal. The mechanisms underlying the NVC are complex, which warrants a model-based analysis of data. We have previously developed a mechanistically detailed model for the NVC, and others have proposed detailed models for cerebral metabolism. However, existing metabolic models are still not fully utilizing available magnetic resonance spectroscopy (MRS) data and are not connected to detailed models for NVC. Therefore, we herein present a new model that integrates mechanistic modelling of both MRS and BOLD data. The metabolic model covers central metabolism, using a minimal set of interactions, and can describe time-series data for glucose, lactate, aspartate, and glutamate, measured after visual stimuli. Statistical tests confirm that the model can describe both estimation data and predict independent validation data, not used for model training. The interconnected NVC model can simultaneously describe BOLD data and can be used to predict expected metabolic responses in experiments where metabolism has not been measured. This model is a step towards a useful and mechanistically detailed model for cerebral blood flow and metabolism, with potential applications in both basic research and clinical applications.
Author summary
The neurovascular and neurometabolic couplings are highly central for several clinical imaging techniques since these frequently use blood oxygenation (the BOLD signal) as a proxy for neuronal activity. This relationship is described by the highly complex neurovascular and neurometabolic couplings, which describe the balancing between increased metabolic demand and blood flow, and which involve several cell types and regulatory systems, which all change dynamically over time. While there are previous works that describe the neurovascular coupling in detail, neither we nor others have developed connections to corresponding mechanistic models for the third aspect, the metabolic aspect. Furthermore, magnetic resonance spectroscopy (MRS) data for such modelling readily is available. In this paper we present a minimal mechanistic model that can describe the metabolic response to visual stimuli. The model is trained to describe experimental data for the relative change in metabolic concentrations of several metabolites in the visual cortex during stimulation. The model is also validated against independent validation data, that was not used for model training. Finally, we also connect this metabolic model to a detailed mechanistic model of the neurovascular coupling. Showing that the model can describe both the metabolic response and a neurovascular response simultaneously.
Citation: Sundqvist N, Sten S, Thompson P, Andersson BJ, Engström M, Cedersund G (2022) Mechanistic model for human brain metabolism and its connection to the neurovascular coupling. PLoS Comput Biol 18(12): e1010798. https://doi.org/10.1371/journal.pcbi.1010798
Editor: Anders Wallqvist, US Army Medical Research and Materiel Command: US Army Medical Research and Development Command, UNITED STATES
Received: May 25, 2022; Accepted: December 7, 2022; Published: December 22, 2022
Copyright: © 2022 Sundqvist et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Data Availability: All supplementary material required for reproducing the results is provided in the supplementary material and is also available at: https://gitlab.liu.se/nicsu70/nvc_metabolismmodel.
Funding: This work was supported by the Swedish Research Council (2018-05418 and 2018-03319, GC; 2018-03391, ME). Additional support came from CENIIT (15.09, GC), from the Swedish foundation for strategic research (ITM17-0245, GC), from SciLifeLab and KAW (2020.0182, GC), from the H2020 project PRECISE4Q (777107, GC), from the Swedish Fund for Research without Animal Experiments (F2019-0010, GC), from the Swedish Brain Foundation (ME), from ELLIIT (GC), from VisualSweden (GC), and from VINNOVA (2020-04711, GC). The funders have had no role in the study design, data collection and analysis, decision to publish, or preparing the manuscript.
Competing interests: The authors have declared that no competing interests exist.
1 Introduction
The brain is the central organ of the nervous system in humans and it has many complex functions. Despite that it only constitutes approximately 2% of the body weight of an adult individual, the brain is responsible for 20–25% of the body’s overall energy consumption [1,2]. The brain’s main source of energy comes from oxidation of glucose. Therefore, the brain requires a continuous supply of glucose and oxygen and while the brain can metabolise stored glycogen to produce glucose, part of the required glucose and all of the oxygen is delivered via the blood [3–5]. For these reasons, the cerebral metabolic glucose consumption is directly linked to the regional cerebral blood flow [6,7]. The regional cerebral blood flow is, in turn, tightly coupled to the neuronal activity [8]. An increase in neuronal activity will lead to an increased regional cerebral blood flow, which provides an increased supply of glucose. These couplings between the neuronal, hemodynamic, and metabolic activity are commonly known as the neurovascular coupling (NVC) and the neurometabolic coupling (NMC).
More specifically, the NVC describes how neuronal activity affects cerebral blood volume (CBV) and cerebral blood flow (CBF), while the NMC describes how neuronal activity affects the cerebral metabolic activity such as the cerebral metabolic rate of oxygen (CMRO2) [9–11]. The NVC is a cornerstone of functional magnetic imaging (fMRI), which is a key method for understanding how different stimuli affect regional neuronal activity [12]. Data from fMRI is based on analysis of a blood-oxygenation-level-dependent (BOLD) signal [12,13]. The BOLD signal is governed by the regional balance between the hemodynamic and metabolic activity i.e., by the balance between oxygen supply and oxygen consumption. Both the supply of and the demand for oxygen are regulated by the neuronal activity. This connection between oxygen regulation and neuronal activity is the basis for why the BOLD signal is used as a proxy for the neuronal activity [14,15]. However, using the BOLD signal to draw correct conclusions is not straightforward since the crosstalk between hemodynamic, metabolic, and neuronal activity is nonlinear, time-varying, and partially unknown. One approach to deal with this complexity is mathematical modelling.
Mathematical modelling of the NVC has evolved over time. Initial mathematical modelling of the NVC relied on statistical interpretations of the canonical hemodynamic response function [16–18] (Fig 1A). These models provide estimates of the location and timing of neuronal activity, given a measured BOLD response, but fail to provide any mechanistic insight into the NVC, and the interplay between BOLD, CBV, CMRO2, etc. To remedy this, later mathematical modelling included some mechanistic insight and resulted in the Balloon model, developed by Buxton et al. 1998 [19–21]. The Balloon model describes the interplay between neuronal activity, CBF, CMRO2, and the BOLD response, but does so using phenomenological expressions, and thus lack mechanistic interplays between specific substances and intermediaries (Fig 1A). To take NVC modelling one step further, we have therefore created a new type of systems biology modelling approach, which incorporates such mechanistic interplays into the models (Fig 1B, top). This modelling approach first demonstrated that the BOLD response is not controlled by a negative feedback but by a positive feed-forward mechanism [22]. This initial paper was then extended to also explain the negative BOLD response [23]. Both these papers are based exclusively on human data. Lately, our models have also added cell-specific contributions, by also incorporating both rodent optogenetic and primate data [24,25]. All these more mechanistic NVC models still have a quite simple description of the metabolism involved during neuronal activity. However, regarding metabolism modelling without NVC, there are some more detailed models that describe for instance the compartmentalization of the metabolism between different neurons [26–29]. Furthermore, there are works that integrate the metabolic response with NVC, but these works rely on the phenomenological Balloon model [30–33]. It is fairly common that the metabolic parts of these models utilize magnetic resonance spectroscopy (MRS) data [34–36] (Fig 1B, bottom). However, there are still important MRS data that have not been incorporated in any metabolic model (Fig 1B, middle). Also, no such metabolic model has been integrated into a mechanistically detailed NVC model (Fig 1A).
A. A table summary of different models and what aspects of the neurovascular coupling they cover [16–19,21,23–33,55–57]. B. A schematic overview of how this work connects pre-existing models for the NVC with a description for the cerebral metabolism and how this new interconnected model can be used for informative simulations. C: A detailed illustration of the metabolism model precented in this work. Neuronal activity triggers increased consumption of glucose, which triggers downstream signaling cascades of different metabolites, which can be captured using MRS. D. A schematic illustration of the modelling cycle used to develop a minimal model.
Herein, we utilize both modelled and non-modelled MRS data to construct a new minimal model for NMC (Fig 1C). We use a minimal modelling approach to gain mechanistic insight into the NMC without introducing unneeded complexity. Statistical tests confirm that this model can describe both estimation data (Figs 2 and 3) and independent validation data (Figs 4 and 5), not used for model training. The model is also connected to our previously developed detailed NVC model [23], and the combined model can still describe the previous BOLD-data. The combined model can be used to obtain simulated predictions for the metabolic activity for complex paradigms, involving multiple stimulations, where MRS is currently not available (Fig 6).
The data gathered from Lin et al. 2012 [34] shows the change in metabolic concentrations as a response to two different visual stimulus paradigms. The first paradigm (A-D) consisted of a single visual stimulus lasting for 13.6 minutes, indicated by the black bar at the bottom of the plots. The second paradigm (E-H) consists of two 9.9-minute periods of stimulus, separated by a 9.9-minute rest period for each. For both stimulation paradigms four metabolites were analysed: lactate (A and E), glutamate (B and F), glucose (C and G), and aspartate (D and H). The data is here illustrated as a relative change from a baseline, with a mean value ± SEM (red error bars). The model estimation is illustrated as the blue shaded areas.
The model estimation of the BOLD response. A. The model estimation (blue shaded area) of the positive BOLD response fitted to data (red error bars) to a 0.5 second stimulus. B. The model estimation (blue shaded area) of the negative BOLD response fitted to data (red error bars) to a 0.5 second stimulus. The data is shown as mean values ± SEM and is taken from Sten et al. [23].
The model predictions (blue shaded areas) of the changes in metabolic concentrations for lactate (A) and glutamate (B) as a response to a double stimulation paradigm. The stimulation consisted of two 5-minute periods of stimulation (black bars) interwoven by a 5-minute rest period. The data is illustrated as mean values ± SEM (red error bars) and was gathered from Schaller et al. 2013 [36].
The model predictions (blue shaded areas) of the changes in metabolic concentrations for lactate (A), Glutamate (B), Glucose (C), and aspartate (D) as a response to a double stimulation paradigm. The stimulation consisted of two 5.3-minute periods of stimulation (black bars) interwoven by a 5.3-minute rest period. The data is illustrated as mean values ± SEM (red error bars) and was gathered from Bednařík et al. 2015 [35].
A. An illustration of how the maximum amplitude of the metabolic response (grey shaded areas) changes with the length of stimulation. B. Shows a model prediction (blue shaded area) of a BOLD response for two short 0.5 seconds stimuli 4 seconds apart and corresponding data (red error bars). The data is gathered from Lundengård et al. [22]. C-F. shows the model predicted metabolic response to the same short double stimuli. C. shows the metabolic response for lactate; D. for glutamate; E. for glucose; and F. for aspartate.
2 Methods
2.1 Model formulation
The model presented herein is formulated using ordinary differential equations (ODEs) and can in general be denoted as
(1A)
(1B)
(1C)
where x is a vector describing the model states e.g., the concentrations of the different metabolites, signalling intermediates, etc; where
describe the derivative of these states with respect to time, denoted t; where θ is a vector of the constant unknown parameters; where u is the model input, here representing a visual stimulus corresponding to the experimental setup; where x0 is the state values corresponding to an initial time point t = 0; where
are the observed model properties, here the percental difference in metabolite concentrations as well as the BOLD response; and where f and g are smooth nonlinear functions. This general description of an ODE, as well as the general modelling methods in Sections 2.2, 2.4, and 2.6, are described in a similar way as in many previous papers e.g., [23].
2.2 Development of a minimal model: general workflow
The model presented herein was developed according to an iterative modelling cycle based on hypothesis testing and rejection (Fig 1D). To start with, a falsifiable hypothesis regarding the relevant mechanisms is formulated based on literature and prior knowledge. This hypothesis is then combined with existing experimental data, for the system of interest, to formulate a mathematical model. These mechanistic models are developed as minimal models, which means that they are as simple as possible while still explaining relevant experimental data. Once formulated, each mathematical model is fitted to the experimental data and evaluated with respect to how well the model can describe the data. Such evaluation is done via qualitative inspection of the model simulation but also through quantitative statistical tests; here a χ2-test. If a model structure cannot satisfactorily describe the data, the model and the corresponding hypothesis is rejected and needs to be revised and then fitted again. If the model structure is not rejected, further model analysis is performed, to determine model identifiability and to generate well-determined predictions, which we call core-predictions [37]. Such core-predictions are ideally generated such that testing them experimentally will generate further biological insights into the performance of the model. Thus, model predictions can act as the basis for further experiment design that in turn generates or find more relevant experimental data. Finally, if this new data agrees with the model predictions, it serves as a support, or partial validation, of the model structure, since such agreement indicates that the model can explain more properties and states than just the experimental data used for model training. A thus validated model can then be used to generate additional predictions, within similar applications as the original development and testing of the model. In this project, we have developed the model in this iterative fashion, and below we present the final accepted model.
2.3 Model structure
2.3.1 Model structure for the metabolism model.
All events in the model starts with a stimulation. This stimulation represents a visual stimulus and is implemented in the form of a simple step function [23,24,38,39] formulated as:
(2)
where t is the time; ton is the time stimulation is turned on and toff is the time the stimulation ends.
To turn this step function into a softer response function, we introduce a second stimulus variable that increases depending on the step function u(t), and has a built-in decrease proportional to the level of activation, i.e.:
(3)
where Stimulus is a state that represent the effect of the visual stimuli; kstim1 describes the rate of stimulus activation and kstim2 describes the rate of stimulus decay.
The electrical activity in the neurons triggers a metabolic response involving many different metabolic reactions. Herein, we consider a simplified network of metabolites and reactions, for which we have available measurements. More specifically, the simplified metabolic network described by this model consists of around a dozen reactions. These reactions represent the metabolic pathways in and directly around the tricarboxylic acid cycle (TCA cycle). At the top of this network, glucose diffuses across the blood brain barrier from blood vessels into the neuronal tissue (Gluct). The rate of this reaction is described by the variable VGluc, which is assumed to be irreversible and governed by mass action kinetics. Following the inflow, glucose is consumed and we here assume that this consumption occurs only via glycolysis. The rate of glycolysis (Vglycolysis) is assumed to have a basal rate and a Stimulus activated rate, both governed by mass action kinetics. With these assumptions, the ODE for Gluct is given by:
(4)
where Gluct describes the concentration of glucose in the neuronal tissue and where the kinetic rate parameters for inflow and consumption of Gluct are described by
, respectively.
In the model, we assume that the plasma glucose concentration close to the activation site can be modified i.e., that Glucc is a time-dependent state. In contrast, we assume a constant inflow of glucose to the local circulation since the visual stimulus is assumed to only impact the local metabolism.
(5)
where Glucc describes the concentration of glucose in the local circulation, and where Glucoseblood represents a constant inflow of glucose from the wider circulation.
Glucose in the tissue is metabolised in intracellular glycolysis which consists of several reactions. In this model, glycolysis is simplified into a single reaction that converts glucose to pyruvate. The rate expression for this is described by Vglycolysis, which was included in Eq (4). Pyruvate in turn is then either converted into lactate or sent to the TCA cycle. The rate expression for conversion to lactate is described by VLDH, which for simplicity is given by mass action kinetics. The steps leading to the TCA cycle can occur either via pyruvate carboxylase (PC) or via pyruvate dehydrogenase (PDH). The rate expression for PC is given by VPC and assumes Michaels-Menten kinetics. The second reaction going into the TCA cycle, PDH, also assumes Michaels-Menten kinetics but is given by a more complex rate expression, dependent on two substrates. The reason for this complexity is that the pyruvate that enters the TCA cycle via the PDH pathway is combined with oxaloacetic acid (OAA), which then, through several reactions, is turned into 2-oxoglutarate (OG). These multiple rection steps are in the model simplified into a single reaction (with reaction rate VTCA1), which describes the first part of the TCA cycle (Fig 1C). The rate expression for VTCA1 is thus dependent on two substrates and is assumed to have a saturation with respect to both Pyruvate and OAA. With all these reactions included the ODE for pyruvate becomes:
(6)
where Pyr describe the concentrations of pyruvate;
,
, and
are kinetic rate parameters; and where
,
and
are Michaelis-Menten constants. Note that
,
describe the maximum rate of VTCA1 and VPC, respectively, while
describes the rate of conversion from pyruvate to lactate.
Lactate is not considered in detail in this model and is simply described by a conversion from pyruvate (VLDH) and a clearance term (Vclear1). The differential equation for lactate is given by:
(7)
where Lac is the concentration of lactate and k1 is the rate parameter describing the degradation rate of lactate.
As described earlier, the model condenses the multiple reactions of the TCA-cycle that converts oxaloacetic acid (OAA) and pyruvate to 2-oxoglutarate (OG) into a single reaction TCA1 (reaction rate: VTCA1). Further, the remaining reactions that returns OG to OAA is also simplified into a single reaction TCA2 (reaction rate: VTCA2). Thus, the reactions TCA1 and TCA2 create a cyclical relationship between the model states of OAA and OG, which together make out the TCA cycle. Additionally, OAA can undergo transamination to form aspartate, this is described in the model as VGOT. With these reactions, the full differential expression for OAA is thus described as:
(8)
where OAA is the concentration of oxaloacetic acid in the tissue and where
, and
are the kinetic rate parameters describing the reaction rates of VTCA2 and VGOT, respectively.
Analogously, OG is synthesised via TCA1 and is consumed by TCA2, but OG can also be transaminated to form glutamate. This reaction from OG to glutamate is denoted XM (reaction rate: Vxm). The full ODE for OG is thus given by:
(9)
where OG is the concentration of 2-oxoglutarate, and
is the kinetic rate parameter for the transamination of OG.
Lastly, the model describes the concentrations of three metabolites on the periphery of the TCA cycle that are vital to cerebral metabolism. The first of these metabolites is aspartate, which is converted from OAA via VGOT. The subsequent metabolism of aspartate is implemented as a simple reaction denoted Vclear2
(10)
where Asp is the concentration of aspartate and
is the kinetic rate parameter describing the rate of further metabolization of aspartate, not described in detail by this model.
The second of the peripheral metabolites is glutamate, which as described is transaminased from OG. Glutamate is then converted into glutamine, the third peripheral metabolite, via glutamine synthetase (rate expression: VGS). Reversely, glutamine is converted back to glutamate via phosphate-activated glutaminase (rate expression: VPAG). Finally, glutamate is decreased by Vclear3, which is a simplified description for further metabolism and utilization of glutamate in the neuronal tissue. All these reactions and simplifications means that the ODE for glutamate is implemented as:
(11)
where Glut is the concentrations of glutamate;
and
are the kinetic rate parameters for the forward and reverse conversion of glutamine to glutamate, respectively; and
is the kinetic rate parameter of any further metabolism of glutamate, not described in detail in this work.
The last metabolite described in this model is glutamine, and the only glutamine dynamics described in this model is governed by the already described VGS and VGLUT rate expressions, i.e.
(12)
where Gln is the concentration of glutamine.
2.3.2 Model structure for the NVC model.
The model structure for the NVC model was taken from the model presented by Sten et al. 2017 [23]. This model gives a detailed physiological description of how neuronal activity affects the CBF, CBV, and the cerebral metabolic activity to generate a BOLD-response. For a detailed description and motivation of this model, please see the supplementary S1 Table or the original work [23], and all reactions in the complete model, along with simulation files are specified in the S1 Supporting Information file.
2.3.3 Connection between NVC model and metabolism model.
The two models–the new metabolism model in 2.2.1 and the Sten et al. 2017 model–can be considered as two separate standalone models. However, we here also present and analyse a connected version. In this connected version, Eqs (3), (4), and (5) were removed from the metabolic model since corresponding reactions already exist in the Sten et al. 2017 model. Further, the reaction named Vglycolysis in Eq (6) was replaced with the reactions baseMet (reaction rate: VbaseMet) and stimMet (reaction rate: VstimMet), which are a part of the Sten et al. 2017 model [23]. The first of these reactions, baseMet describes the basal non-stimulated metabolism of glucose and oxygen, and stimMet describes the increase in metabolism caused by neuronal activity. With these alterations implemented, the ODE for pyruvate is described by:
(13)
where kbasalMet, kprop1 and kprop2 are kinetic rate parameters; where GlucoseA and O2A are the respective concentrations of glucose and oxygen in the neuronal tissue; where DelayM represents an intermediate state of the stimulus effect on the metabolism. The variable GlucoseA describes the same quantity as the state Gluct described in the removed Eq (4). The reason for this change of name in the combined model is that GlucoseA was an already implemented state in the NVC model. The parameters and states of the NVC model are described in further detail in Sten et al. 2017 [23], and both that and the combined models are described in the supplementary S1 Table.
In this combined model, as in the two constituent sub-models, the observed model properties, or model observables, consist of four metabolic states: Lac, Glut, Gluct, and Asp. Additionally, a fifth model observable is the BOLD-signal defined in Sten et al. 2017 as a signal proportional to the amount of deoxy haemoglobin (dHb). The BOLD-signal is formally defined as:
(14)
where dHb is the amount of deoxy haemoglobin and ky is a scaling parameter.
More specifically, each of the five model observables are defined as a percentage difference from a steady-state baseline. In practice, this is achieved by dividing the observable value with its corresponding steady-state value, withdrawing one and multiplying with 100%:
(15)
where Xi is one of the observable properties Lac, Glut, Gluct, Asp, or BOLD; where ssXi is the corresponding steady-state value, and where kyi is the ith scaling parameter.
2.4. Model evaluation
2.4.1. Parameter estimation.
The model parameters (θ) are estimated from experimental data (D), by minimizing the negative log-likelihood function J:
(16)
where θ is the model parameters; ne is the number of experiments e,
is the number of observations in experiment e;
is the number of samples at observation o in experiment e;
is the measured data sample point;
is the model simulation for the corresponding sample point; and
is the standard deviation of the sample point.
The function J(θ) is implemented as the objective function for a suitable optimization algorithm (see section 2.6) and optimized with respect to the parameters θ.
2.4.2 Identifiability analysis.
For the model presented herein two types of identifiability analysis were performed. A structural identifiability analysis to evaluate if the model structure resulted in any unidentifiable model states or parameters, and a practical identifiability analysis to evaluate what simulated model states and parameters can be accurately estimated from the experimental data. For the structural identifiability analysis, a local identifiability analysis was performed on the metabolism model described in section 2.3.1. An implementation of the recently developed methods described by Thompson et al., 2022 [40] was used. These methods extend the algorithm presented by Sedoglavic et al. in 2002 [41] to determine local identifiability of states and parameters when different numbers of derivatives of signals are available. This analysis showed that all model states and parameters were identifiable when at least eight derivatives of the measurement signals was available. A complete report of these results is presented in the supplementary S2 Table.
For the practical identifiability analysis, the uncertainties of the simulated model states and the parameters were estimated through a Markov chain Monte Carlo (MCMC) sampling approach. A posterior distribution of the parameter values was generated, using 105 samples, and all parameter sets found acceptable according to a χ2-test was collected. The uncertainty is then determined by the confidence interval:
(17)
Where α is the confidence level; Δα(χ2) is the α quantile of the χ2 statistic [42,43]; DoF is the degrees of freedom; and θ* are the optimal parameters. In this work, the DoF is equal to the number of model parameters (n = 58). The optimal parameter values and the boundaries of the parameter distributions can be found in the supplementary S3 Table.
2.5 Experimental data
In this work, experimental data from previously published studies have been used for model evaluation. Mainly, H-MRS data from human subjects exposed to visual stimuli have been considered for the model evaluation presented herein. Here follows a condensed description of the data, for further details please see the original publications for each dataset. When necessary, data were extracted from publication figures, using the WebPlotDigitizer tool [44].
2.5.1 Metabolic estimation data.
In their work, Lin et al. 2012 [34] present metabolite data gathered from the visual cortex of human subjects using functional MRS. Data was gathered form ten healthy subjects (7 males, 3 females, 25±3 years old) and the data used herein consists of two sets of four metabolites. The first data set describes the relative change in metabolic concentrations in the primary visual cortex for lactate, glutamate, glucose, and aspartate, for a single visual stimulus. The second data set describes the relative change in metabolic concentrations for the same four metabolites but for a double visual stimuli paradigm. The visual stimulus consisted of contrast-defined wedges, moving towards and from a central point. The single stimulus paradigm consisted of a 6.6-minute baseline followed by a 13.2-minute exposure to the visual stimulus, followed in turn by a 19.8-minute recovery. In the double stimulus paradigm, the baseline was established for 9.9 minutes followed by a two 9.9-minute periods of exposure to the visual stimulus, with a 9.9-minute rest period in between.
2.5.2 BOLD estimation data.
The data used to train the BOLD-response is gathered from Witt et al. 2016 [45] and consists of fMRI data from a visual-motor task. The data was gathered from 11 healthy subjects (5 males, 6 females, 21–28 years of age) that was tasked to push a button in response to a certain visual stimulus. The visual stimulus was shown for 0.5 seconds, and a BOLD-response was extracted from each subject. The positive BOLD responses were gathered from the bilateral primary cortex, while negative BOLD-responses were gathered from the posterior cingulate cortex. The BOLD data is normalized as a percentage signal change.
2.5.3 Metabolic validation data.
The model presented herein was validated with respect to independent validation data i.e., data that was not been used for model training. This data was taken from two independent studies. The first such study is the work presented by Schaller et al. 2013 [36]. In this work, ten healthy subjects (9 males, 1 female, 20–28 years of age) were exposed to a visual stimuli and changes in metabolic concentrations for lactate, glutamate was gathered via MRS. The data was sampled at 75 time points over a 25-minute period. This period consists of a 5-minute baseline period followed by two periods of exposure to visual stimuli for 5-minutes each. These exposure periods were interwoven by a 5-minute rest period and followed by a 5-minute recovery period.
The second study from which validation data were gathered is the work presented by Bednařík et al. 2015 [35]. In this work, the authors present MRS data that shows how metabolic concentrations of glutamate, aspartate, lactate, and glucose change as a response to a visual stimulus. The data is gathered form fifteen healthy subjects (7 males, 8 females, 33±13 years of age). The visual stimulus consisted of a red and black flickering checkerboard. For the metabolite measurements the stimulation sequence consisted of a 5.3-minute baseline acquisition, followed by a stimulation-rest-stimulation-rest sequence were each block lasted 5.3 minutes. The metabolic concentration data was sampled with a 2.7-minute resolution, for a total of nine sample points for each metabolite. For the model validation the model was used to simulate the metabolic response that corresponded to the experimental setting of the validation data, described above. The parameters used for this simulation were the optimal parameters, with respect to the estimation data, obtained from the parameter estimation described in Section 2.4.1.
2.6 Model implementation
Model implementations and simulations were done in the “Advanced Multi-language Interface to CVODES and IDAS” (AMICI) toolbox [46–48]. The implementation of the equations described in section 2.3 and the subsequent model analysis were done in MATLAB by The MathWorks, Inc. releases 2017b and 2020a. Model parameter estimation were done using the “Metaheuristics for systems biology and bioinformatics global optimization” (MEIGO) toolbox [49] with the enhanced scatter search (ess) algorithm.
3 Results
We have developed a new model for the central metabolism in the visual cortex and connected this model to our existing model for the NVC [23]. The model is described in detail in Materials and Methods and all simulation files are available in the supplementary information (S1 Supporting Information). The model has been developed using experimental data from two different studies, Lin et al. 2012 [34] (Fig 2) and Sten et al. 2017[23] (Fig 3). This data was used to estimate the model parameters. Additionally, the model has been evaluated with experimental data gathered from Schaller et al. 2013 [36] (Fig 4) and Bednařík et al. 2015 [35] (Fig 5). This data was used as independent validation data to assess the model’s ability to explain data for which it had no prior knowledge. This provides some support for the included mechanisms and predictive capability of the model. Finally, we demonstrate the usefulness of the model by showing how one can provide reasonable simulation results of metabolism in experiments where metabolism has not been measured. In these simulations, we also used the model to predict how long a stimulation would need to be to generate a measurable metabolic response (Fig 6).
3.1 Metabolic estimation results
As mentioned, part of the dataset used to train the model was gathered from Lin et al. 2012 [34], described in further detail in Section 2.5.1. Briefly, this data describes the relative change in metabolic concentrations in response to either a single stimulus (Fig 2A–2D, error bars) or a double stimulus (Fig 2E–2H, error bars). More specifically, the single stimulation consisted of a sustained visual stimulus for 13.6 minutes (Fig 2A–2D, black bar). In contrast, the double stimuli paradigm consists of two periods of sustained visual stimuli for 9.9 minutes, separated by a 9.9-minute rest period (Fig 2E–2D, black bars). The experimental data is presented as mean ± SEM (Fig 2, red error bars).
The model was simultaneously fitted to both the metabolic data in response to the single and double stimulation paradigm in Lin et al. 2012 [34], as well as to the BOLD data in Sten et al. [23] (presented in Section 3.2). The resulting agreement between data (error bars) and the simulation (shaded areas) is shown in Figs 2 and 3. As can be seen, for the single stimulus, the model clearly describes the increase in lactate (Fig 2A) and glutamate (Fig 2B) and the decrease in glucose (Fig 2C) and aspartate (Fig 2D), in response to the visual stimuli. Similarly, the double stimuli paradigm shows a model behaviour that is consistent with the experimental data (Fig 2E–2H): during the visual stimulation the lactate and glucose levels increase (Fig 2E and 2F) and the glucose and aspartate levels decrease (Fig 2G and 2H), relative to their baseline. These visual assessments are also supported by a statistical χ2-test (J(θ*) = 213.44, threshold: χ2(α = 0.05, DoF = 191) = 224.24) (Materials and Methods).
3.2 The model simultaneously agrees with BOLD data
The data used to train the model’s BOLD-response originates from Witt et al. [45] and was used to train the model in Sten et al. [23]. The data consists of a positive BOLD-response (Fig 3A) and a negative BOLD-response (Fig 3B). This data is described in more detail in Section 2.5.2. The stimulus used is a short (0.5 s) visual stimuli at t = 0.1 (Fig 3, black line), and the data is presented as a percental mean ± SEM deviation from the baseline (Fig 3, red error bars). As described in the previous section, the model was simultaneously fitted to both these BOLD data and the MRS data (Fig 2). The agreement between the BOLD data and simulation (blue shaded areas) is presented in Fig 3. As can be seen, the model displays all the expected features of an archetypical BOLD response: a small initial dip around 2 seconds after the stimulation, a larger main response at around 5–7 seconds and a post-peak undershoot between around 10–16 seconds [22,23]. This response is in good agreement with the data, which again is formally supported by a χ2-test (J(θ*) = 213.44, threshold: χ2(α = 0.05, DoF = 191) = 224.24).
3.3 Model validation against additional MRS data
Once the model had been trained on the MRS data gathered from Lin et al. 2012 [34] and BOLD data from Sten et al. [23], the model was validated against the data gathered from Schaller et al. [36] (Fig 4) and Bednařík et al. [35] (Fig 5). The MRS data gathered from Schaller et al. describes the change in metabolic concentrations of lactate and glutamate in response to two visual stimuli, with a short break in between. The data is herein presented as a relative change in mean concentration with the SEM displayed as error bars (Fig 4, red error bars). The model predicted dynamics for the concentrations of lactate (Fig 4A) and glutamate (Fig 4B) which are displayed as the blue shaded area, and which is in very good agreement with the data. For changes in lactate concentration, the model provides a good prediction of how lactate increases as a response to the visual stimulus i.e., the same increase can be seen in the data (Fig 4A). For the changes in glutamate, the model predicts a slight increase when the visual stimulus is turned on. This is in accordance with the data. However, when the visual stimulus is turned off (at 600 s and 1200 s) the model describes a decrease in glutamate levels, but this decrease is slower than the corresponding decrease seen in the data. This makes the model prediction slightly too high for all points following the first stimulus period (Fig 4B). All in all, given that the model is not trained for any of these data, the agreement between model and data is good, both qualitatively and quantitative, and comparable to the size of the experimental variations and uncertainties.
Further, the model was simultaneously validated using the data gathered from Bednařík et al. [35]. This data also describes changes in metabolic concentrations as a response to a sequence of visual stimuli, further described in Section 2.5.3. This data set contains data for lactate (Fig 5A), glutamate (Fig 5B), glucose (Fig 5C), and aspartate (Fig 5D), and again the data is presented as the mean concentrations ± SEM indicated by the red error bars. Once again, the model predicts the changes in metabolic concentrations with good accuracy (Fig 5, blue shaded areas). Similarly, to previous results, the model predictions for lactate and glutamate show an increase when the visual stimulus is turned on and as previously this is supported by the data (Fig 5A and 5B). Also, similarly to before, the predicted recovery rate for glutamate appears slightly too slow when compared to data (Fig 5B). Further, for glucose and aspartate, the model-predicted levels are slightly lower than what is indicated by the data. For glucose, the model prediction is consistently below the data points after the first period of stimulus (Fig 5C). For aspartate, the model prediction is consistently below the data points following the start of the first stimulus period (Fig 5D). All in all, the overall agreement between model predictions and validation data is qualitatively correct, and quantitatively the simulation-data differences are comparable to the experimental uncertainty, at least for lactate, glutamate, and glucose.
3.4 Using the final model: simulating non-measured variables
Once the model had been validated, we used the model to predict various variables and conditions for which we have no data. In particular, we investigated how long a visual stimulation needs to be to generate a measurable response in metabolic concentrations. This was achieved by evaluating the amplitude of the metabolic responses for glucose, lactate, aspartate, and glutamate, for different stimulation lengths (Fig 6A, grey shaded areas). The length of the stimulations was varied from 0.5 s to 10 000 s. As to be expected, the amplitude of the metabolic response increases as the length of the stimulation increases. Furthermore, our simulations suggests that the metabolic response reaches a maximal amplitude, after a few hundred seconds. For lactate and aspartate, a stimulation of around 200 s or longer yields this maximal response, for glutamate, the maximal response is reached for stimulations around 300 s, and for glucose the stimulation length needs to be around 600 s to reach a maximal response (Fig 6A).
Further, the model can be used to make predictions for non-measured variables in an existing experiment. In Lundengård et al. [22], there are BOLD response data for a series of two short 0.5 seconds stimuli, with 4 seconds in-between, but no metabolic data were collected. Our model can be used to assess what the metabolism probably looked like. The model can accurately predict this BOLD response (Fig 6B), but it also predicts the corresponding metabolic response for the same stimulation paradigm. These predictions shows that the metabolic response is slower than the BOLD response: the BOLD signal is completely over after 15–20 seconds, even though the metabolism still is maximal (Fig 6C–6F). Furthermore, for such short stimuli, the amplitude of the metabolic response is considerably lower than in the previous examples above (Figs 2, 4, and 5). This lower amplitude is consistent with the results in Fig 6A. This means that corresponding metabolic changes would not be experimentally detected, which makes sense since metabolic experiments usually have longer paradigms. The simulations, however, provides an estimation of MRS experiments with higher signal-to-noise-ratio.
4 Discussion
Mathematical modelling provides a great potential to understand the complex interplay between cerebral metabolic, hemodynamic, and neuronal activity. Herein, we present a mathematical model that can explain previously unmodelled MRS data for the central cerebral metabolism. Further, we present a first connection between a minimal metabolic model and a mechanistically detailed NVC model. We show that this combined model can accurately describe MRS and fMRI data simultaneously and that the model can also describe independent validation data, not used for model training. This MRS data, used to train the model, consisted of time series that show how concentrations for glucose, lactate, aspartate, and glutamate change as a response to visual stimuli (Fig 2). The fMRI data consisted of a positive and negative BOLD-response (Fig 3), also from a visual stimulus. Similarly, the independent MRS data, used for model validation, also consists of time-series that show the changes in metabolic concentrations as a response to visual stimuli (Figs 4 and 5). Finally, we show that this model can be used to make predictions of non-measured variables and conditions, and that it can be used to assess e.g., how long a stimulation needs to be to generate a measurable change in metabolic concentrations, and what the likely shape of the metabolism is in cases of shorter stimulations, when MRS cannot detect any signal (Fig 6).
In this study we have based our modelling of the cerebral metabolism on the MRS data gathered by Lin et al. 2012 [34], Schaller et al. [36], and Bednarík et al. [35]. However, MRS is by no means the only measuring technique employed to investigate the cerebral metabolic pathways. One more accurate, but also more demanding, approach to determine cerebral metabolic pathways is modelling combined with 13C MRS data. In this technique researchers insert isotopically labelled substrates into the metabolic system and measure the fractional enrichment of these isotopes in metabolic products [50,51]. Mathematical modelling based on this type of data then allows for accurate quantification of different metabolic pathways during different conditions e.g., before and after stimulation [52]. However, these models require specific knowledge regarding the relevant atom transitions of the isotopically labelled substance. This means that such atom-based models are more difficult to scale up, and that they mainly are used to investigate specific metabolic pathways in isolation. Furthermore, such methods are mostly applied to animal and in vitro experimental system.
In this study we also investigated how the model predicts that the length of the stimulations affects the metabolic response. Most studies that investigate the cerebral metabolism via MRS uses very long periods of visual stimulation. For instance, the study from which we have gathered data [34] uses a period of continuous visual stimuli that lasts almost ten minutes. While our model suggests that this is the time required for the maximal signal amplitude and thus the largest signal-to-noise ratio (SNR), our simulations also suggest that a detectable signal-to-noise ratio could be reached with shorter stimulations, such as 2 minutes (Fig 6A). However, it should be noted that such shorter stimulation protocols still need a longer data-collection period, since the dynamics of the metabolism is slow, and probably continues long after short stimulations are over (Fig 6). These statements are based on our interconnected model that can be used to make predictions of probable behaviours of variables in cases where data is not available. Here, this new possibility was illustrated by doing model predictions for the metabolic responses, in an experimental scenario where only BOLD responses were measured (Fig 6). This shows that two 0.5 s stimulations with 4 s in-between will give rise to two, in principles, measurable peaks within 10 s in the BOLD signal, but to slower and lower responses in the metabolites: <5% response in glucose and lactate, and <2% response in glutamate and aspartate. These expected changes are smaller than the usual experimental variability (Figs 2, 4, and 5), implying that those predictions would probably not be experimentally detectable. However, these predictions of short stimulation metabolic response are done simultaneously with a prediction of the BOLD-response (Fig 6B) which is validated by experimental data which increases the confidence that the prediction for the metabolic response is also reasonable. This likelihood is of course also supported by the fact that the model can both describe and predict other types of metabolic data. This shows that we can do a prediction of a situation that likely cannot be experimentally measured. As such, there is potential for our model to be used as a tool for experiment-design of an experiment that has shorter stimulations than the original studies, but that still would give measurable responses.
There are a couple of limitations regarding this work that should be mentioned. As is the case with all mathematical models, the model presented herein is the result of several simplifications and assumptions. For instance, several previous works have shown that aspects such as compartmentalisation between neuronal and glial tissue and glycogen metabolism have important roles in the cerebral metabolism [26–29,32,52,53]. However, adding additional complexity to the model, for instance in the form of multiple compartments, does not necessarily improve the model’s ability to explain the data we have considered in the manuscript, but it does reduce the model’s parameter identifiability. This means that the estimated confidence intervals for the parameter values might increase by an order of magnitude, or more, by adding complexity to the model. This is illustrated by an example case presented in supplementary file S1 Text, where the parameter identifiability of a two-compartment version of the model presented herein is analyzed. Consequently, the model gives less accurate estimates of the parameters that are included while not improving the explanation of the given experimental data. The reason for this reduced identifiability is that the complexity of the model will increase beyond the complexity of the experimental data. For instance, the data we used to train the model (Lin et al., 2012) [34]cannot distinguish which compartment the change in a metabolite’s concentration occurs in. Similarly, adding additional metabolic pathways such as an inflow of glycogen to the model would make the production of glucose less identifiable i.e., more uncertain. Therefore, the source of glucose could be either from the blood or from glycogen degradation. Given the considered data it would not be possible to specify the source of glucose. For these reasons we choose to maintain the minimal model approach even though this means that the model does not include certain previously modelled aspects of the cerebral metabolism.
Another such simplification is the current interpretation of the oxygen consumption as implemented in Sten et al. 2017 [23]. More specifically, this implementation means that the only part of the metabolism that affects the CMRO2 is the glycolysis step i.e., the conversion form glucose to pyruvate as described by Eq (13). This interpretation is a simplification, as oxygen is consumed both in the glycolysis and in the mitochondrial oxidative phosphorylation [54]. However, the effect of oxygen in the glycolysis seen in Eq (13) has no corresponding effect on the steps that represents the TCA-cycle (VTCA1 and VTCA2 in e.g., Eq (8)). Such an effect has not been necessary to accurately describe the data presented in this work. The inclusion of a more detailed description of the oxygen metabolism would produce a more physiologically accurate interpretation of the cerebral metabolism but also would require more data than we had access to in this work.
Finally, the metabolic model presented herein makes simplifications with respect to the units of the different states defined in the model. The metabolic concentrations have been modelled in arbitrary units and are scaled by unknown scaling parameters as defined in Eq (15). Despite these simplifications, the model is still able to accurately explain both the metabolic and hemodynamic dynamic activity that is seen in experimental data, as well as predict independent validation data.
Furthermore, the presented model uses a simplistic expression for calculating the BOLD-signal (see Eq 14). This simple expression is taken directly from the Sten et al. 2017 [23] model and does not directly take into account the effects of variables such as CBV, CBF and CMRO2. However, these variables do indirectly affect the calculated BOLD-signal as they are incorporated into the expression for dHb. There are models that gives a more detailed description of the BOLD-signal; for instance, Havlicek et.al. 2015 [55] present an expression where the BOLD signal constitutes a combination of the extravascular and intravascular signals. Further, one other more detailed version of calculating the BOLD signal, not included if in our model, is to consider the compartmentalization of the intravascular component as is done in e.g. Griffeth and Buxton 2011 [56] and Kim and Ress 2016 [57]. In these papers, the BOLD signal is calculated as a sum of contributions from arterial, capillary, and venous vessels, as well as an extravascular compartment. A more comprehensive implementation of the BOLD-signal would likely lead to a model that could more accurately describe a versatile set of BOLD-responses possibly improving the prediction of the BOLD-response seen in Fig 6B, to be more in line with the experimental data. Such an implementation is included in Sten et al. 2021 [25] where multiple sets of experimental data are used to identify key mechanisms of the BOLD-response that are preserved across different species and time-scales. The reason why such an expression for the BOLD -signal was not implemented in this work is that the model developed herein focuses on the metabolic responses and the connection to the already existing Sten et al. 2017 NVC model.
Our new combined NVC-metabolism model opens the door to several important applications in the future. First, one such application is model-based experiment-design to e.g., design new shorter stimulation paradigms, as discussed above. Second, another application is to integrate a wider variety of data into the same model-based analysis. In other words, this model could incorporate data from future experiments where e.g., both BOLD and metabolism has been measured simultaneously. Third, an extension of this application is to incorporate data for CBV and/or CBF. Since all these variables—BOLD, CBV, CBF, and metabolism—are interconnected, they should be analysed using one and the same model, to fully exploit all the information that is contained in such a joint dataset. However, for that to be possible, a more advanced interconnected model than the one presented herein is needed, since this model probably does not provide a good enough description of CBV and CBF dynamics. Nevertheless, once such an interconnected model is in place, it can be used to provide a new updated ability to use a model to infer CMRO2, compared to the highly simple model that is used presently [56]. Fourth, such integrated analyses could also be useful in clinical contexts, since multiple variables analysed together would provide an integrated understanding of the brain, with the potential for new biomarkers and mechanistic insights regarding patient screening, stratification, and monitoring [58]. Finally, a future metabolic-NVC model can also be used in multi-organ digital twin models, that also can be used for a variety of applications. All in all, our new integrated model is the first to include intracellular details regarding both metabolism and NVC (Fig 1A), and even though it still is just a first step in this new direction, it points to way towards many important applications in the future.
Supporting information
S1 Fig. An illustration of the two-compartment version of the model presented in Fig 1C.
Illustrating the division of the metabolic reactions into a glial and a neutron compartment with some cross over reactions.
https://doi.org/10.1371/journal.pcbi.1010798.s001
(TIF)
S2 Fig. Comparison of parameter distributions from the identifiability analysis.
The box chart illustrates the parameter distributions for four parameters (,
,
,
). With the distribution of the parameters in the single compartment model shown in blue and distributions of the corresponding parameters for the two-compartment model shown in red and yellow respectively.
https://doi.org/10.1371/journal.pcbi.1010798.s002
(TIF)
S1 Table. A complete list of the implemented model equations form the Sten et al.
2017 model [23]. The equations are taken with the consent of the authors of the original work and are presented here as for the convenience of the reader. Please see the original publication for a detailed description of the model equations.
https://doi.org/10.1371/journal.pcbi.1010798.s003
(DOCX)
S2 Table. The results from the structural identifiability analysis that shows how many derivatives of the measurement equations are sufficient for each model state and parameter to be identifiable.
The calculations are based on an implementation of the methods presented by Thompson et al., 2022 [40].
https://doi.org/10.1371/journal.pcbi.1010798.s004
(DOCX)
S3 Table. A list of the parameter values for the presented model.
The parameter values correspond to the values obtained for the best fit to the estimation data. The upper and lower bounds correspond to the range of acceptable parameter values obtained from the uncertainty analysis.
https://doi.org/10.1371/journal.pcbi.1010798.s005
(DOCX)
S1 Text. An example case of the problems with parameter identifiability.
This example illustrates the problems with parameter identifiability that arise when additional model complexity is introduced. In this example a two-compartment version of the metabolic model presented here is fitted to the data for the metabolic responses and the resulting parameter identifiability is analysed.
https://doi.org/10.1371/journal.pcbi.1010798.s006
(DOCX)
S1 Supporting Information. The zip-file that contains all the necessary files required to reproduce the results presented in this paper.
To ensure compatibility please ensure that all the requirements described in the provided read me file are fulfilled.
https://doi.org/10.1371/journal.pcbi.1010798.s007
(ZIP)
Acknowledgments
Preliminary analysis leading up to these results presented herein have been done by Sofie Johansson, who worked on this project in a student internship.
References
- 1. Raichle ME. The restless brain: how intrinsic activity organizes brain function. Philos Trans R Soc Lond B Biol Sci. 2015 May 19;370(1668). pmid:25823869
- 2.
Magistretti PJ, Allaman I. Brain Energy Metabolism. In: Pfaff DW, editor. Neuroscience in the 21st Century: From Basic to Clinical. New York, NY: Springer; 2013. p. 1591–620.
- 3. Brain energy metabolism by B. K. Siesjo, John Wiley & Sons, New York, 1978. Ann Neurol. 1979;5(3):308–308.
- 4. Öz G, Seaquist ER, Kumar A, Criego AB, Benedict LE, Rao JP, et al. Human brain glycogen content and metabolism: implications on its role in brain energy metabolism. Am J Physiol-Endocrinol Metab. 2007 Mar;292(3):E946–51. pmid:17132822
- 5. Rich LR, Harris W, Brown AM. The Role of Brain Glycogen in Supporting Physiological Function. Front Neurosci. 2019; 13. pmid:31749677
- 6. Sokoloff L. Local cerebral energy metabolism: its relationships to local functional activity and blood flow. Ciba Found Symp. 1978 Mar;(56):171–97. pmid:97062
- 7. Roland PE, Eriksson L, Stone-Elander S, Widen L. Does mental activity change the oxidative metabolism of the brain? J Neurosci. 1987 Aug 1;7(8):2373–89. pmid:3612246
- 8. Roy CS, Sherrington CS. On the Regulation of the Blood-supply of the Brain. J Physiol. 1890 Jan;11(1–2):85–158.17. pmid:16991945
- 9. Attwell D, Buchan AM, Charpak S, Lauritzen M, MacVicar BA, Newman EA. Glial and neuronal control of brain blood flow. Nature. 2010 Nov;468(7321):232–43. pmid:21068832
- 10. Cauli B, Hamel E. Revisiting the role of neurons in neurovascular coupling. Front Neuroenergetics. 2010; 2. pmid:20616884
- 11. Iadecola C. Neurovascular regulation in the normal brain and in Alzheimer’s disease. Nat Rev Neurosci. 2004 May;5(5):347–60. pmid:15100718
- 12. Ogawa S, Lee TM, Kay AR, Tank DW. Brain magnetic resonance imaging with contrast dependent on blood oxygenation. Proc Natl Acad Sci U S A. 1990 Dec;87(24):9868–72. pmid:2124706
- 13. Hillman EMC. Coupling mechanism and significance of the BOLD signal: a status report. Annu Rev Neurosci. 2014;37:161–81. pmid:25032494
- 14. Logothetis NK, Pauls J, Augath M, Trinath T, Oeltermann A. Neurophysiological investigation of the basis of the fMRI signal. Nature. 2001 Jul;412(6843):150–7. pmid:11449264
- 15. Shmuel A, Augath M, Oeltermann A, Logothetis NK. Negative functional MRI response correlates with decreases in neuronal activity in monkey visual area V1. Nat Neurosci. 2006 Apr;9(4):569–77. pmid:16547508
- 16. Lindquist MA, Meng Loh J, Atlas LY, Wager TD. Modeling the hemodynamic response function in fMRI: Efficiency, bias and mis-modeling. NeuroImage. 2009 Mar 1;45(1, Supplement 1):S187–98. pmid:19084070
- 17. Lindquist MA, Wager TD. Validity and power in hemodynamic response modeling: A comparison study and a new approach. Hum Brain Mapp. 2007;28(8):764–84. pmid:17094118
- 18. Loh JM, Lindquist MA, Wager TD. RESIDUAL ANALYSIS FOR DETECTING MIS-MODELING IN fMRI. Stat Sin. 2008;18(4):1421–48.
- 19. Buxton RB, Wong EC, Frank LR. Dynamics of blood flow and oxygenation changes during brain activation: The balloon model. Magn Reson Med. 1998;39(6):855–64. pmid:9621908
- 20. Buxton RB, Frank LR. A Model for the Coupling between Cerebral Blood Flow and Oxygen Metabolism during Neural Stimulation. J Cereb Blood Flow Metab. 1997 Jan 1;17(1):64–72. pmid:8978388
- 21. Friston KJ, Mechelli A, Turner R, Price CJ. Nonlinear Responses in fMRI: The Balloon Model, Volterra Kernels, and Other Hemodynamics. NeuroImage. 2000 Oct 1;12(4):466–77. pmid:10988040
- 22. Lundengård K, Cedersund G, Sten S, Leong F, Smedberg A, Elinder F, et al. Mechanistic Mathematical Modeling Tests Hypotheses of the Neurovascular Coupling in fMRI. PLoS Comput Biol. 2016;12(6):e1004971. pmid:27310017
- 23. Sten S, Lundengård K, Witt ST, Cedersund G, Elinder F, Engström M. Neural inhibition can explain negative BOLD responses: A mechanistic modelling and fMRI study. NeuroImage. 2017 Sep 1;158:219–31. pmid:28687518
- 24. Sten S, Elinder F, Cedersund G, Engström M. A quantitative analysis of cell-specific contributions and the role of anesthetics to the neurovascular coupling. NeuroImage. 2020 Jul 15;215:116827. pmid:32289456
- 25. Sten S, Podéus H, Sundqvist N, Elinder F, Engström M, Cedersund G. A multi-data based quantitative model for the neurovascular coupling in the brain. bioRxiv. 2021 Mar 26;2021.03.25.437053.
- 26. Aubert A, Costalat R. Compartmentalization of brain energy metabolism between glia and neurons: Insights from mathematical modeling. Glia. 2007 Sep;55(12):1272–9. pmid:17659526
- 27. Gruetter R, Seaquist ER, Ugurbil K. A mathematical model of compartmentalized neurotransmitter metabolism in the human brain. Am J Physiol-Endocrinol Metab. 2001 Jul;281(1):E100–12. pmid:11404227
- 28. Duarte JMN, Gruetter R. Glutamatergic and GABAergic energy metabolism measured in the rat brain by 13C NMR spectroscopy at 14.1 T. J Neurochem. 2013;126(5):579–90. pmid:23745684
- 29. Mangia S, Simpson IA, Vannucci SJ, Carruthers A. The in vivo neuron-to-astrocyte lactate shuttle in human brain. J Neurochem. 2009 May;109(Suppl 1):55–62.
- 30. Aubert A, Costalat R, Valabrègue R. Modelling of the Coupling between Brain Electrical Activity and Metabolism. Acta Biotheor. 2001;49(4):301–26. pmid:11804241
- 31. Winter F, Bludszuweit-Philipp C, Wolkenhauer O. Mathematical analysis of the influence of brain metabolism on the BOLD signal in Alzheimer’s disease. J Cereb Blood Flow Metab. 2018 Feb 1;38(2):304–16. pmid:28271954
- 32. Aubert A, Pellerin L, Magistretti PJ, Costalat R. A coherent neurobiological framework for functional neuroimaging provided by a model integrating compartmentalized energy metabolism. Proc Natl Acad Sci. 2007 Mar 6;104(10):4188–93. pmid:17360498
- 33. DiNuzzo M, Gili T, Maraviglia B, Giove F. Modeling the contribution of neuron-astrocyte cross talk to slow blood oxygenation level-dependent signal oscillations. J Neurophysiol. 2011 Dec;106(6):3010–8. pmid:21917999
- 34. Lin Y, Stephenson MC, Xin L, Napolitano A, Morris PG. Investigating the metabolic changes due to visual stimulation using functional proton magnetic resonance spectroscopy at 7 T. J Cereb Blood Flow Metab. 2012;32(8):1484–95.
- 35. Bednařík P, Tkáč I, Giove F, Dinuzzo M, Deelchand DK, Emir UE, et al. Neurochemical and BOLD responses during neuronal activation measured in the human visual cortex at 7 Tesla. J Cereb Blood Flow Metab. 2015;35(November 2014):601–10. pmid:25564236
- 36. Schaller B, Mekle R, Xin L, Kunz N, Gruetter R. Net increase of lactate and glutamate concentration in activated human visual cortex detected with magnetic resonance spectroscopy at 7 tesla. J Neurosci Res. 2013;91(8):1076–83. pmid:23378234
- 37. Cedersund G. Conclusions via unique predictions obtained despite unidentifiability–new definitions and a general method. FEBS J. 2012 Sep;279(18):3513–27. pmid:22846178
- 38. Silfvergren O, Simonsson C, Ekstedt M, Lundberg P, Gennemark P, Cedersund G. Digital twin predicting diet response before and after long-term fasting. PLOS Comput Biol. 2022 Sep;18(9):1–26.
- 39. Brännmark C, Lövfors W, Komai AM, Axelsson T, El Hachmane MF, Musovic S, et al. Mathematical modeling of white adipocyte exocytosis predicts adiponectin secretion and quantifies the rates of vesicle exo- and endocytosis. J Biol Chem. 2017 Dec 8;292(49):20032–43. pmid:28972187
- 40. Thompson P, Andersson BJ, Cedersund G. A new method for a priori practical identifiability. bioRxiv; 2022. p. 2022.10.20.511900.
- 41. Sedoglavic A. A Probabilistic Algorithm to Test Local Algebraic Observability in Polynomial Time. J Symb Comput. 2002 May 1;33(5):735–55.
- 42. Kreutz C, Raue A, Kaschek D, Timmer J. Profile likelihood in systems biology. FEBS J. 2013;280(11):2564–71. pmid:23581573
- 43. Maiwald T, Hass H, Steiert B, Vanlier J, Engesser R, Raue A, et al. Driving the Model to Its Limit: Profile Likelihood Based Model Reduction. PLOS ONE. 2016 Sep 2;11(9):e0162366. pmid:27588423
- 44.
Rohatgi A. WebPlotDigitizer [Internet]. Pacifica, CA, USA; 2020. Available from: https://automeris.io/WebPlotDigitizer/.
- 45. Witt ST, Warntjes M, Engström M. Increased fMRI Sensitivity at Equal Data Burden Using Averaged Shifted Echo Acquisition. Front Neurosci. 2016;10:544. pmid:27932947
- 46. Fröhlich F, Weindl D, Schälte Y, Pathirana D, Paszkowski Ł, Lines GT, et al. AMICI: High-Performance Sensitivity Analysis for Large Ordinary Differential Equation Models. Bioinforma Oxf Engl. 2021 Apr 2; pmid:33821950
- 47. Fröhlich F, Kaltenbacher B, Theis FJ, Hasenauer J. Scalable Parameter Estimation for Genome-Scale Biochemical Reaction Networks. PLOS Comput Biol. 2017 Jan 23;13(1):e1005331. pmid:28114351
- 48. Fröhlich F, Theis FJ, Rädler JO, Hasenauer J. Parameter estimation for dynamical systems with discrete events and logical operations. Bioinformatics. 2017 Apr 1;33(7):1049–56. pmid:28040696
- 49. Egea JA, Henriques D, Cokelaer T, Villaverde AF, MacNamara A, Danciu DP, et al. MEIGO: an open-source software suite based on metaheuristics for global optimization in systems biology and bioinformatics. BMC Bioinformatics. 2014 May 10;15(1):136. pmid:24885957
- 50. Sonnay S, Gruetter R, Duarte JMN. How Energy Metabolism Supports Cerebral Function: Insights from 13C Magnetic Resonance Studies In vivo. Front Neurosci. 2017 May 26; 11. pmid:28603480
- 51. Rothman DL, De Feyter HM, de Graaf RA, Mason GF, Behar KL. 13C MRS studies of neuroenergetics and neurotransmitter cycling in humans. NMR Biomed. 2011 Oct;24(8):943–57. pmid:21882281
- 52. Sonnay S, Duarte JM, Just N, Gruetter R. Compartmentalised energy metabolism supporting glutamatergic neurotransmission in response to increased activity in the rat cerebral cortex: A 13C MRS study in vivo at 14.1 T. J Cereb Blood Flow Metab. 2016 May;36(5):928–40.
- 53. Duarte JMN, Lanz B, Gruetter R. Compartmentalized Cerebral Metabolism of [1,6-13C]Glucose Determined by in vivo 13C NMR Spectroscopy at 14.1 T. Front Neuroenergetics. 2011 Jun 6; 3. pmid:21713114
- 54. Schurr A. Glycolysis Paradigm Shift Dictates a Reevaluation of Glucose and Oxygen Metabolic Rates of Activated Neural Tissue. Front Neurosci. 2018 Oct 10;12. pmid:30364172
- 55. Havlicek M, Roebroeck A, Friston K, Gardumi A, Ivanov D, Uludag K. Physiologically informed dynamic causal modeling of fMRI data. NeuroImage. 2015 Nov 15;122:355–72. pmid:26254113
- 56. Griffeth VE, Buxton RB. A theoretical framework for estimating cerebral oxygen metabolism changes using the calibrated-BOLD method: modeling the effects of blood volume distribution, hematocrit, oxygen extraction fraction, and tissue signal properties on the BOLD signal. NeuroImage. 2011 Sep 1;58(1):198–212. pmid:21669292
- 57. Kim JH, Ress D. Arterial impulse model for the BOLD response to brief neural activation. NeuroImage. 2016 Jan 1;124(0 0):394–408. pmid:26363350
- 58.
Lundengård K. Mechanistic modelling—a BOLD response to the fMRI information loss problem [PhD Thesis]. [Linköping, Sweden]: Linköping University; 2017.