The Balloon model was originally designed by Buxton et al. (1998). It describes the changes in the BOLD signal of a tissue region, often called region of interest (ROI), as a function of normalized CBF (f_in). According to this model, the BOLD signal corresponds to the sum of the extravascular and intravascular signal resulting from the normalized total deoxyhemoglobin content (q) and the normalized venous volume fraction (v). The normalized venous volume fraction is described as a balloon that expands with increasing inflow and slowly recovers after a stimulus. The normalized deoxyhemoglobin content is determined by the dynamics of the volume fraction and the blood oxygen extraction fraction (E), whose behavior is based on the oxygen limitation model (Buxton and Frank, 1997).

Friston et al. (2000) extended the Balloon model so that it can be used to simulate BOLD signals using network models. The extension included a neurovascular coupling component that links the normalized CBF of the Balloon model to simulated neuronal activity. Based on this extension, the normalized CBF is modeled as a damped oscillator that is stimulated by neuronal activity. This extension allows the Balloon model to be used to simulate a change in the BOLD signal due to a change in some type of simulated neuronal activity (hereafter, more generally referred to as input signal). In this form, the Balloon model has been used in several studies to compute simulated BOLD signals from network models (Friston et al., 2003; Smith et al., 2011; Deco and Jirsa, 2012; Van Hartevelt et al., 2014; Bennett et al., 2015; Maith et al., 2021). The individual components of the extended Balloon model and their dynamics following a rectangular input signal change are shown in Figure 1A.

Different values for the parameters and even variations of some equations of the model can be found in the literature. We use a version from Stephan et al. (2007) with a non-linear BOLD equation with revised coefficients for our default BOLD monitor. The other versions of Stephan et al. (2007) are also implemented in ANNarchy and available as alternatives. All equations are summarized in Supplementary Section 2. The implementation of the default model in ANNarchy is described in Section 3.4.

In the classic Balloon model, CBF and CMRO2 are tightly coupled. The greater increase in CBF compared to CMRO2 in response to a stimulus is explained by the oxygen limitation model (Buxton and Frank, 1997). This model is based on the assumptions that oxygen coming from the capillaries is completely metabolized in the tissue and that all brain capillaries are perfused at rest. As a consequence, an increase in CMRO2 would only be possible by increasing the transport of oxygen from the capillaries to the tissue, and an increase in CBF would be accompanied by an increase in capillary blood velocity. Because an increase in CBF increases the available oxygen in the capillaries, but also decreases the fraction of oxygen extracted from the capillaries, an increase in CMRO2 (i.e., oxygen transport from the capillaries to the tissue) requires a disproportionate increase in CBF (for further details, see Buxton and Frank, 1997).

However, in recent years, it has been proposed that CBF and CMRO2 are driven in parallel by different sources rather than being tightly coupled (Buxton, 2012, 2021; Buxton et al., 2014). Recently, Buxton (2021) has put forward a new theory, based on the thermodynamics of metabolism, that could explain why CBF needs to increase more than CMRO2 in response to a stimulus and has proposed that CBF and CMRO2 are both driven in parallel in a feed-forward manner. The open question here is by which neural signals CBF and CMRO2 are driven. One suggestion is that CMRO2 is tightly coupled to the energy consumption of neurons, whereas CBF is controlled by vasodilatory signals. These vasodilatory signals are not necessarily coupled to energy consumption and are caused, for example, by activated astrocytes (Buxton, 2012; Howarth et al., 2021). Network models, in which a wide variety of populations can be simulated and manipulated in a controlled manner, may be useful in investigating this question. Therefore, not only the classic Balloon model with tightly coupled CBF and CMRO2 can be used in our BOLD monitor, but also user-defined BOLD models, potentially using more than one input signal from the network model.

We demonstrate how to define BOLD models with multiple input signals for the BOLD monitor in ANNarchy by implementing a modified version of the Balloon model where CBF and CMRO2 are driven in parallel by separate input signals (hereafter referred to as two-input Balloon model). For simplicity, in the two-input Balloon model, we describe both, the normalized CBF and CMRO2, as damped oscillators similar to the normalized CBF in the classic Balloon model version of Friston et al. (2000). Equal input signals elicit responses with equal amplitudes for the normalized CBF and CMRO2. Thus, the coupling between CBF and CMRO2 is determined by the coupling of the two input signals. Figure 1B demonstrates how the individual components of the two-input Balloon model change during stimulation. Compared to the normalized CBF, the normalized CMRO2 responds faster to a changing input signal and without an overshoot or undershoot. The faster response allows for an initial dip in the BOLD signal. For the transformation from normalized CBF and CMRO2 to q and v, we use the Balloon model equations from Buxton et al. (2004). This is a slightly modified version of the classic Balloon model, which additionally considers viscoelastic effects causing the venous volume fraction to lag behind its steady-state relation with the outflow during transient changes. Thus, a post-stimulus undershoot in the BOLD signal is caused by the undershoot of the CBF (based on the damped oscillator modeling approach) as well as by the slow recovery of the venous volume fraction (based on the viscoelastic effects). Finally, the change in the BOLD signal is computed by the non-linear BOLD equation with revised coefficients from Stephan et al. (2007). A more detailed description including the equations of the two-input Balloon model summarized here can be found in Supplementary Section 4.2.

The ANNarchy neural simulator is intended for the simulation of network models at the single-unit level using rate-coded and spiking neuron models. The equation-based interface of ANNarchy allows a flexible and easy implementation of network models by defining equations describing the dynamics of specific neuron types in so called neuron models and equations defining synaptic transmission dynamics (e.g., plasticity) in so called synapse models (Vitay et al., 2015). For efficiency, the model description is transformed into optimized C++ code, optionally using parallel programming frameworks such as openMP for multi-core CPUs or CUDA for GPUs (Dinkelbach et al., 2019). An earlier version of the BOLD monitor in ANNarchy relied on the normalization of pre-synaptic activity and was used in Maith et al. (2021). This implementation was limited to one specific BOLD model and allowed only a few parameter variations, unlike the version presented here. All the simulations in this work use the version 4.7.0.1 of the neural simulator ANNarchy. All references to neurons, populations, synapses, BOLD signals and other neural quantities and data in the following sections refer to simulated values from a network model.

BOLD models, for example the Balloon model (Buxton et al., 1998) or the Davis model (Davis et al., 1998), are based on signals that characterize the dynamics of an entire ROI, such as the change in the normalized CBF or CMRO2 (Figure 2, f_in, r). To combine such a BOLD model with a network model, it is necessary to bridge the gap between these ROI-wide signals and the individual components of the ROI in the network model (e.g., multiple populations, individual neurons). In this section, we will focus on the processes necessary to obtain the input signals for a BOLD model from a ROI that represents part of a network model consisting of multiple populations. Figure 2 shows the general functionality of the BOLD monitor in ANNarchy.

First, the populations of the network model that are part of the ROI for the BOLD computation have to be specified and instantiated. In the example shown in Figure 2, the ROI consists of two populations labeled pop1 and pop2. From the definition of their neuron models, variables must be selected or defined (hereafter referred to as source variables) which will be used to derive the input signals of the BOLD model. In Figure 2, two different source variables are defined: one variable var_CBF that causes the input signal of the CBF (I_CBF) and one variable var_CMRO2 that causes the input signal of the CMRO2 (I_CMRO2). These source variables can correspond to any variables or combinations of variables present in the neuron models (membrane potential, firing rate, etc.).

After defining the ROI and the mapping between source variables in the neuron models and the input variables of the BOLD model, the BOLD monitor implements four preprocessing steps. First, the source variables are averaged over all the neurons for each population of the ROI, resulting in only one signal per population and source variable. This averaging is followed by an optional population-wide normalization that computes the relative deviation of the signal from a baseline value. The baseline corresponds to the mean of the raw averaged source variable signal calculated over a specified initial period. This normalization is useful when deviations from the resting-state are required as input signal in the BOLD model. After the optional normalization, the signals are scaled per population. By default, the signals of each population are scaled based on the ratio between the size of the population and the total number of neurons in the ROI. Thus, the larger a population, the greater its influence on the input variables of the BOLD model. Finally, the population signals are summed across all populations of the ROI, resulting in one input signal for each input variable of the BOLD model.

This section describes a minimal example demonstrating the use of the BOLD monitor in the ANNarchy framework. ANNarchy modules and the BOLD extension must first be imported:

The evaluation of equations is performed with the forward Euler numerical method using a fixed time grid of step dt (in ms):

Two populations, both composed of 100 Izhikevich spiking neurons are then created (line 4, 5). The Izhikevich neuron model is part of the standard models pre-implemented in ANNarchy, with equations and parameters derived from Izhikevich (2003). Initially, the baseline activity in both populations is defined by setting their noise variables to 5.0 (line 7). The term noise refers to an internal variable of the pre-implemented Izhikevich neuron model in ANNarchy which simply determines a baseline current in the membrane potential equation.

To keep the example simple and still have a modulation in the source variable of the BOLD monitor, the baseline activity (the noise variable) is varied during the simulation to mimic the effect of external inputs. The mean-firing rate r of the individual neurons is used as the source variable for the computation of the BOLD signal. As the computation of this value requires an additional overhead, it must be enabled explicitly. The time window for the averaged activity is set to 100 ms:

The BOLD monitor is then created and initialized (line 10–16). The populations in the ROI have to be assigned in the populations argument in form of a list of or a single population (line 11). The desired BOLD model can be optionally defined in the argument bold_model by assigning the corresponding BOLD model object (line 12). The BOLD model can be either one of the built-in BOLD models provided by the module or user-defined as we will demonstrate in Section 3.4. The default BOLD model is the built-in implementation balloon_RN containing the Balloon model with revised coefficients and a non-linear BOLD equation (described in Section 2, implementation shown in Section 3.4).

The mapping between the source variables of the populations (here mean-firing rate r) and the input signals of the BOLD model (referred to as input variables, here I_CBF) has to be defined in the mapping argument by providing a dictionary for each input variable-source variable pair (line 13).

A time window relevant to the normalization of the source variables can be optionally defined (in ms, line 14), whose purpose we explain in Section 4.2. By default, no normalization is performed.

Finally, the variables of the BOLD model which should be recorded during the simulation can be optionally assigned in the recorded_variables argument (line 15) as a string or list of strings. All variables of the BOLD model can be recorded. By default, the output variable of the BOLD model (here the variable BOLD) defined in the BOLD model implementation (see Section 3.4) is recorded.

The C++ code representing the model (network model and BOLD monitor) can now be generated and compiled:

The last part of this section describes a sample simulation to demonstrate the BOLD recording on our simple example. A short simulation period (1,000 ms, line 19) ensures that the network reaches a stable state, which is necessary for a meaningful baseline calculation (required for the normalization outlined in Section 4.2). The recording of BOLD signals is started (line 22) and the simulation is run for 5 s (line 25). After this, the baseline activity (noise variable) of half of the recorded neurons (one population, pop0) is increased for 5 s (lines 26, 27) and afterwards set back to the previous value (line 28, 29).

This leads to an increased mean firing rate in the recorded area and consequently to a BOLD signal response as depicted in Figure 3. The figure shows that the increase of the noise variable in pop0 leads to an increase in the mean-firing rate, which is the source variable for the BOLD monitor (Figure 3A, blue line). This increase of activity results in an increase of the input signal (input variable I_CBF) of the BOLD model depicted in Figure 3B, consequently leading to an increase of the BOLD signal depicted in Figure 3C. After resetting the noise variable, the firing rates of both populations reach again the same level, which reduces the input signal of the BOLD model as well as the resulting BOLD signal.

In the previous example, the default Balloon model (balloon_RN) was used as the BOLD model, but ANNarchy allows users to create their own BOLD model by defining a BoldModel object representing the desired equations. We describe the definition of a BoldModel object using the BOLD model balloon_RN (described in Section 2, applied in Section 3.3) as an example. This BOLD model is implemented as follows:

A BoldModel object requires a parameters argument (line 2), which is a string defining all constants of the BOLD model in a key-value pair notation, i.e., a parameter name on the left and the initialization value on the right side of the assignment operator.

The equations argument (line 10) describes all time-dependent variables, defined either by regular equations or ordinary differential equations evaluated on a fixed time grid. Note that parameters resulting from the combination of other parameters can also be defined here (in this example k_1, k_2, k_3). In the case of a regular equation, the variable name is on the left side and the update performed in each step on the right side. If the update is defined by a differential equation, the left side needs to contain a d[var]dt symbol. To limit the range of values taken by a variable, the min and max keywords can be used. The initial value for variables is 0.0 by default, but it can be changed by providing an init keyword.

The inputs argument (line 30) specifies which input signals are expected by the BOLD model. It consists of a single string or a list of strings. These variables can be accessed in the BOLD model definition by using sum(NAME) in the equations argument, where NAME corresponds to the name of the variable (here I_CBF, line 12).

Finally, in the output argument (line 31), one output variable of the BOLD model is defined, which is automatically recorded by the BOLD monitor. In the following implementation example and all other BOLD models implemented in ANNarchy mentioned in this work, this default output variable corresponds to the BOLD signal (variable BOLD), which is also the default value for the output argument (here only defined for demonstration purposes).

The balloon_RN model is one of the four pre-implemented BOLD models (balloon_RN, balloon_RL, balloon_CN and balloon_CL, Stephan et al., 2007) and therefore does not need to be defined by the user (but its parameters can be changed dynamically). With the BoldModel object, the user can implement new models with the same equation-based interface. For example, a user might want to additionally implement the Davis model (Davis et al., 1998) described by Equation 1 to calculate the change of the BOLD signal ΔBOLD from normalized CBF f and CMRO2 r.

Here, M, α, and β are additional parameters of the Davis model. In the BoldModel above, the normalized CBF is already defined (f_in). Thus, only the calculation of the normalized CMRO2 (r) must be added. This could be done with the term fin·EE0 (see also Buxton et al., 2004). The following code demonstrates how the previous BoldModel could be extended to additionally compute the normalized CMRO2 (r, line 34) and the Davis model BOLD signal (line 35) in the equations argument:

This way, a custom BoldModel is obtained, where the BOLD signal is additionally calculated according to the Davis model and the modified signal (BOLD_Davis) can additionally be recorded. In addition, the parameters of the Davis model would have to be added to the parameters argument, which we have not shown explicitly (but see Supplementary Section 4.3 for a full implementation).

In this section, we implement a simple network model of a cortical microcircuit (hereafter referred to as microcircuit model) to further demonstrate use cases of the BOLD monitor. The microcircuit model consists of a population of excitatory neurons and a population of inhibitory interneurons. As neuron models, we use a regular spiking cortical neuron model for the excitatory population (corE) and a fast-spiking cortical interneuron model for the inhibitory population (corI), both introduced in Izhikevich (2007). The two populations receive excitatory inputs from another population whose neurons randomly emit spikes such that their inter-spike intervals correspond to a Poisson process (hereafter referred to as Poisson neurons). The structure of the microcircuit model is shown in Figure 4A. The projections of the microcircuit model include feed-forward excitation (Poisson neurons → corE), feed-forward inhibition (Poisson neurons → corI → corE), and feedback inhibition (corE → corI → corE). The ratio between excitatory neurons and inhibitory interneurons is 4:1, as found, for example, for the visual cortex (Beaulieu et al., 1992; Potjans and Diesmann, 2014). The equations and parameters of the microcircuit model can be found in Supplementary Section 3.

Each neuron receives synaptic input from 10 random neurons in the pre-synaptic population for each projection. Following our previous modeling approaches (Baladron et al., 2019; Goenner et al., 2021; Maith et al., 2021), we model synaptic inputs as conductance-based synapses in our neuron models. Therefore, the synaptic currents (which drive the membrane potential of the neurons) are proportional to the product of a voltage difference (between the membrane potential and the synaptic reversal potential) and a conductance variable representing the spike input of the corresponding synapse (see Supplementary Section 3 for equations). We model only two different types of conductance-based synapses, excitatory synapses (AMPA) and inhibitory synapses (GABA). The conductance variables of the synapses are instantaneously increased by a fixed value (by the weight of the synaptic connection) for each incoming spike and otherwise decay exponentially to zero with a time constant of 10 ms.

A conductance greater than zero causes a synaptic current that drives the membrane potential toward the reversal potential associated with the synapse (0 mV for AMPA synapses and −90 mV for GABA synapses). All synaptic weights are drawn from a log-normal distribution and scaled by a factor for each projection during model initialization. The weights and scaling factors were optimized to replicate distributions from excitatory post-synaptic potentials (Song et al., 2005) and firing rates (Buzsáki and Mizuseki, 2014) with the microcircuit model (see Figure 4B). Further details about obtaining the distributions and optimizing the parameters can be found in Supplementary Section 3.3.

Although the use of neuron models mimicking spiking patterns of real cortical neurons and tuning the parameters to replicate experimental data can provide more realistic network models (see e.g., Humphries et al., 2006; Günay et al., 2008; Pospischil et al., 2008; Goenner et al., 2021), the microcircuit model presented here only aims at demonstrating the application of the BOLD monitor and not at replicating any particular experimental data. To keep the model simple, we chose a network model with two spiking populations and multiple excitatory and inhibitory projections. No particular functional processing takes place in this microcircuit model, as it consists of only two small homogeneous populations, the connectivity is random and synaptic plasticity, important neurotransmitters such as NMDA, the effect of neuromodulators and potential dynamic changes in activity were not taken into account during construction. However, the applicability of the BOLD monitor to larger-scale network models is demonstrated in Section 4.4.

We first demonstrate the effect of baseline normalization in the BOLD monitor using the microcircuit model. To do so, we simulate a brief stimulus presentation corresponding to studies of the event-based BOLD response (Glover, 1999; Serences, 2004) by briefly increasing the mean firing rate of the Poisson neurons and meanwhile recording the BOLD response.

All simulations start with an initialization period of 2 s to allow the microcircuit model to enter its steady-state. After that, the recordings are started. A 10-s resting-period is simulated, after which the mean firing rates of the Poisson neurons are increased by a factor of five for 100 ms (hereafter referred to as stimulus pulse). Finally, another post-stimulus resting-period is simulated until a total simulation time of 25 s. This procedure is performed for 40 different random microcircuit model initializations (each with different seeds producing different synaptic contacts, weights, and mean firing rates of Poisson neurons). Additionally, we run 40 simulations without a stimulus pulse, in which only a 25 s resting-period is simulated for comparison.

The BOLD response is recorded simultaneously using two differently initialized BOLD monitors. Both BOLD monitors use the default BOLD model (balloon_RN) shown in Section 3.4 and determine the BOLD signal of the ROI which comprises both the corE and corI populations. The source variable for the BOLD monitor is the synaptic activity of the neurons normalized by the number of afferent connections, which has already been used and described in Maith et al. (2021). The key difference between the two BOLD monitors is the baseline normalization. One BOLD monitor uses no baseline normalization and the other BOLD monitor uses a baseline computed over the first 5 s after the 2-s initialization period.

Figure 5 shows the recorded variables of the BOLD model: the input variable (I_CBF) of the BOLD model and the resulting BOLD signal. Although the response of the microcircuit model to the stimulus pulse can be clearly seen in the I_CBF of both BOLD monitors, an important difference is that the I_CBF of the BOLD monitor without baseline normalization has an offset greater than zero, while the I_CBF with baseline normalization fluctuates around zero. It is also noticeable that I_CBF with baseline normalization is zero in the first 5 s. This is because the input variable for the BOLD model is not calculated during the time in which the baseline for normalization is determined. In the normalized CBF signal and the BOLD signal, one can clearly see the effect of baseline normalization on the Balloon model dynamics. The response to the stimulus pulse is much more pronounced for the BOLD monitor with baseline normalization. Without baseline normalization, the normalized CBF signal and BOLD signal at rest have an offset greater than zero, whereas with baseline normalization, the signals fluctuate around one and zero, respectively.

The CBF and BOLD signals are defined in the Balloon model relative to their value at rest (normalized CBF and relative change of BOLD). Therefore, the normalized CBF signal or the BOLD signal should only deviate from one or zero, respectively, when the underlying system deviates from its resting-state. For models with resting-state activity, we recommend using the baseline normalization of the BOLD monitor when using the Balloon model.

One important motivation for developing the BOLD monitor is to provide a simple way to flexibly adjust both the source variables and the BOLD model itself. In Section 3.4, we have already shown how to implement a user-defined BOLD model. Here, we also want to show the possibility to use different variables of the neurons as source variables. The underlying neural processes influencing CBF and CMRO2, and thus ultimately the BOLD signal, are still rather unclear (Howarth et al., 2021). Many different hypotheses and modeling approaches can be found in the literature (Smith et al., 2011; Van Hartevelt et al., 2014; Bennett et al., 2015; Heikkinen et al., 2015; Schmidt et al., 2018). The flexible BOLD monitor in ANNarchy allows us to easily create and compare BOLD models implementing different hypotheses on spiking or rate-coded network models. In this section, we demonstrate this by implementing six different hypotheses using our microcircuit model. For each hypothesis, we add a different BOLD monitor to the microcircuit model, each with different source variables. The six different BOLD monitors are summarized in Table 1. The source code for adding them to the microcircuit model can be found in Supplementary Section 4. Note that the simulated BOLD signals are not compared with experimental data, so we do not make any statements about the validity of the hypotheses. Such an analysis would require an extensive underlying network model, tailored to the brain region under investigation.

We again use the stimulus pulse simulation from Section 4.2 to compare the different BOLD signal responses (see Figure 6). The first three hypotheses are based on previous studies that used the classic Balloon model. Thus, we also use the classic Balloon model (BOLD model balloon_RN) for the BOLD calculation, which includes a single CBF-driving input signal (see Figure 1A) whose source variable we vary for each hypothesis. The first hypothesis we implement is that the CBF or the BOLD signal is driven by the total synaptic activity of the neurons (as in Van Hartevelt et al., 2014; Schmidt et al., 2018; Maith et al., 2021). To implement this, we use the normalized synaptic activity as the source variable of the BOLD monitor (BOLD monitor A), as previously in Section 4.2. The second hypothesis we implement is that the CBF or the BOLD signal is driven only by the excitatory (glutamatergic) synaptic activity (similar to Heikkinen et al., 2015). For this, we use the conductance variable of the excitatory synapses of the neurons as source variable for the BOLD monitor (BOLD monitor B). The third hypothesis is that the CBF or the BOLD signal is driven by the neuronal output of the neurons, for example, the mean firing rate (as in Smith et al., 2011; Bennett et al., 2015). Thus, for this BOLD monitor (BOLD monitor C), we use the mean firing rate of the neurons as source variable, as in Section 3.3. Figures 6A–C shows that the normalized CBF and BOLD responses vary for these three BOLD monitors with different source variables. The response based on the mean firing rates (Figure 6C) is the strongest, because the firing rates change more relatively to the resting-state compared to the two other source variables. However, the shape of the responses is almost identical.

As mentioned in Section 2.2, it has also been proposed that the CBF and CMRO2 are driven in parallel in a feed-forward manner. Therefore, for the following three BOLD monitors, we use the two-input Balloon model defined in Section 2.2 (balloon_two_inputs), which requires two input signals (I_CBF, I_CMRO2, see Figure 1). The source variables used to obtain I_CBF and I_CMRO2 can be freely chosen from the neuron models of the corE and corI populations.

The first hypothesis considering CBF and CMRO2 being driven in parallel proposes that the CMRO2 is driven only by excitatory synaptic processes and that the CBF is driven by both excitatory and inhibitory synaptic processes (Buxton, 2012, 2021). To implement this hypothesis in BOLD monitor D, we define the current caused by AMPA synapses (I_AMPA) as the CMRO2-driving source variable, and the sum of I_AMPA and the current caused by GABA synapses (I_GABA) as the CBF-driving source variable. These source variables have to be additionally defined in the neuron models of the neurons of the corE and corI populations (see Supplementary Section 4.1).

The next hypothesis is similar, but additionally states that in inhibitory interneurons, energy consumption, and thus CMRO2, is driven by neuronal output rather than synaptic input (in contrast to excitatory neurons) (Howarth et al., 2021). To implement this in BOLD monitor E, the mean firing rate of the neurons rather than I_AMPA is defined as the CMRO2-driving source variable for the inhibitory interneurons of the corI population. For the excitatory neurons of the corE population, the same source variables are used as in the previous BOLD monitor.

Figures 6D,E show that the normalized CBF, CMRO2, and BOLD (relative change) responses of these two BOLD monitors are significantly different from the previous ones (with a single input). The BOLD signal shows a much stronger initial dip as CMRO2 increases much faster than CBF. There is little difference between the responses of the two BOLD monitors. The CMRO2 of BOLD monitor E is slightly lower because the firing rate of the inhibitory interneurons increases less than their synaptic current caused by AMPA synapses. However, because the inhibitory interneurons only contribute one-fifth to the input signal of the BOLD monitor (due to the ratio between corE and corI sizes), there is only a small difference from BOLD monitor D to E.

In the last BOLD monitor (Figure 6F), we use almost the same source variables as in BOLD monitor D. We only introduce an additional non-linear operation for the source variable driving CMRO2 by defining the current caused by the AMPA synapses, to the power of one third, as the source variable (instead of the current itself). As a result, energy consumption or CMRO2 no longer increases linearly with the current. Thus, we are still basically following the same general hypothesis (CMRO2 driven by AMPA synaptic processes, CBF driven by AMPA and GABA synaptic processes), but assuming different mathematical relationships for CMRO2. This change causes CMRO2 to increase much less due to the stimulus pulse, as shown in Figure 6F. Thus, the initial dip in the BOLD response is smaller than for the BOLD monitors D & E.

In summary, the BOLD monitor allows users to determine the BOLD signal based on individually chosen source variables. Without much effort, we can define different source variables and even compare different BOLD models (e.g., a model with two input variables). With the classic Balloon model, the BOLD response for our microcircuit model hardly differs for different source variables. Since all variables in the microcircuit model increase similarly in response to the stimulus pulse, the BOLD response also looks similar and only differs in amplitude. When driving CBF and CMRO2 in parallel with different source variables, the choice of the source variable is much more important, because the relationship between them critically affects the shape of the BOLD response not only the amplitude. Nevertheless, the effect of changing the source variable on the resulting BOLD signal may be different for other underlying network models with different dynamics of the different variables (e.g., synaptic currents, mean firing rate, etc.), even when the classic Balloon model is used.

In a second experiment, we perform a simulation with sustained stimulation (longer stimulus pulse) with our six different BOLD monitors. The firing rate of the Poisson neurons is increased by a factor of 1.2 for 20 s. Like in the stimulus pulse simulations, the responses of the first three BOLD monitors (A-C) are very similar and only differ in amplitude (Figure 7, left). The CBF or BOLD responses show a slight initial overshoot, then reach a plateau, and finally, show a slight post-stimulus undershoot. The three BOLD monitor variants with two input signals (D-F) again show significant differences from the three BOLD monitors using the classic Balloon model. The BOLD monitors D and E showed almost identical responses consisting of an initial undershoot a negative plateau and a post-stimulus over- and undershoot (for results of BOLD monitor E see Supplementary Figure 1). It is particularly noticeable that the plateau is negative for the BOLD monitors D & E but not for BOLD monitor F because only for BOLD monitor F, the CBF increases more than the CMRO2. This again illustrates how critical the choice of source variables is when CBF and CMRO2 are driven in parallel by them.

In this section, we study the additional computational time introduced by the BOLD monitor (hereafter referred to as computational overhead). We use a scaled version of the microcircuit model described in Section 4.1, by incrementally increasing the number of neurons for the populations and leaving the number of synaptic inputs for a neuron fixed to 10 connections (from 10 different neurons of the pre-synaptic population) per projection. Table 2 shows an overview of the total number of neurons and connections for each network model instance.

Figure 8 depicts the single thread computational time in seconds as a function of the number of recorded neurons with (blue line) and without (orange line) BOLD recording. For each configuration, we performed 10 runs, each simulating 25 s biological time and we measured the elapsed real time with the Python time module. The relative standard deviation was in the range of 0.55% to 2.53% which is too small to be depicted meaningfully in the graph and was therefore omitted.

For all simulated configurations, the computational time with and without BOLD recording is globally similar (between 1% and 8% of overhead depending on the model's size). The relative computational overhead (visualized as gray bars) is larger for small network models but shrinks when the model size increases. Therefore, if network models get more complex, in the sense of number of neurons, complexity of neuron models and the number of connections, one can expect that the share of the computational overhead will shrink accordingly. Overall, the computational time is dominated by the complexity of the network model and the BOLD recording plays a minor role, especially in complex network models.