Mechanistic model for human brain metabolism and its connection to the neurovascular coupling

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.


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][4][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 (CMRO 2 ) [9][10][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][17][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, CMRO 2 , etc. To remedy this, later mathematical modelling included some mechanistic insight and resulted in the Balloon model, developed by Buxton et al. 1998 [19][20][21]. The Balloon model describes the interplay between neuronal  [16][17][18][19]21,[23][24][25][26][27][28][29][30][31][32][33][55][56][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. https://doi.org/10.1371/journal.pcbi.1010798.g001 PLOS COMPUTATIONAL BIOLOGY activity, CBF, CMRO 2 , 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][27][28][29]. Furthermore, there are works that integrate the metabolic response with NVC, but these works rely on the phenomenological Balloon model [30][31][32][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).
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).

Model formulation
The model presented herein is formulated using ordinary differential equations (ODEs) and can in general be denoted as _ x ¼ f ðxðtÞ; y; uðtÞ; tÞ ð1AÞ where x is a vector describing the model states e.g., the concentrations of the different metabolites, signalling intermediates, etc; where _ x 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 x 0 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]. 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. https://doi.org/10.1371/journal.pcbi.1010798.g002

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 χ 2test. 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.  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.  [23,24,38,39] formulated as: where t is the time; t on is the time stimulation is turned on and t off 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.: where Stimulus is a state that represent the effect of the visual stimuli; k stim1 describes the rate of stimulus activation and k stim2 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 (Gluc t ). The rate of this reaction is described by the variable V Gluc , 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 (V glycolysis ) is assumed to have a basal rate and a Stimulus activated rate, both governed by mass action kinetics. With these assumptions, the ODE for Gluc t is given by: where Gluc t describes the concentration of glucose in the neuronal tissue and where the kinetic rate parameters for inflow and consumption of Gluc t are described by k max gluc c and k max gluc t , respectively.
In the model, we assume that the plasma glucose concentration close to the activation site can be modified i.e., that Gluc c 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.
where Gluc c 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 V glycolysis , 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 V LDH , 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 V PC 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 V TCA1 ), which describes the first part of the TCA cycle ( Fig 1C). The rate expression for V TCA1 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: where Pyr describe the concentrations of pyruvate; k max PO , k max Pyr , and k max Pyr2 are kinetic rate parameters; and where K M Pyr , K M OAA and K M Pyr2 are Michaelis-Menten constants. Note that k max PO , k max Pyr2 describe the maximum rate of V TCA1 and V PC , respectively, while k max Pyr 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 (V LDH ) and a clearance term (V clear1 ). The differential equation for lactate is given by: where Lac is the concentration of lactate and k 1 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: V TCA1 ). Further, the remaining reactions that returns OG to OAA is also simplified into a single reaction TCA2 (reaction rate: V TCA2 ). 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 V GOT . With these reactions, the full differential expression for OAA is thus described as:

PLOS COMPUTATIONAL BIOLOGY
where OAA is the concentration of oxaloacetic acid in the tissue and where k max OG1 , and k max OAA are the kinetic rate parameters describing the reaction rates of V TCA2 and V GOT , 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: V xm ). The full ODE for OG is thus given by: where OG is the concentration of 2-oxoglutarate, and k max OG2 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 V GOT . The subsequent metabolism of aspartate is implemented as a simple reaction denoted V clear2 where Asp is the concentration of aspartate and k max Asp 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: V GS ). Reversely, glutamine is converted back to glutamate via phosphate-activated glutaminase (rate expression: V PAG ). Finally, glutamate is decreased by V clear3 , 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: where Glut is the concentrations of glutamate; k max Gln and k max Glut1 are the kinetic rate parameters for the forward and reverse conversion of glutamine to glutamate, respectively; and k max Glut2 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 V GS and V GLUT rate expressions, i.e.
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.  [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:

Connection between NVC
where k basalMet , k prop1 and k prop2 are kinetic rate parameters; where Glucose A and O 2A are the respective concentrations of glucose and oxygen in the neuronal tissue; where Delay M represents an intermediate state of the stimulus effect on the metabolism. The variable Glucose A describes the same quantity as the state Gluc t described in the removed Eq (4). The reason for this change of name in the combined model is that Glucose A 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, Gluc t , 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: 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%: where X i is one of the observable properties Lac, Glut, Gluc t , Asp, or BOLD; where ssX i is the corresponding steady-state value, and where ky i is the i th scaling parameter. 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 10 5 samples, and all parameter sets found acceptable according to a χ 2 -test was collected. The uncertainty is then determined by the confidence interval: 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.

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]. [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.

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.

Model implementation
Model implementations and simulations were done in the "Advanced Multi-language Interface to CVODES and IDAS" (AMICI) toolbox [46][47][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.

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).

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).

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).

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.

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.

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 13 C 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-tonoise 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][27][28][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 CMRO 2 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 (V TCA1 and V TCA2 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 BOLDsignal (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 CMRO 2 . 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 BOLDresponses 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 BOLDresponse 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.

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 CMRO 2 , 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.  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]. (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. (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. (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. (ZIP) 32. Aubert