Resting state simulation using neural mass models

This tutorial will introduce you to simulating resting state brain dynamics using Neuroblox. We will be using the FitzHugh-Nagumo model as a building block. The FitzHugh-Nagumo model is described by the follwoing equations:

In this tutorial you will learn to:

  • Build a whole-brain resting state network from an empirical connectivity matrix.
  • Use SDEProblem to simulate noisy (stochastic) neural dynamics with the EulerHeun solver.
  • Extract voltage timeseries from multiple blox using voltage_timeseries.
  • Compute functional connectivity from simulated data and compare it to the structural connectivity.

\[ \begin{align} \dot{V} &= d \, \tau (-f V^3 + e V^2 + \alpha W - \gamma I_{c} + \sigma w(t) ) \\ \dot{W} &= \dfrac{d}{\tau}\,\,(b V - \beta W + a + \sigma w(t) ) \end{align}\]

We start by building the resting state circuit from individual Generic2dOscillator Blox

using Neuroblox
using CSV
using DataFrames
using StochasticDiffEq
using Random
using CairoMakie
using Statistics
using HypothesisTests
using Downloads
using SparseArrays

Random.seed!(123); # Set a fixed RNG seed for reproducible results

# download and read connection matrix from a file
weights = CSV.read(Downloads.download("raw.githubusercontent.com/Neuroblox/NeurobloxDocsHost/refs/heads/main/data/weights.csv"), DataFrame)
region_names = names(weights)

wm = Matrix(weights) ## transform the weights into a matrix
N_bloxs = size(wm)[1]; ## number of blox components

You can visualize the sparsity structure by converting the weights matrix into a sparse matrix

wm_sparse = SparseMatrixCSC(wm)
76×76 SparseArrays.SparseMatrixCSC{Float64, Int64} with 1560 stored entries:
⎡⡛⣌⡀⣑⣀⣉⡐⣒⣒⣁⣙⣒⣗⣀⢀⠧⡜⡇⡀⠀⢄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎤
⎢⣭⢬⣯⢍⠠⣬⢑⣛⣿⠆⣭⣥⣅⡀⡉⡅⣬⠅⡅⠀⠀⠑⢄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠤⣬⢥⡤⢡⣭⢨⣥⣬⡄⣬⢥⣮⢡⣭⠀⡤⡄⡄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠰⢾⣷⢰⠆⣿⣵⣿⣧⡆⢈⢱⣾⡿⠿⡷⢶⡦⡇⠀⠀⠀⠀⠀⠀⠁⢀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠐⢹⠿⠿⠀⠿⠹⠟⠿⠷⠾⠿⠿⠿⠳⠯⠼⠍⠇⠀⠀⠀⠀⠀⠀⠀⠀⠑⢄⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⣻⢻⡻⣗⠘⣿⢰⣆⣿⡇⣿⣷⣷⣶⢶⣷⣿⡇⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠁⠄⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠹⢹⠉⠹⠈⣿⣾⣷⣿⡷⢹⠹⣱⣶⣾⡏⠹⠁⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠐⠄⠀⠀⠀⠀⠀⎥
⎢⠤⡴⠄⢤⠃⢽⢻⡧⡽⠄⢤⢄⡌⠑⠛⣤⣤⡤⡄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠁⠀⠀⠀⠀⎥
⎢⠿⡯⠮⠵⠀⠯⢸⡆⠾⠆⠽⠿⡗⠀⠠⡿⢾⢇⡆⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢀⠀⎥
⎢⠀⠉⠀⠉⠀⠉⠈⠉⠉⠁⠈⠉⠉⠉⠀⠉⠉⠉⠁⣤⠠⠀⢄⠀⠤⢀⣀⣀⠄⢤⣀⣄⠀⠀⡄⢠⡄⠁⎥
⎢⠀⠑⢄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠦⠳⡦⠶⠒⠲⢆⣶⣶⡒⠶⠖⠗⠒⠴⠍⠣⠇⠆⎥
⎢⠀⠀⠀⠑⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠛⠹⠟⠑⠌⠿⠰⠖⠻⠁⠻⠟⡳⠆⠦⠃⠛⠁⠃⎥
⎢⠀⠀⠀⠀⠀⠀⠄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢉⣻⣝⢋⡘⣿⢜⣿⡟⡃⠻⢝⣻⣼⣿⣄⣋⡃⡇⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠐⢄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢈⢽⣿⣼⠁⣿⢿⣿⣿⣇⣰⣼⣿⣯⣍⡯⢹⠯⡇⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠑⠄⠀⠀⠀⠀⠀⠀⠀⠀⣤⣼⣭⣍⢠⣭⢈⡁⣭⡍⣭⣍⣍⣉⣈⣍⣭⡅⡅⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠁⢀⠀⠀⠀⠀⠀⠀⢾⢼⠮⢷⠠⣿⣸⣗⣿⣇⢿⢿⢟⣛⣹⡿⢿⠇⡇⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠁⠄⠀⠀⠀⠀⠈⢘⠀⠈⡄⢿⣿⡟⢿⠋⠘⠈⠺⢟⣿⠃⠈⠀⠃⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⣭⡯⡡⢝⠀⡽⢸⡋⣫⡁⢽⣵⣇⠀⠀⣿⣻⡏⡃⎥
⎣⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠐⠄⠉⠯⠉⠭⠀⠭⠸⠧⠭⠅⠩⠭⠧⠤⠈⠯⠽⠵⠇⎦

After the connectivity structure, it's time to define the neural mass components of our model and then use the weight matrix to connect them together into our final system.

# create an array of neural mass models
@graph g begin
    @nodes blox = [Generic2dOscillator(bn=sqrt(5e-4)) for i in 1:N_bloxs]
end

# add connections between the neurons according to `wm`

create_adjacency_edges!(g, wm)

To solve the system, we first create an Stochastic Differential Equation Problem and then solve it using a EulerHeun solver. The solution is saved every 0.5 ms. The unit of time in Neuroblox is 1 ms.

tspan = (0.0, 6e5)

u0map = Iterators.flatten([(osc.V => rand(-2:0.1:4), osc.W => rand(-2:0.1:4)) for osc ∈ blox])

prob = SDEProblem(g, u0map, tspan, [])
# `EulerHeun()` is a fixed-step stochastic solver suitable for Stratonovich SDEs. It requires an
# explicit `dt` (step size). For Itô SDEs you would use `EM()` (Euler-Maruyama). See the
# DifferentialEquations.jl solver guide: https://docs.sciml.ai/DiffEqDocs/stable/solvers/sde_solve/
sol = solve(prob, EulerHeun(), dt=0.5, saveat=5);

Let us now plot the voltage potential of the first couple of components. voltage_timeseries(blox, sol) is a convenience function that extracts the membrane voltage (V) state variable of a given blox from the solution object and returns it as a plain Julia vector. It is equivalent to calling state_timeseries(blox, sol, "V"), but is more readable when you specifically want voltage. Here we call it for the first two oscillators in our network.

v1 = voltage_timeseries(blox[1], sol)
v2 = voltage_timeseries(blox[2], sol)

fig = Figure()
ax = Axis(fig[1,1]; xlabel = "time (ms)", ylabel = "Potential")
lines!(ax, sol.t, v1)
lines!(ax, sol.t, v2)
xlims!(ax, (0, 1000)) ## limit the x-axis to the first second of simulation
fig
Example block output

To evaluate the connectivity of our simulated resting state network, we calculate the statistically significant correlations

step_sz = 1000
time_steps = range(1, length(sol.t); step = step_sz)

cs = Array{Float64}(undef, N_bloxs, N_bloxs, length(time_steps)-1)

# First we get the voltage timeseries of all blox in time bins
for (i, t) in enumerate(time_steps[1:end-1])
    # Get the voltage timeseries of all blox within a time window in a matrix
    V = voltage_timeseries(blox, sol; ts = t:(t + step_sz))
    # Calculate the correlation matrix of all timeseries per time window
    cs[:,:,i] = cor(V)
end

p = zeros(N_bloxs, N_bloxs)
for i in 1:N_bloxs
    for j in 1:N_bloxs
        # Perform a t-test on the voltage correlations across all time windows
        p[i,j] = pvalue(OneSampleTTest(cs[i,j,:]))
    end
end

fig = Figure()
ax = [
    Axis(fig[1,1]; title = "Correlation p-values", aspect = AxisAspect(1)),
    Axis(fig[1,2]; title = "Adjacency matrix", aspect = AxisAspect(1))
]
# Use logical indexing to only plot the statistically significan p-values (p .< 0.05)
heatmap!(ax[1], log10.(p) .* (p .< 0.05))

heatmap!(ax[2], wm)
fig
Example block output

Notice how the correlation heatmap qualitatively matches the sparsity structure that we printed above with wm_sparse.

This page was generated using Literate.jl.