|Home | About | Journals | Submit | Contact Us | Français|
Conceived and designed the experiments: KL ME. Performed the experiments: AS FL. Analyzed the data: KL SS AS FL. Contributed reagents/materials/analysis tools: ME GC. Wrote the paper: KL SS FE AS FL GC ME. Conceived and designed the study: ME GC KL FE.
Functional magnetic resonance imaging (fMRI) measures brain activity by detecting the blood-oxygen-level dependent (BOLD) response to neural activity. The BOLD response depends on the neurovascular coupling, which connects cerebral blood flow, cerebral blood volume, and deoxyhemoglobin level to neuronal activity. The exact mechanisms behind this neurovascular coupling are not yet fully investigated. There are at least three different ways in which these mechanisms are being discussed. Firstly, mathematical models involving the so-called Balloon model describes the relation between oxygen metabolism, cerebral blood volume, and cerebral blood flow. However, the Balloon model does not describe cellular and biochemical mechanisms. Secondly, the metabolic feedback hypothesis, which is based on experimental findings on metabolism associated with brain activation, and thirdly, the neurotransmitter feed-forward hypothesis which describes intracellular pathways leading to vasoactive substance release. Both the metabolic feedback and the neurotransmitter feed-forward hypotheses have been extensively studied, but only experimentally. These two hypotheses have never been implemented as mathematical models. Here we investigate these two hypotheses by mechanistic mathematical modeling using a systems biology approach; these methods have been used in biological research for many years but never been applied to the BOLD response in fMRI. In the current work, model structures describing the metabolic feedback and the neurotransmitter feed-forward hypotheses were applied to measured BOLD responses in the visual cortex of 12 healthy volunteers. Evaluating each hypothesis separately shows that neither hypothesis alone can describe the data in a biologically plausible way. However, by adding metabolism to the neurotransmitter feed-forward model structure, we obtained a new model structure which is able to fit the estimation data and successfully predict new, independent validation data. These results open the door to a new type of fMRI analysis that more accurately reflects the true neuronal activity.
Functional magnetic resonance imaging (fMRI) is a widely used technique for measuring brain activity. However, the signal registered by fMRI is not a direct measurement of the neuronal activity in the brain, but it is influenced by the interplay between the metabolism, blood flow and blood volume in the active area. This signal is called the blood-oxygen-level dependent (BOLD) response and occurs when the blood supply to the active area increases in response to neuronal activity. The mechanisms that the cells use to influence the blood supply are not fully known, and therefore it is difficult to know the true neuronal signalling only from inspection of the fMRI signal. In this article, we present a new mathematical model built on the physiological mechanisms thought to underlie the BOLD response. We could successfully fit the model to data and predict the activity caused by new stimuli. By using the validated model we investigated physiological mechanisms that cause different parts of the BOLD response.
Functional magnetic resonance imaging (fMRI) measures brain activity by detecting associated changes in blood oxygenation through the blood-oxygen-level dependent (BOLD) response. The BOLD response is caused by time-dependent changes in deoxyhemoglobin (dHb) concentration . Although the BOLD response reflects neuronal activity through the neurovascular coupling , the mechanisms causing the response are still not fully understood. Here, we investigate these mechanisms by mathematical modeling using a systems biology approach.
The BOLD response endures approximately 15 seconds after a short neural stimulus and it has several characteristic features (Fig 1A) : (i) During the first couple of seconds a shallow undershoot, referred to as the initial dip, is sometimes observed in activated areas of the brain . The initial dip is hypothesized to reflect an increased cerebral metabolic rate of oxygen (CMRO2) that is followed by an increase of dHb content in the blood. (ii) At 6–8 s after the stimulus, the BOLD response peaks as a result of increased cerebral blood volume (CBV) and/or increased cerebral blood flow (CBF). (iii) After the peak, the BOLD response decays and shows a post-peak undershoot before returning to baseline. The mechanisms controlling these processes (i-iii) remain unresolved, and there are at least three different approaches to understand these mechanisms.
The first approach is centered around mathematical modeling. One of the most common approaches is to model the hemodynamic response function (HRF) using the so-called Balloon model , which has been of paramount importance in the development of fMRI image analysis . The Balloon models describe the interplay between CMRO2, CBV, and CBF. The dynamics of these three entities are described in part by purely phenomenological descriptions, such as convolutions with the covariate gamma functions, and in part by physical models e.g. describing the dynamics between CBV and CBF in an expanding balloon. In other words, these Balloon models typically do not incorporate intracellular and biochemical mechanisms involved in cell metabolism or intra-cellular signaling processes related to the BOLD response. Nevertheless, there do exist mathematical models that also incorporate intracellular metabolism , but ultimately even these models explain the actual BOLD response via the Balloon model, which appears as a sub-model in the complete model. Other models, which are not extensions of the Balloon model, describe e.g. spatiotemporal properties of the BOLD response as hemodynamic traveling waves  or oxygen transport in the brain by modeling CBF with a linear flow model and CMRO2 using a gamma function .
The second approach to understanding the BOLD response is centered around the so-called metabolic feedback hypothesis. According to this classical hypothesis (Fig 1B, left), the BOLD response is the result of a tight connection between glucose metabolism and blood flow; when the brain is activated, the neurons consume more energy, resulting in decreased blood glucose and oxygen levels , which trigger a feedback signal increasing CBF to meet metabolic demands. In other words, the metabolic hypothesis is centered around a feedback control to keep glucose level constant.
The third and final approach relevant to this paper is referred to as the neurotransmitter feed-forward hypothesis. This hypothesis is reviewed in e.g , and it is today more actively discussed than the metabolic feedback hypothesis. The neurotransmitter feed-forward hypothesis (Fig 1B, right) suggests sequential feed-forward signaling where neurotransmitters, especially glutamate, cause neurons and astrocytes to activate a chain of intracellular events, involving the release of nitric oxide (NO) or arachidonic acid (AA) metabolites, which in turn control constriction and dilation of the blood vessels. In this way, the feed-forward system “anticipates” the increased need, and goes directly from increased neural activity to increased blood supply.
Of the three approaches mentioned above, only the first involves mathematical modeling and these models are focused mainly on the phenomenological description of the HRF. However, although the metabolic feedback and the neurotransmitter feed-forward hypotheses have been extensively studied through purely experimental approaches, these two hypotheses have never been implemented as mathematical models. It is therefore not known whether the proposed mechanisms of the metabolic feedback and the neurotransmitter feed-forward hypotheses actually would produce a BOLD response or not.
Model-based testing of mechanistic hypotheses has been done in biological research for many years, and has gained increased interest through the rise of systems biology. As mentioned above, mathematical models are already standard when analyzing fMRI, but these models are partially phenomenological and have not been developed to test intracellularly centered hypotheses such as the metabolic feedback and the neurotransmitter feed-forward hypotheses. In contrast, intracellular mechanistic models are the main focus in systems biology, and here hypotheses are formulated as direct representations of the assumed biochemical reactions . This formulation allows for a new type of data analysis, which revolves around two steps: (i) rejections and (ii) uniquely identified core predictions (Fig 1C). This model-based approach provides a more comprehensive, correct, and verifiable analysis, compared to analyses based on inspection and reasoning. In other words, while it sometimes may seem logical to draw a certain conclusion based on visual inspection of some given data, we and others have repeatedly shown that such manual inspections often lead to incorrect, or at the very least incomplete, interpretations of the data [19, 20].
In this paper, we provide a first systems biology analysis of the metabolic feedback and the neurotransmitter feed-forward hypotheses (Fig 1B) with regards to their ability to describe the BOLD response. We show that neither of the two hypotheses alone can satisfactorily describe the response. In contrast, a feed-forward mechanism with added oxygen metabolism can provide a satisfactory explanation of the BOLD response measured in fMRI.
Mechanistic modeling has within systems biology evolved into an iterative process, which alternates between model-based data analysis and the collection of new experimental data. This process is outlined in Fig 1C. In Step 1, existing hypotheses are reformulated into a set of mathematical equations. In this paper, the two main hypotheses are the metabolic feedback and the neurotransmitter feed-forward hypotheses. In Step 2, data are collected, which in this paper correspond to BOLD responses during visual stimulation. The model is then fitted to the data by optimization of the model parameters. Step 3 involves the formulation of biological criteria which might not be found in the collected data, but which have been previously reported in the literature and which the model must fulfil. Step 1–3 compose the hypothesis testing, which includes formulation, fitting and testing of the models. Model testing investigates which sets of equations can and cannot explain the given data, and whether or not these explanations are biologically realistic. This analysis leads to either rejection, which bring the process back to the model formulation step, or acceptance and further analysis. Step 4, minimization, is a simplification of the model where as many states and parameters as possible are removed in order to identify key mechanisms of the system and facilitate computation and overview of the model. In Step 5, further analysis of the explanations consists of the identification of relevant core predictions, i.e. uniquely identified predictions with uncertainty . These predictions may sometimes be suitable for experimental testing, and this leads to collection of new data, Step 6, which in turn leads to the final testing and analysis of the model. Sometimes, interesting model behaviors are discovered in this final step, which might lead back to the hypothesis testing step for further investigation. In this way, systems biology modeling has the potential to be a never-ending cycle, but for each step that is passed, new information about the system is obtained. The iterations end when the model is satisfyingly detailed or no more suitable data can be gathered.
The models herein are formulated using ordinary differential equations (ODEs), which have the following general structure
where x are the states, describing the concentration or amount of various substances; where represents the derivative of the states with respect to time; where f and g are non-linear smooth functions; where px are the parameters used to calculate f, here kinetic rate constants; where u is the input, here the visual stimuli given to the subjects; where x(0) contains the values of the states at time t = 0, and where these values are described by the parameters x0; where are the simulated model outputs corresponding to the measured experimental signals, here the BOLD response; and where py are parameters only appearing in the measurement equations, here scaling parameters. Recall that there are three types of parameters, p, with potentially unknown values,
How these parameters are determined and evaluated is described below. Note that x, u, and y depend on t, but that the notation is dropped unless the time-dependence needs to be especially stated, as in Eq (2). All symbols in Eqs (1), (2) and (3) are vectors.
In the formulation of mechanistic hypotheses into ODEs, there are three levels, which are distinguished using the following notation. The hypothesis x is denoted Mx, where x here has the values m, n, and nm, corresponding to the metabolic feedback, the neurotransmitter feed-forward, and the extended neurotransmitter feed-forward hypothesis hypotheses, respectively. Each of these hypotheses have been implemented using alternative sets of equations, corresponding to further specifications and assumptions, and these alternatives are tagged using additional numbers. Such a set of equations is usually referred to as a model structure . Finally, a model structure is referred to as a model if a set of specific parameter values has been chosen, and these parameters are specified with a final bracket. For example, denotes the 4th model structure implementing the neurotransmitter feed-forward hypothesis, which should be analyzed using the parameters in .
The metabolic model is centered around a feedback control loop, where the rate of blood flow is altered to keep the blood glucose or oxygen levels constant. A schematic overview of this model is shown in Fig 2A. A specific implementation of this hypothesis, Mm3, is plotted in Fig B in S1 Appendix. Fig B in S1 Appendix is an interaction graph, which means that it directly visualizes the interactions included in the model structure, as described in e.g. . More specifically, this means that the non-regulated rates are given by mass-action kinetics, and each differential equation is given by the sum of the in- and out-going reactions, weighted by the stoichiometric matrix. There are three exceptions to this interpretation. First, the blood flow is a variable, not a state, and it affects all four states oxyhemoglobin (oHb), dHb, glucose, and molecular oxygen (O2) by transporting the species in and out of the studied vessel. Second, the stoichiometries of the metabolism in the stimulated and the basal states are different and not specified in Fig B in S1 Appendix. Finally, the delay boxes means that intermediate states have been introduced for e.g. glucose and its influence on the blood flow. This means that the ODE for e.g. oHb is given by
where k1f and k1b are reaction rate constants; where vflow denotes the blood flow, and where [oHb]basal denotes the concentration of oHb in the blood flowing in to the studied area. Similar equations describe states and reactions seen in Fig 2A. All equations and model parameters are specified in detail in S1 Appendix, where also all scripts used for the analyses in the paper can be found. The underlying assumptions of the model are discussed in Section “Assumptions and Limitations”.
The neurotransmitter hypothesis is centered around a feed-forward signaling system, where the regulation of the blood flow is given by the balance of positive and negative regulations. The first model structure corresponding to this hypothesis, Mn1, is depicted in Fig 2B. The ODEs are directly specified by the interaction graph in Fig C in S1 Appendix, using standard and already mentioned conventions. As can be seen in both figures, the input signal triggers release of glutamate into the synaptic cleft. Glutamate increase triggers calcium-permeable non-methyl D aspartate (NMDA) activated channels in neurons and astrocytes to open, letting calcium flow into the cells. In neurons, the calcium influx leads to an increase in the concentration of nitric oxide (NO) that stimulates the production of vasodilating cyclic guanosine monophosphate (cGMP). In the astrocyte, calcium ions increase the production of arachidonic acid (AA), whose metabolites epoxyeicosatrienoic acids (EET), prostaglandins (PG) and hydroxyeicosatetraeonic acid (20-HETE) effect blood vessel radii. The interaction graph in Fig C in S1 Appendix is based on a similar figure in Attwell et al. .
Once the model structure has been formulated (Fig 1C, Step 1 in the modeling workflow) and data has been collected (Step 2, data acquisition described below), the parameters, p, need to be determined. The parameter evaluations are centered on the following cost function
where y(t) are the measured data points at time t; where are the corresponding simulated data points at time t; where N is the number of time-points; and where σ(t) is the measurement uncertainty at time t. The summation in Eq (6) sums the squared and normalized residuals, indicating how far the simulations are from the data y. The additional terms (also referred to as “punishments”, “weights” or “penalties”) are included only when additional requirements are needed. Such requirement might be the presence of an initial dip in the BOLD response (as described in the evaluation criteria in Section “Model evaluation”). If the additional requirements are fulfilled, the additional terms equal zero. If the additional requirements are not fulfilled, the additional terms are increased to force the optimization away from such parameter sets.
In practice, the parameters are determined by optimizing χ2(p) over p, using the function “simannealingSBAO” in the Systems Biology toolbox for Matlab . In other words, the optimal parameters are given by
The final part of the hypothesis testing in the model construction work flow is to check whether the obtained model fulfills the set criteria of the study. In this study the model structure must be able to:
The first criterion is tested using the χ2 test. This test is based on the observation that the cost function in Eq (6) follows a χ2 distribution if the measurement noise follows a normal distribution, with standard deviation σ. Therefore, the resulting cost is compared to the inverse of a cumulative χ2 distribution, where the degrees of freedom are given by the number of data points minus 1. In practice, one decides whether the cost is acceptable or not by comparing the cost with the cut-off value . In our case this means that the cut-off have been 49.8 (χ2, α = 0.05, df = 35) for the estimation data. The second criterion is evaluated via simple simulations, and their fulfillment is ensured by adding these criteria as additional terms in Eq (6).
It is necessary that the model output can fit the data (Criterion 1), but this criterion is not enough to test the underlying mechanism. Therefore, the second and third criteria are added. Most notably, the initial dip is often absent in collected data, although it has been proposed to carry important information about the underlying neuronal activity (discussed further in Discussion: The Initial Dip). The second criterion ensures that the model has the mechanisms required to simulate an initial dip if it is fitted to a data set where it is expressed.
The third criterion is more vaguely phrased, as biological plausibility depends on the mechanisms present in each specific model structure. The third criterion is therefore developed and tested separately for each model structure.
Mechanistic models describing biological systems often become complex with several interacting states and parameters. Such a model, although good for illustrating the biological hypothesis, is computationally heavy and poses difficulties during overview and analysis. A strategy to analyze how the model describes the main mechanisms behind the BOLD response is model minimization (Step 4 in Fig 1C). Model minimization is done by excluding parts of the original model and optimizing it anew to the data. If the reduced model is still able to describe the data, yet another state or a group of states can be excluded. The minimal model has been reached when no more states or reactions can be excluded without the model loosing the ability to fit to the data.
Reducing the number of parameters in the model often makes it harder for the model to fit the data, and thereby the cost of the model will increase compared to the original version of the model. In order to see if the model fit is good enough to compensate for the increased cost, a likelihood ratio test can be performed. If
then the minimized model is not significantly worse at describing the data than the original model, despite its reduced number of parameters. The symbol nmi is the number of parameters in the model structure i and is the lowest cost that the optimization has found for the same model structure.
The final model analysis is referred to as core prediction analysis, seen in Step 5 in Fig 1C. Core prediction analysis is making predictions with uncertainties  and , which are then tested towards validation data (Step 6 in Fig 1C). More precisely, the predictions are analyzed for not only one, but for all parameter sets that pass a χ2 test and thereby can describe the estimation data. In practice, we do not analyze the complete set of acceptable parameters. Instead, we create an approximation of this set, where only a limited number of parameters are studied. This limited set is obtained by saving those acceptable parameters that are encountered during the original optimization and during the subsequent analysis. This corresponds to the hypothesis testing in Fig 1, and to Step 1 in .
Time series data of the BOLD response were collected from the visual cortex of 13 healthy subjects. The subjects were instructed not to consume caffeine, alcohol, or use nicotine on the day of examination. All subjects gave their written informed consent and the study was approved by the Regional Ethical Review Board in Linköping (M74-06).
Individuals failing to fill the criteria set in the standard MR safety screening form were excluded from participating in the study. Maximum allowed translational head movement was limited by a cut-off value of 3.5 mm (the one-dimensional size of one voxel). After reviewing the head movement of each subject, one subject (out of 13) was excluded from further analysis. Thus, 12 subjects (mean age = 23.5 years SD: 4.4, range = 19–35 years) remained in the study. Seven subjects were men and five were women.
The study was designed in the form of two experiments, the intensity and the frequency experiment (Fig 3A and 3B, respectively). These experiments were based on brief visual stimulation using a sparse event-related design in order to isolate individual BOLD responses. The principal visual stimulus consisted of filled white or grey circles shown for 500 ms on a black background. Each stimulus was followed by a randomly jittered intertrial interval (ITI) to reduce adaptation effects. During the ITI a grey focus cross was presented against black background.
Nine subjects performed the intensity experiment first and the frequency experiment last. Four subjects performed the experiments in the opposite order. In both experiments the visual stimuli were presented using video goggles (VisuaStimDigital, Resonance Technology Inc., USA) with a built in screen and correction lenses. The experimental paradigm was presented using the software package SuperLab 4.5 (Cedrus Corporation, San Pedro, CA, USA) using a Windows XP computer. The presentation of the stimuli was randomized using SuperLab’s randomize function, randomizing once per group of participants. Thus, the stimulus onset times did not vary between individuals.
The intensity experiment and the frequency experiment had three stimuli each. In both experiments, the single, bright white circle shown for 500 ms, was used as primary stimulus. The time course of the primary stimulus of each experiment was used as estimation data when the models were fitted.
In the intensity experiment (Fig 3A), the color of the circle was white, light grey, or dark grey. Each stimulus was presented a total of 9 times each with an ITI of 18 to 20 seconds. Thus, the experiment contained 27 trials with a total runtime of approximately 10 minutes.
The frequency experiment (Fig 3B) included the bright white circle on black background presented in three different frequency modes. The first mode consisted of the primary stimulus, which was the same as in the intensity experiment. In the other modes, the stimuli were paired, with an inter-pair interval (IPI) of 1 or 4 seconds, respectively. Each stimulus was presented 8 times each with an ITI of 19 to 21 seconds. The experiment contained 24 trials with a total runtime of approximately 10 minutes.
All experiments were performed with a Philips Ingenia 3 T MR scanner and a 24-channel head coil. BOLD-images were acquired using a gradient echo sequence sensitive for the BOLD contrast using the following parameters: repetition time (TR) = 500 ms, echo time (TE) = 30 ms, resolution = 3.5 mm isotropic, field of view = 224 mm × 196 mm × 35 mm, flip angle = 60 degrees, echo planar imaging (EPI) factor = 29, sense factor = 2.2. Ten axial slices were collected, oriented from the calcarine sulcus to the cingulate gyrus. The number of acquired volumes (number of dynamics) was 1160 per experiment. T1-weighted (T1W) scans were obtained for each individual, as a basis for co-registration of the BOLD-images to high-resolution anatomical images. The following parameters were used for T1W imaging: field of view = 240 mm × 240 mm × 180 mm, voxel size = 0.5 mm × 0.5mm × 0.6 mm, TR = 13 ms, and TE = 6.3 ms.
Images from each subject were preprocessed using SPM8 (www.fil.ion.ucl.ac.uk/spm). All images were re-aligned to the first image in the time series to correct for motion during scanning. Thereafter the images were co-registered to the T1W anatomical reference and normalized to the standard template in MNI (Montreal Neurological Institute) space. The normalized images were smoothed with 7 mm Gaussian kernel to reduce noise and ameliorate inter-subject differences in brain anatomy during the voxel wise group analyses.
BOLD-images from all individual subjects were analyzed using the canonical hemodynamic response function implemented in SPM8. A one-sample t-test was used to identify the peak activation in the visual cortex of the study group. This result was used to guide the extraction of the BOLD responses in individual subjects. The BOLD responses were extracted from the un-smoothed images of each individual in native space using a sphere with radius 5 mm around the individual subject’s peak activation in the visual cortex. The individual peak activation was conjointly estimated for both experiments and all stimulations. MarsBaR  was used to create the spherical masks. BOLD time series data were normalized to baseline by subtracting the signal value of the last time point before the stimulation from the values of the entire BOLD time series corresponding to each stimulus. BOLD responses from each stimulus were first averaged over each individual and thereafter over the group.
Visual stimulation during both the intensity and the frequency experiment elicited significant activation in bilateral primary visual cortex, p < 0.05 family wise error (FWE) corrected for multiple comparisons (Fig 3C). FWE is a Bonferroni-correction applying the random-field theory (RFT) to control the FWE rate by assuming that the data follow certain specified spatial patterns . The Montreal Neurological Institute (MNI) co-ordinates of the activation peaks were: [-2, -96, 4] and [-10, -84, 2] for the intensity and the frequency experiments, respectively.
The mean BOLD response to the primary visual stimuli in both experiments had a characteristic response peak at approximately 6 s after the stimuli (Fig 4). Peak amplitude was 23.5 (2.12% signal change) in the intensity experiment and 19.8 (1.95% signal change) in the frequency experiment. We also observed a post-peak undershoot, but neither of the resulting BOLD responses in any subject displayed a clear initial dip.
When implementing the metabolic feedback model initial attempts were made with a model structure, Mm1, where oxygen levels controlled the blood flow, but such a model structure could neither fulfill the criteria of displaying a post-peak undershoot nor fit the data, as can be seen in Fig A in S1 Appendix. The model structure was therefore rejected at Steps 2 and 3 in the modeling workflow (Fig 1C). If glucose instead of oxygen was selected to control blood flow, it was possible to simulate the shape of the BOLD response by varying the degree of aerobic vs. anaerobic metabolism, as is done in model structure Mm2.
Further tests of Mm2 in the hypothesis testing revealed an important insight regarding a necessary difference between the basal metabolism and the stimulated metabolism. This difference concerns the relationship between oxygen and glucose consumption. If the metabolism is continuously aerobic, i.e. if the ratio of oxygen and glucose consumption is equal during basal state and stimulation, the effect of increased metabolism is oxygen level reduction (, Fig 5B, first three seconds). This oxygen reduction persists until the blood flow is sufficiently up-regulated to normalize oxygen levels (at approximately 10 seconds, Fig 5B). In other words, in this situation, there is an initial dip but no peak.
Conversely, if the stimulated metabolism is completely anaerobic, i.e. if glucose but not oxygen metabolism is increased during stimulation, the effect would be an increase in blood flow to compensate for decreased glucose levels. Oxygen levels will then increase as the increased blood flow brings more oxygen to the capillaries, but this extra oxygen is not metabolized. In this situation (which can be seen in , simulated in Fig 5C) there is a peak but no initial dip.
A combination, , where the stimulated metabolism is adjusted to be more anaerobic than the basal metabolism, both an initial dip and a following peak is produced (Fig 5D). In practice, the proportion of aerobic and anaerobic metabolism was varied by changing the number of oxygen molecules used during glucose metabolism. The parameter sets where the proportion parameters are changed can be seen in section 1.2.3 in S1 Appendix. We did not find any versions of the model that did not require a more anaerobic metabolism during stimulation, and we therefore conclude that, if the metabolic feedback hypothesis is true, the metabolism must be more anaerobic during stimulation compared to the metabolism during basal state.
The metabolic feedback model structure Mm2 assumes glucose regulation of blood flow and a more anaerobic process during stimulation. Mm2 has a statistically acceptable fit to the intensity data primary stimulus according to χ2 goodness-of-fit test (cost = 45.5, cut-off = 49.8). However, as can be seen in Fig 5D, the glucose state decreases very slowly. This means that glucose minimum occurs simultaneously with the BOLD response peak, i.e. simultaneously with the oxygen peak. This behavior contradicts the basic principle of the metabolic hypothesis and is caused by delay states inserted between the neuronal signal and the glucose metabolism in the model structure (dashed square in Fig 5A). These delay states between the stimulus and the glucose metabolism contribute to the shape of the initial dip, but have no biological interpretation.
A new model structure, Mm3, was constructed. In the model structure Mm3 (Fig 5F), the delay states are placed between the metabolism and the blood flow, and represent the action of smooth muscle controlling the radius of the blood vessels. The fit of () and () is shown in Fig 5F and 5G. As can be seen, the model structure Mm3 is incapable of displaying both the initial dip and a post-peak undershoot simultaneously. Furthermore, both fits of Mm3 entail problems with the predicted glucose levels. The graphs depicting glucose levels (Fig 5G and 5I) show that the glucose levels decrease to almost zero (< 5% of the original value) within the first 2 seconds after the stimulus, and in the case of an initial dip (Fig 5I), remain low until the peak of the BOLD response has passed. Even though the exact stimulated glucose dynamics is unknown, we consider this predicted behavior unrealistic.
In summary, the metabolic feedback hypothesis Mm can describe the experimental data, but is still rejected during the hypothesis testing for two reasons. Firstly, it cannot produce both an initial dip and a post-peak undershoot in the same simulation. Secondly, and most importantly, the tested metabolic feedback model structures are not biologically plausible, because they predict depletion of glucose levels and unrealistically fast and/or slow time course of stimulated glucose dynamics.
The neurotransmitter feed-forward model can fit the estimation data from the primary stimulus. As can be seen in Fig 6A, the neurotransmitter feed-forward model clearly displays the typical peak in the BOLD response and the characteristic initial dip and post-peak undershoot. Since the hypothesis presented in Attwell et al.  focuses on the neurovascular control of the blood vessels, there is no metabolism, and under this assumption, the blood flow can be used as a direct proxy for the output, i.e. for the oxygen, dHb, and oHb levels as the output of the model structure Mn1.
The neurotransmitter feed-forward model structure Mn1 is complex with several states, which includes most of the signaling molecules mentioned in the review by Atwell et al. . In order to investigate the key mechanisms of the neurotransmitter feed-forward hypothesis, Mn1 was minimized. Fig 6B shows the minimized neurotransmitter feed-forward model that consists of the simplest combination of states which can still describe the typical BOLD response (fit to data shown in Fig E in S1 Appendix). Mn2 shows that the main mechanism in the neurotransmitter feed-forward model is described by one vasoconstricting and one vasodilating process, shown by the hatched and filled lines respectively in Fig 6C. In order to model a BOLD response with initial dip, peak and post-peak undershoot, it is necessary to have a constricting output term that rises early, but to a lower amplitude than the dilating output term, and falls slowly back to baseline. The dilating output term has a later but quicker rise to high amplitudes, and falls more quickly back to baseline compared to the dilating term.
The neurotransmitter feed-forward model can display all characteristic features of the BOLD response i.e., the initial dip, the peak, and the post-peak undershoot. can also fit the estimation data. However, when investigating the biological mechanisms, we observed that the initial dip, according to the model, is caused by a constricting output term that rises earlier than the dilating term. That is to say, the neurotransmitter feed-forward model predicts vasoconstriction to cause the initial dip. As previous research suggests that the initial dip most probably is related to oxygen metabolism , we reject the neurotransmitter feed-forward model. However, as both the existence of and the mechanisms behind the initial dip is debated, this issue is further addressed in the Discussion.
In summary, the neurotransmitter feed-forward model structure Mn1 can describe the experimental data, and it can produce an initial dip, peak, and post-peak undershoot. However, it is rejected during the hypothesis testing, because the hypothesis is not biologically plausible, as the initial dip is caused by vasocontraction instead of oxygen metabolism.
Results described above show that the increased oxygen metabolism in the metabolic feedback hypothesis can produce an initial dip and that the neurotransmitter feed-forward hypothesis can give a realistic description of the blood flow increase during the BOLD response. Therefore, the neurotransmitter model structure was extended with a metabolic module, Mnm1 (Fig 7). In Mnm1, the neuronal activity increases metabolism in the metabolic module and at the same time triggers glutamate release in the neurotransmitter module. The levels of dHb and oHb are controlled by the metabolic module and the blood flow is controlled by the neurotransmitter feed-forward module. Mnm1 has no feedback control of the blood flow. The output of this final model structure is the ratio of dHb and oHb.
The model structure Mnm1 bridges the gap between cellular action and the regulation of the BOLD response on a vascular level. The model structure Mnm1 can fit data from the primary stimuli in both experiments, see Fig 8A. As can be seen in Fig 8B, the final model has an early and moderate glucose metabolism (the glucose level drops about 5% during the first second). In Fig 8B, it can also be noted that oxygen drops during the first second while in Fig 8C the blood flow is stable during the first seconds, indicating that oxygen metabolism causes the initial dip. The blood flow causing the peak and undershoot is controlled by the neurotransmitter feed-forward module (Fig 8C, peak at 6 s).
The extended model structure Mnm1 represents the system described in Atwell et al.  and contains the key biological elements described there. However, the model structure of Mnm1 can be simplified without loosing the essential mechanisms of the neurovascular coupling. Therefore, Mnm1 was minimized to the minimal model structure Mnm2 (seen in Fig H in S1 Appendix). The minimization was done by removing the parallell pathways controlling vasodilation and other key biological elements, with the criterion that Mnm2 pass both the the likelihood ratio and the χ2 test. If any more such states were removed, the cost of the model fit increased drastically and the model passed neither the likelihood ratio nor the χ2 test. The 49 parameters of Mnm1 were reduced to 27 in Mnm2. The lowest cost of any parameter set found for Mnm1 was 11.2 for the primary stimulus of the intensity experiment and 7.2 for the primary stimulus of the frequency experiment. The lowest cost for the minimal model structure Mnm2 in the same experiments were 27.5 and 18.8, respectively. Cutoff for the likelihood ratio test was 33.9 (df = 22), and thus the model structure Mnm2 passes the likelihood ratio test despite the decreased number of parameters. The minimal model structure Mnm2 does not separate between neurons and astrocytes, but retains only the principle of a dilating and a constricting arm controlling the blood flow. Just as in Mnm1 before the minimization, increased oxygen metabolism causes the initial dip in Mnm2, while the blood flow remains stable during the first seconds after the stimulus.
In order to further test the model structures Mnm, we made core predictions of the BOLD responses to all stimuli in the intensity and frequency experiments according to the model construction work flow described in Fig 1C, Step 5, and validated the core predictions with the new data (Fig 1C, Step 6). Fig 9 shows that the intensity experiment validates the predictions of both the extended model and the minimized model . BOLD responses from the intensity experiment were simulated by changing the amplitude of the input. As the true experimental stimulation intensity for the three types of stimulation (white, light grey and dark grey circles) were not known, only the qualitative behavior of the BOLD response was predicted. Both model structures predict decreased BOLD response peak amplitudes in response to a decreased input signal (Fig 9B and 9C, middle and right panels). Experimental data (Fig 9B and 9C, left panels) confirmed this prediction and indicated that visual stimuli with lower intensity had lower amplitudes of the BOLD response. However, this decrease was not statistically significant for our small sample, p = 0.058 in repeated measures ANOVA (Fig 9, right panels). The white stimulus resulted in the highest mean amplitude, 25.5 au (sd = 7.4) and the dark grey stimulus in the lowest mean amplitude, 20.1 au (sd = 8.2).
In the frequency experiment, the amplitude of the input signal is constant but the stimuli are repeated with 1 s or 4 s IPI. Here, quantitative traits of the data were also predicted. To account for changing basal conditions, the model was re-optimized and the time course from the primary stimulus in the frequency experiment (Fig 10A, left panel) was used as estimation data, yielding and .
For paired stimuli with 1 s IPI, both models predicted a single BOLD response peak with higher amplitude compared to the single stimulus (Fig 10B, middle and right panel). For paired stimuli with 4 s IPI, the extended model Mnm1 predicted a BOLD response doublet with approximately the same amplitude as the single stimulus (Fig 10C, middle panel), while the minimized model Mnm2 predicted that the second peak should be higher than the first (Fig 10C, right panel).
The core predictions were then compared to experimental data. The paired stimulus with 1s IPI resulted in a BOLD response with one peak, while the paired stimulus with 4 s IPI produced a BOLD response doublet where the peaks were of approximately equal hight (Fig 10B and 10C, left panels). The 1 s IPI stimulus resulted in a BOLD response with a significantly greater mean amplitude, 33.1 au (sd = 12.8) than the mean amplitude of the single stimulus response, 22.7 au (sd = 7.8), p = 0.005 in two-tailed paired t-test. Statistically, both models were able to predict the 4 s IPI data, but as can be seen in Fig 10C, the predictions of the extended model Mnm1 has a shape more similar to the validation data compared to Mnm2.
We have presented mathematical modeling of the mechanisms underlying the BOLD response in fMRI, based on the metabolic feedback and the neurotransmitter feed-forward hypotheses, extensively discussed in the literature . These hypotheses describing the fundamental mechanisms behind the BOLD response have, to our knowledge, not been mechanistically modelled before using systems biology approaches. Such approaches provide new tools to evaluate the influence of different hypotheses of the neurovascular coupling causing the BOLD response.
The metabolic feedback model structures Mm1, Mm2, and Mm3 have problems fitting the data (Fig 11). The model structures also have problems with simulating the characteristic traits of the BOLD response, (i) the initial dip, (ii) the peak, and (iii) the post-peak undershoot, at the same time. Further, none of the metabolic feedback model structures could describe a biologically plausible time course of the glucose state (Fig 5). The neurotransmitter feed-forward model structure could describe the data and the characteristic traits, but lacked mechanisms for a correct description of oxygen metabolism as the driving force of the initial dip. Of the model structures tested herein, only the model structures Mnm1 and Mnm2 that combine neurotransmitter control over the blood flow with glucose and oxygen metabolism can fully describe the BOLD response. The model structure Mnm1 could also predict the BOLD responses triggered by stimuli in the intensity and frequency experiments that were not present in the estimation data. Based on these results, we argue that an addition of metabolism to the neurotransmitter feed-forward hypothesis explains necessary mechanisms generating the BOLD response.
We will now discuss the biological interpretation and the underlying assumptions of the model structure Mnm1.
The initial dip is stated as a biological criteria that the model must be able to perform in order to be accepted. The existence of the initial dip is debated ; it has been observed in some studies , but not in others, at least not in all subjects . Several factors can explain the absence of an initial dip in experimental data: (1) The initial dip is reported to be only 1–2% of the baseline signal ; in other words, given the low signal to noise ratio in fMRI, the shallow dip could easily be undetected. (2) The intersubject variability is considerable in fMRI , and thus, the subject selection could be decisive for observing an initial dip or not. (3) The observance of an initial dip could be dependent on the experimental design. For instance, Hu et al.  found that the magnitude of the dip was reduced for brief stimuli; the minimum stimulus duration in their study was 1.5 s, and we used a duration of only 0.5 s, which could possibly explain the absence of the initial dip in our study. Even though the data largely lacked the initial dip, we included the initial dip as a constraint in our models.
The leading hypothesis of the mechanism underlying the initial dip is an uncoupling of the oxygen metabolism from the CBF , where the increased stimulus-induced oxygen metabolism leads to increased dHb levels and following decreased early-phase BOLD response. This hypothesis has been described previously by a model using the gamma variate curve , and is also supported by optical imaging studies showing early stimulus-related dHb increases . Based on these previous studies and the simulations of the metabolic feedback model structure in our work, the oxygen metabolism (but not the feedback) of the metabolic feedback hypothesis is necessary, but not sufficient, to reproduce the shape of the BOLD response. However, the metabolic feedback model structure has difficulties explaining the initial dip in combination with a post-peak undershoot (see Fig 5). The neurotransmitter feed-forward model structure, on the other hand, predicts an initial dip caused by initial vasoconstriction, a prediction that is not supported by any previous data (e.g. ) leading to rejection of that model structure. According to our extended model structure, Mnm1, the initial dip is caused by changes in dHb/oHb ratio due to increased oxygen metabolism occurring during a CBF delay period, as described above.
The most noticeable feature of the BOLD response is the large stimulus-induced peak of the fMRI signal. Originally it was suggested that reduced blood-oxygen levels increased the CBF , and consequently the fMRI signal. This suggestion formed the basis for the model structure Mm1. However, we showed that the oxygen-triggered feedback cannot cause such large overshoot (Fig 5A). Later on, Fox et al.  showed that CBF and CMRO2 correlate strongly in the brain at rest, but not in response to stimuli; they found that CBF increased by 50%, CMRglu by 51%, but CMRO2 only by 5% when the brain is activated. They also concluded that the brain metabolism is aerobic during rest and anaerobic in response to stimuli, to cover the excess energy needed, which is reflected later in the astrocyte-neuron-lactate-shuttle model proposed by Magistretti and Pellerin . Fox et al.  also showed that regional CBF is not driven by oxidative metabolism, but that stimulus-induced CBF is driven by increased glucose demand (see review by Paulson et al. ). Prichard et al.  hypothesized that aerobic glycolysis might be close to its maximum capacity in the resting brain. Therefore, stimulus-induced activity requires quick energy increases via anaerobic glycolysis, causing the uncoupling of glucose and oxygen metabolism and CBF.
In line with previous experimental research , we show that if metabolism is the driving agent for stimulus-induced CBF then blood flow needs to be controlled by glucose, and stimulated metabolism must be partly anaerobic to obtain both an initial dip and a peak in the BOLD response. However, we also show that the metabolic feedback model structures Mm2 and Mm3 overstate predicted glucose reduction in response to stimuli, as they predict almost total depletion of glucose to trigger CBF feedback leading to rejection of the metabolic feedback model structure, Mm. This suggets that the metabolic feedback hypothesis plays a limited, if any, role in the neurovascular response, a conclusion supported by results from Lindauer et al.  and Wolf et al.  who showed that an increased CBF response will still occur even when hemoglobin is fully oxygenated and that CBF remains unchanged at hypoglycemia.
The neurotransmitter feed-forward model structure Mn suggests that CBF is regulated in response to neuronal signaling itself and determined by the intensity and duration of the input signal. In our work, Mn predicts all characteristic features of the BOLD response and has acceptable fit to data. By minimizing the model, we could show that the main mechanism in Mn is a balance between vasoconstriction and vasodilation. One caveat with the neurotransmitter feed-forward model is that it predicts vasoconstriction to cause the initial dip, as discussed above.
According to the final model structure Mnm1, which has neurotransmitter control of the blood flow combined with metabolism of glucose and oxygen, CBF is primarily controlled by neurotransmitters that initiate processes in neurons and astrocytes causing release of vasoactive agents. The final model also predicts glucose response with similar shape as post-stimulus glucose levels measured by optical methods in rat  (Fig 8).
The final characteristic of the BOLD response is the post-peak undershoot. According to our final model structure, Mnm1, the post-peak undershoot is dependent solely on neurotransmitter-triggered changes in CBF. This result is supported by previous literature that suggests that the post-peak undershoot is caused by a post-stimulus CBF undershoot ( reviewed in ). In our data, the post-peak undershoot appears in most of the individual data sets. However, the undershoot seems to become deeper as the peak grows higher in the intensity data sets, and the models tested in the current work cannot predict this behavior. Furthermore, the model cannot describe the deeper undershoot in the 1 s IPI dataset, although it can predict the increased amplitude of the peak.
There are other hypotheses of the mechanisms of the post-peak undershoot that are not investigated here. For example, it is suggested that the post-peak undershoot is caused by slow post-stimulus baseline return of CMRO2-related oxygenation or venous CBV . Mandeville and coworkers  found slow recovery of CBV that matched the post-peak undershoot duration suggesting a biomechanical rather than a metabolic effect . It is worth noting that the balloon model explains the post-peak undershoot as a slow CBV recovery . It has also been suggested that the post-peak undershoot is modulated by post-stimulus neural activity . Future studies incorporating models for CBV changes and/or post-stimulus neural activity might clarify the neurovascular mechanisms behind the post-stimulus undershoot.
The neurotransmitter feed-forward model with metabolism Mnm1 contains several experimentally undetermined parameters, such as glutamate and glucose levels. In future studies, the model parameters can be evaluated and optimized using magnetic resonance spectroscopy (MRS) in combination with BOLD-fMRI. In proton MRS, it is possible to obtain time-dependent variations of the neurotransmitters glutamate and GABA and metabolites such as glucose and lactate . In these recent high-field (7 T) MRS studies, it has been shown that glutamate, GABA, and lactate levels increase during visual stimulation and motor activation, whereas the glucose levels decreases during the same period. If volume was added to the model, a closer estimation of some parameters would also be possible, using experimental values from current literature. This would open the door to prediction of parameter values, in addition to the current predictions of model behaviour.
As with all models, the model structures evaluated in this work depend on a number of underlying assumptions, which in turn limit the conclusions that can be drawn from the results. Nevertheless, such assumptions are essential in order to build a comprehensible model. We are well aware that there are several different hypotheses of the mechanisms behind specific features of the BOLD response, of which some are described above. In this work we chose to focus on two fundamental hypotheses.
In the final model structure, Mnm1, the metabolic and the neurotransmitter module run in parallel. One of the consequences of this structure is that the metabolism is directly controlled by the input signal (see Fig 7). A more physiologically relevant model would be to integrate the metabolism into neurons and astrocytes. That is to say, to model the glycolysis and oxidative metabolism as actually occurring in the neuronal cells in response to stimuli.
Another part of the model structure that lacks physiological details is the action of smooth muscle cells and effects of cortical vessel elasticity. These mechanisms are in our model represented by delay states, which will not accurately represent the possible non-linearities of receptor actions and muscular contraction or relaxation. In addition, the current model does not differentiate between capillaries and arterioles. Recent research has found that cerebral hemodynamics has a spatiotemporal dependence related to the effective blood viscosity and cortical vessel stiffness  and mechanical restrictions on blood vessels depending on cortical depth . A compartmentalized model that takes spatiotemporal hemodynamics into account would provide a physiologically more accurate description of the BOLD response.
Finally, the output signal of the model is oHb/dHb, a simplification suggested by Ogawa et al. . However, the BOLD signal equation provides a more correct description of the output signal.
where S0 is the MR signal at baseline and ΔS is the BOLD signal change with activation. is the difference in transversal relaxation rate between the activated state and baseline. is linearly related to dHb concentration. Following ideas from Davis et al. , several improvements of the BOLD signal description have been published . In the current work, volumes are not included in the model and therefore it is not possible to implement an output dependent on dHb concentration. However, by dividing the model into a tissue compartment and a blood compartment, as has been done in e.g. , a more realistic expression for the output BOLD signal can be obtained.
When comparing models with each other, it is important to keep track of model complexity and watch out for potential problems with overfitting. We approach these issues first by choosing a model framework not designed to be as flexible as possible, but instead based on the actual mechanisms believed to be present in the system. Second, we do visual inspection of the plots comparing data and model fits (Figs (Figs9A9A and 10A). As can be seen in both Fig 9A and 9B, in the time-window 12–18 seconds the mean values in the data show minor fluctuations, which most likely are noise. The model is not following these minor variations, which argues that we do not have problems with overfitting. Nevertheless, in the early time-points (t = 0–3 s), the extended model does an initial dip, stays down a while, and then rises. Since a similar delayed rise can be seen in the data, this could in principle be a sign of overfitting. However, our third approach to test for overfitting—core prediction analysis—argues against that. In the core prediction analysis, we study a representative sub-set of all parameters that describe the data in a statistically acceptable way. In other words, since the core prediction analysis includes both optimal and less optimal parameters, it does not matter if there are some parameters that are overfitted, as long as parameters that are not overfitted are included in the prediction uncertainty analysis. Furthermore, this core prediction analysis shows that all found parameters show an initial dip and a delay (Figs (Figs9A9A and 10A, middle columns), arguing that this property is a necessary consequence of the model structure and the data, i.e. a uniquely identified core prediction. Our fourth approach for checking for unnecessary over-parametrization is model minimization. This punishes for unnecessary complexity in the sense of parameters that can be removed without significantly worsening the agreement with the estimation data (Figs (Figs9A9A and 10A, middle and right columns). Finally, the perhaps most important approach to check for overfitting is to use independent validation data. As can be seen in Figs 9B, 9C and 10B, 10C, both the extended model (middle columns) and the minimized model (right columns), agree with this independent data, to which they have not been fitted. Furthermore, as can be seen in e.g. Fig 10C, the original extended model actually agrees slightly better with the data than the minimized model. All these facts argues that our models—although over-parametrized in the sense that many parameters have non-unique values—still are based on realistic biological mechanisms that capture the main features seen in the data, and not on too flexible model structures that are fitting to a specific noise realization.
Although the BOLD response has been extensively studied and systems biology is a well established method, no one has so far investigated the BOLD response using this type of modeling. In this article, a model based on current physiological hypotheses of the mechanisms behind the BOLD response is presented. The model structures Mnm1 and Mnm2 can describe the time course of the BOLD response in the human visual cortex and correctly predict the response to several variations of the original stimulus not present in the estimation data. In contrast, the individual hypotheses Mm and Mn cannot alone describe the BOLD response in a realistic manner.
Systems biology opens the door to a new type of fMRI analysis, which is firmly based in the physiology of the neurovascular coupling behind the measured signal. Systems biology also gives us the opportunity to obtain information about neural activation beyond what we can measure and may thereby help deepen our understanding of the complex system that is the brain.
It also contains graphs showing the fit of the models Mm1 and Mn2 as well as the simulated glucose behaviour in the model Mnm1.
We thank Dr Suzanne T. Witt for valuable discussions regarding fMRI image analysis and the test subjects for donating their time in the scanner.
This work was supported by the Swedish Research council (2014-6249): http://www.vr.se/; Knut and Alice Wallenbergs foundation, KAW (2013.0076): https://www.wallenberg.com/kaw/; the Research council of Southeast Sweden (FORSS-481691):http://www.fou.nu/is/forss; and Linköping University local funds: http://www.fou.nu/is/lio. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.
All relevant data are within the paper and its Supporting Information files.