You signed in with another tab or window. Reload to refresh your session.You signed out in another tab or window. Reload to refresh your session.You switched accounts on another tab or window. Reload to refresh your session.Dismiss alert
# Collect the graph delays and create a DDEProblem. This will all be zero in this case.
51
+
final_delays =graph_delays(g)
52
+
sim_dur =1000.0# Simulate for 1 second
53
+
prob =DDEProblem(final_system_sys,
54
+
[],
55
+
(0.0, sim_dur),
56
+
constant_lags = final_delays)
57
+
58
+
# Select the algorihm. MethodOfSteps will return the same as Vern7() in this case because there are no non-zero delays, but is required since this is a DDEProblem.
59
+
alg =MethodOfSteps(Vern7())
60
+
sol_dde_no_delays =solve(prob, alg, saveat=1)
61
+
62
+
63
+
# Example of delayed connections
64
+
65
+
# First, recreate the graph to remove previous connections
66
+
g =MetaDiGraph()
67
+
add_blox!.(Ref(g), blox)
68
+
69
+
# Now, add the edges with the specified delays. Again, if you prefer, there's a version using adjacency and delay matrices to assign all at once.
# This is a tutorial for how to run a whole-brain simulation using the Larter-Breakspear model.
2
+
# This approach is the same as that taken in Antal et al. (2023) and van Nieuwenhuizen et al. (2023), with the exception that a different DTI dataset is used.
3
+
# The data used in those papers is available upon request from the authors where it was validated in Endo et al. (2020).
4
+
# As we'd like to make this tutorial as accessible as possible, we'll use a different dataset that is publicly available.
5
+
# The dataset used here is from Rosen and Halgren (2021) and is available at https://zenodo.org/records/10150880.
6
+
7
+
# References:
8
+
# 1. Antal et al. (2023). DOI: 10.48550/arXiv.2303.13746.
9
+
# 2. van Nieuwenhuizen et al. (2023). DOI: 10.1101/2023.05.10.540257.
10
+
# 3. Endo et al. (2020). DOI: 10.3389/fncom.2019.00091.
11
+
# 4. Rosen and Halgren (2021). DOI: 10.1523/ENEURO.0416-20.2020.
12
+
13
+
using Neuroblox # Core functionality
14
+
using CSV, MAT, DataFrames # Import/export data functions
15
+
using MetaGraphs # Set up graph of systems
16
+
using DifferentialEquations # Needed for solver TODO: make a new simulate that can handle system_from_graph
17
+
using StatsBase # Needed for rescaling the DTI matrix
18
+
#using Plots # Only uncomment if you actually want to do plotting, otherwise save yourself the overhead
19
+
20
+
# A note on the data: this is a 360-region parcellation of the brain. For this example, we'll extract and rescale the left hemisphere default mode network.
21
+
# If you want a true rescaling method, you'd need to do some additional corrections (e.g., volume correction) for the number of streamlines computed.
22
+
# For this tutorial though, the rescaling is approximate to just give us a working example.
23
+
# See Endo et al. for more details on how to do a better DTI rescaling.
24
+
25
+
# Load the data
26
+
# As mentioned above, this data is from Rosen and Halgren (2021). You can download the original at https://zenodo.org/records/10150880.
27
+
# For convenience, we have included the average connectivity matrix (log scaled probability from streamline counts) available from that link within this repository.
28
+
data =matread("averageConnectivity_Fpt.mat")
29
+
30
+
# Extract the data
31
+
adj = data["Fpt"]
32
+
adj[findall(isnan, adj)] .=0# Replace NaNs with 0s
33
+
adj = StatsBase.transform(StatsBase.fit(UnitRangeTransform, adj, dims=2), adj) # Equivalent of Matlab rescale function. Resamples to unit range.
34
+
35
+
# For the purpose of this simulation, we want something that will run relatively quickly.
36
+
# In this parcellation, there are 40 regions per hemisphere in the default mode network, so we'll extract the left hemisphere DMN.
37
+
# Original indices with networks are listed in allTables.xlsx at the same download link from the Rosen and Halgren (2021) dataset listed above.
# This is optional but gives you a sense of what the overall connectivity looks like.
47
+
# If you want to do plotting remember to uncomment the Plots import above.
48
+
#heatmap(adj)
49
+
50
+
# Number of regions in this particular dataset
51
+
N =40
52
+
53
+
# Create list of all blocks
54
+
blocks =Vector{LarterBreakspear}(undef, N)
55
+
56
+
for i in1:N
57
+
# Since we're creating an ODESystem inside of the blox, we need to use a symbolic name
58
+
# Why the extra noise in connectivity? The DTI scaling is arbitrary in this demo, so adding stochasticity to this parameter helps things from just immediately synchronizing.
# Create a graph using the blocks and the DTI defined adjacency matrix
63
+
g =MetaDiGraph()
64
+
add_blox!.(Ref(g), blocks)
65
+
create_adjacency_edges!(g, adj)
66
+
67
+
# Create the ODE system
68
+
# This may take a minute, as it is compiling the whole system
69
+
@named sys =system_from_graph(g)
70
+
sys =structural_simplify(sys)
71
+
72
+
# Simulate for 100ms
73
+
sim_dur =1e3
74
+
75
+
# Create random initial conditions because uniform initial conditions are no bueno. There are 3 states per node.
76
+
v₀ =0.9*rand(N) .-0.6
77
+
z₀ =0.9*rand(N) .-0.9
78
+
w₀ =0.4*rand(N) .+0.11
79
+
u₀ = [v₀ z₀ w₀]'[:] # Trick to interleave vectors based on column-major ordering
80
+
81
+
# Create the ODEProblem to run all the final system of equations
82
+
prob =ODEProblem(sys, u₀, (0.0, sim_dur), [])
83
+
84
+
# Run the simulation and save every 2ms
85
+
@time sol =solve(prob, AutoVern7(Rodas4()), saveat=2)
86
+
87
+
# Visualizing the results
88
+
# Again, to use any of these commands be sure to uncomment the Plots import above.
89
+
# This plot shows all of the states computed. As such, it is very messy, but useful to make sure that no out-of-bounds behavior (due to bad connectivity) occured.
90
+
# plot(sol)
91
+
92
+
# More interesting is to choose the plot from a specific region and see the results. Here, we'll plot a specific region's average voltage.
93
+
# First, confirm the region (left orbitofrontal cortex)
94
+
states(sys)[64] # Should give L_OFC₊V(t)
95
+
# Next plot just this variable
96
+
#plot(sol, idxs=(64))
97
+
98
+
# Suppose we want to run a simulation with different initial conditions. We can do that by remaking the problem to avoid having to recompile the entire system.
99
+
# For example, let's say we want to run the simulation with a different initial condition for the voltage.
100
+
v₀ =0.9*rand(N) .-0.61
101
+
z₀ =0.9*rand(N) .-0.89
102
+
w₀ =0.4*rand(N) .+0.1
103
+
u₀ = [v₀ z₀ w₀]'[:] # Trick to interleave vectors based on column-major ordering
# Now that we've confirmed this runs, let's go ahead and do a 10min (600s) simulation.
111
+
# Takes <1min to simulate this run. If you save more frequently it'll take longer to run, but you'll have more data to work with.
112
+
# Remember to update the dt below for the BOLD signal if you don't save at the same frequency.
113
+
# If you find this is taking an *excessively* long time (i.e., more than a couple minutes), you likely happened upon a parameter set that has put you into a very stiff area.
114
+
# In that case, you can try to re-run the simulation with a different initial condition or parameter set.
115
+
# In a real study, you'd allow even these long runs to finish, but for the sake of this tutorial we'll just stop it early.
116
+
sim_dur =6e5
117
+
prob =remake(prob; tspan=(0.0, sim_dur))
118
+
@time sol =solve(prob, AutoVern7(Rodas4()), saveat=2)
119
+
120
+
# Instead of plotting the data, let's save it out to a CSV file for later analysis
121
+
CSV.write("example_output.csv", DataFrame(sol))
122
+
123
+
# You should see a noticeable difference in speed compared to the first time, and notice you save the overhead of structural_simplify.
124
+
125
+
# Now let's do a simulation that creates an fMRI signal out of the neural simulation run above.
126
+
# This example will use the Balloon model with slightly adjusted parameters - see Endo et al. for details.
127
+
128
+
TR =800# Desired repetition time (TR) in ms
129
+
dt =2# Time step of the simulation in ms
130
+
bold =boldsignal_endo_balloon(sol.t/1000, sol[1:3:end, :], TR, dt) # Note this returns the signal in units of TR
131
+
132
+
# Remember that the BOLD signal takes a while to equilibrate, so drop the first 90 seconds of the signal.
133
+
omit_idx =Int(round(90/(TR/1000)))
134
+
bold = bold[:, omit_idx:end]
135
+
136
+
# Plot an example region to get a sense of what the BOLD signal looks like
0 commit comments