BOLD monitoring - Balloon model

Download JupyterNotebook Download JupyterNotebook

#!pip install ANNarchy

This notebook demonstrates the usage of the BOLD monitoring extension. It has to be explicitly imported:

import numpy as np

import ANNarchy as ann
from ANNarchy.extensions.bold import *
ANNarchy 4.8 (4.8.2) on darwin (posix).

Background

ANNarchy pre-implements some model variants of the BOLD models presented in Stephan et al. (2007) which are variations of the Balloon Model originally designed by Buxton et al. (1998). The four balloon variants model pre-implemented in ANNarchy follow the naming scheme of Stephan et al. (2007). The model name starts with balloon_ followed by either C (i. e. classical coefficient) or R (i. e. revised coefficient) and then followed by either N which means non-linear BOLD equation or L which means linear bold equation.

We only provide here the equations without much explanations, for more details please refer to the literature:

Buxton, R. B., Wong, E. C., and Frank, L. R. (1998). Dynamics of blood flow and oxygenation changes during brain activation: the balloon model. Magnetic resonance in medicine 39, 855–864. doi:10.1002/mrm.1910390602

Friston et al. (2000). Nonlinear responses in fMRI: the balloon model, volterra kernels, and other hemodynamics. NeuroImage 12, 466–477

Buxton et al. (2004). Modeling the hemodynamic response to brain activation. Neuroimage 23, S220–S233. doi:10.1016/j.neuroimage.2004.07.013

Stephan et al. (2007). Comparing hemodynamic models with DCM. Neuroimage 38, 387–401. doi:10.1016/j.neuroimage.2007.07.040

Maith et al. (2021). A computational model-based analysis of basal ganglia pathway changes in Parkinson’s disease inferred from resting-state fMRI. European Journal of Neuroscience. 2021; 53: 2278– 2295. doi:10.1111/ejn.14868

Maith et al. (2022). BOLD Monitoring in the Neural Simulator ANNarchy. Frontiers in Neuroinformatics 16. doi:10.3389/fninf.2022.790966.

Single input Balloon model

This script shows a simple example for the application of one of the default models (the balloon_RN model) on a simple network. The balloon_RN model is described by the following equations:

\tau \cdot \frac{ds}{dt} = \phi \cdot I_\text{CBF} - \kappa \cdot s - \gamma \cdot (f_{in} - 1)

\frac{df_{in}}{dt} = s

E = 1 - (1 - E_{0})^{ \frac{1}{f_{in}} }

\tau \cdot \frac{dv}{dt} = \frac{1}{\tau_0} \cdot (f_{in} - f_{out})

f_{out} = v^{\frac{1}{\alpha}}

\tau \cdot \frac{dq}{dt} = \frac{1}{\tau_0} \cdot ( f_{in} \cdot \frac{E}{E_0} - \frac{q}{v} \cdot f_{out} )

with revised coefficients and non-linear bold equation:

k_1 = 4.3 \, v_0 \, E_0 \, \text{TE}

k_2 = \epsilon \, r_0 \, E_0 \, \text{TE}

k_3 = 1 - \epsilon

BOLD = v_0 \cdot ( k_1 \cdot (1-q) + k_2 \cdot (1 - \dfrac{q}{v}) + k_3 \cdot (1 - v) )

There are two important variables in that model: BOLD which is the output of the model and I_CBF which is the input signal, reflecting the mean firing rate of the input populations.

As the BOLD model reflects the change of activity levels in the recorded area, we will implement the following experiment. We record from two populations with 100 Izhikevich neurons each. After a short period of time we raise the activity level of one population for a fixed time window. For simplicity, we do not use connections between the neurons but simply increase the noise term in the Izhikevich model. This should lead to a noticeable BOLD signal.

Populations

We first create two populations of Izhikevich neurons:

ann.clear()

pop0 = ann.Population(100, neuron=ann.Izhikevich)
pop1 = ann.Population(100, neuron=ann.Izhikevich)

As we will not have any connections between the neurons, we need to increase the noise to create some baseline activity:

# Set noise to create some baseline activity
pop0.noise = 5.0; pop1.noise = 5.0

The mean firing rate of the neurons will be used as an input to the BOLD model. This need to be activated explicitly as the computation of this value is quite expensive. In our example, the mean firing rate should be computed across a time window of 100 ms. The resulting value will be stored in the r attribute of the populations and can be easily recorded.

# Compute mean firing rate in Hz on 100ms window
pop0.compute_firing_rate(window=100.0)
pop1.compute_firing_rate(window=100.0)

# Record the mean firing rate
mon_pop0 = ann.Monitor(pop0, ["r"], start=False)
mon_pop1 = ann.Monitor(pop1, ["r"], start=False)

BOLD Monitor definition

The BOLD monitor expects a list of populations which we want to record (in our case pop0 and pop1). A BOLD model should be specified, here we take balloon_RN which is the default. We then specify the mapping from the population variable that should influence the BOLD signal, in our case the mean firing rate r, to the input variable of the BOLD model I_CBF.

The mean firing rate is normalized over a time window of 2000 ms and we record both input and output signals:

m_bold = BoldMonitor(
    
    populations = [pop0, pop1], # recorded populations
    
    bold_model = balloon_RN(), # BOLD model to use (default is balloon_RN)
    
    mapping = {'I_CBF': 'r'}, # mapping from pop.r to I_CBF
    
    normalize_input = 2000,  # time window to compute baseline.
    
    recorded_variables = ["I_CBF", "BOLD"]  # variables to be recorded
)

Now we can compile and initialize the network:

ann.compile()
Compiling ...  OK

Simulation

We first simulate 1 second biological time to ensure that the network reaches a stable firing rate:

ann.simulate(1000)

We then enable the recording of all monitors:

mon_pop0.start()
mon_pop1.start()
m_bold.start()

We simulate for 5 seconds with lower noise and we increase the noise in pop0 for 5 seconds before decreasing it again:

# We manipulate the noise for the half of the neurons
ann.simulate(5000)      # 5s with low noise
pop0.noise = 7.5
ann.simulate(5000)      # 5s with higher noise (one population)
pop0.noise = 5
ann.simulate(10000)     # 10s with low noise

# Retrieve the recordings
mean_fr1 = np.mean(mon_pop0.get("r"), axis=1)
mean_fr2 = np.mean(mon_pop1.get("r"), axis=1)

input_data = m_bold.get("I_CBF")
bold_data = m_bold.get("BOLD")

Evaluation

We can now plot:

  • the mean firing rate in the input populations.
  • the recorded activity I which serves as an input to the BOLD model.
  • the resulting BOLD signal.
import matplotlib.pyplot as plt

plt.figure(figsize=(15, 5))

# mean firing rate
ax1 = plt.subplot(131)
ax1.plot(mean_fr1, label="pop0")
ax1.plot(mean_fr2, label="pop1")
plt.legend()
ax1.set_ylabel("Average mean firing rate [Hz]")

# BOLD input signal
ax2 = plt.subplot(132)
ax2.plot(input_data)
ax2.set_ylabel("BOLD input I_CBF")

# BOLD output signal
ax3 = plt.subplot(133)
ax3.plot(bold_data*100.0)
ax3.set_ylabel("BOLD [%]")

# x-axis labels as seconds
for ax in [ax1, ax2, ax3]:
    ax.set_xticks(np.arange(0,21,2)*1000)
    ax.set_xticklabels(np.arange(0,21,2))
    ax.set_xlabel("time [s]")

plt.show()