From f92d833acbaad60405c3fcc1d86e481abf1cb3e3 Mon Sep 17 00:00:00 2001 From: haris organtzidis Date: Fri, 11 Oct 2024 17:55:33 +0300 Subject: [PATCH] More docs shape-up (#452) * add Downloads and SparseArrays deps * shape up resting state tutorial * rename `get_bloxs_parts -> get_parts` and add `get_componets` * use new `add_edge!` syntax * fix get_components * fix typo * remove guesses from algebraic states * give download access to image stimulus file --- docs/Project.toml | 2 + docs/src/tutorials/neural_assembly.jl | 4 +- docs/src/tutorials/parkinsons.jl | 70 ++++++++++++--------------- docs/src/tutorials/resting_state.jl | 53 ++++++++++++-------- src/Neurographs.jl | 2 +- src/blox/blox_utilities.jl | 11 +++-- src/blox/connections.jl | 8 +-- src/blox/discrete.jl | 3 +- 8 files changed, 82 insertions(+), 71 deletions(-) diff --git a/docs/Project.toml b/docs/Project.toml index aa9d4b3e..6004c527 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -4,6 +4,7 @@ CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0" DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" +Downloads = "f43a241f-c20a-4ad4-852c-f6b1247861c6" HypothesisTests = "09f84164-cd44-5f33-b23f-e6b0d136a0d5" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306" @@ -12,6 +13,7 @@ Neuroblox = "769b91e5-4c60-41ee-bfae-153c84203cb2" OrderedCollections = "bac558e1-5e72-5ebc-8fee-abe8a469f55d" OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" +SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" StochasticDiffEq = "789caeaf-c7a9-5a7d-9973-96adeb23e2a0" diff --git a/docs/src/tutorials/neural_assembly.jl b/docs/src/tutorials/neural_assembly.jl index f8985ae4..30c5e87b 100644 --- a/docs/src/tutorials/neural_assembly.jl +++ b/docs/src/tutorials/neural_assembly.jl @@ -22,6 +22,7 @@ using Random ## for generating random variables using CairoMakie ## for customized plotting recipies for blox using CSV ## to read data from CSV files using DataFrames ## to format the data into DataFrames +using Downloads ## to download image stimuli files # define a single excitatory neuron 'blox' with steady input current I_bg = 0.5 microA/cm2 nn1 = HHNeuronExciBlox(name=Symbol("nrn1"), I_bg=0.5) @@ -234,8 +235,7 @@ global_namespace=:g # create an image source block which takes image data from a .csv file and gives input to visual cortex -fn = joinpath(@__DIR__, "../data/image_example.csv") ## image data file -image_set = CSV.read(fn, DataFrame) ## reading data into DataFrame format +image_set = CSV.read(Downloads.download("raw.githubusercontent.com/Neuroblox/NeurobloxDocsHost/refs/heads/main/data/image_example.csv"), DataFrame) ## reading data into DataFrame format image_sample = 2 ## set which image to input (from 1 to 1000) ## define stimulus source blox diff --git a/docs/src/tutorials/parkinsons.jl b/docs/src/tutorials/parkinsons.jl index 426245de..09826173 100644 --- a/docs/src/tutorials/parkinsons.jl +++ b/docs/src/tutorials/parkinsons.jl @@ -43,49 +43,41 @@ using CairoMakie @named PY = JansenRit(cortical=true) # default parameters cortical Jansen Rit blox @named II = JansenRit(τ=2.0*τ_factor, H=60/τ_factor, λ=5, r=5) -blox = [Str, GPE, STN, GPI, Th, EI, PY, II] - - #Here, we've created eight Jansen-Rit neural masses representing different brain regions involved in Parkinson's disease. The `τ_factor` is used to convert time units from seconds (as in the original paper) to milliseconds (Neuroblox's default time unit). # ## Building the Circuit -# Now, let's create a graph representing our brain circuit: - -g = MetaDiGraph() -add_blox!.(Ref(g), blox) - -# We've created a MetaDiGraph and added our neural masses as nodes. Next, we'll define the connections between these nodes based on the known anatomy of the basal ganglia-thalamocortical circuit. - -# ModelingToolkit allows us to create parameters that can be passed into the equations symbolically: - -# As an alternative to creating edges with an adjacency matrix (as shown in the previous example), here we demonstrate a different approach by adding edges one by one. In this case, we set the connections specified in Table 2 of Liu et al. [1], although we only implement a subset of the nodes and edges to describe a simplified version of the model: - -## Define connection strength parameters -params = @parameters C_Cor=60 C_BG_Th=60 C_Cor_BG_Th=5 C_BG_Th_Cor=5 - -add_edge!(g, 2, 1, Dict(:weight => -0.5*C_BG_Th)) -add_edge!(g, 2, 2, Dict(:weight => -0.5*C_BG_Th)) -add_edge!(g, 2, 3, Dict(:weight => C_BG_Th)) -add_edge!(g, 3, 2, Dict(:weight => -0.5*C_BG_Th)) -add_edge!(g, 3, 7, Dict(:weight => C_Cor_BG_Th)) -add_edge!(g, 4, 2, Dict(:weight => -0.5*C_BG_Th)) -add_edge!(g, 4, 3, Dict(:weight => C_BG_Th)) -add_edge!(g, 5, 4, Dict(:weight => -0.5*C_BG_Th)) -add_edge!(g, 6, 5, Dict(:weight => C_BG_Th_Cor)) -add_edge!(g, 6, 7, Dict(:weight => 6*C_Cor)) -add_edge!(g, 7, 6, Dict(:weight => 4.8*C_Cor)) -add_edge!(g, 7, 8, Dict(:weight => -1.5*C_Cor)) -add_edge!(g, 8, 7, Dict(:weight => 1.5*C_Cor)) -add_edge!(g, 8, 8, Dict(:weight => 3.3*C_Cor)) -add_edge!(g,1,1,:weight, -0.5*C_BG_Th) -add_edge!(g,1,2,:weight, C_BG_Th) -add_edge!(g,2,1,:weight, -0.5*C_BG_Th) -add_edge!(g,2,5,:weight, C_Cor_BG_Th) -add_edge!(g,3,1,:weight, -0.5*C_BG_Th) -add_edge!(g,3,2,:weight, C_BG_Th) -add_edge!(g,4,3,:weight, -0.5*C_BG_Th) -add_edge!(g,4,4,:weight, C_BG_Th_Cor) +# Now, let's create a graph representing our brain circuit. The nodes on this graph are the neural mass models defined aboe and the edges are the connections between the nodes based on the known anatomy of the basal ganglia-thalamocortical circuit. +# As an alternative to creating edges with an adjacency matrix, here we demonstrate a different approach by adding edges one by one. In this case, we set the connections specified in Table 2 of Liu et al. [1], although we only implement a subset of the nodes and edges to describe a simplified version of the model. +# Our connections share some common parameters which we define here, both as symbols and values, and use them as expressions for the weight of each connection : + +g = MetaDiGraph() ## define an empty graph + +params = @parameters C_Cor=60 C_BG_Th=60 C_Cor_BG_Th=5 C_BG_Th_Cor=5 # define common connection parameters + +## Create connections +add_edge!(g, GPE => Str; weight = -0.5*C_BG_Th) +add_edge!(g, GPE => GPE; weight = -0.5*C_BG_Th) +add_edge!(g, GPE => STN; weight = C_BG_Th) +add_edge!(g, STN => GPE; weight = -0.5*C_BG_Th) +add_edge!(g, STN => PY; weight = C_Cor_BG_Th) +add_edge!(g, GPI => GPE; weight = -0.5*C_BG_Th) +add_edge!(g, GPI => STN; weight = C_BG_Th) +add_edge!(g, Th => GPI; weight = -0.5*C_BG_Th) +add_edge!(g, EI => Th; weight = C_BG_Th_Cor) +add_edge!(g, EI => PY; weight = 6*C_Cor) +add_edge!(g, PY => EI; weight = 4.8*C_Cor) +add_edge!(g, PY => II; weight = -1.5*C_Cor) +add_edge!(g, II => PY; weight = 1.5*C_Cor) +add_edge!(g, II => II; weight = 3.3*C_Cor) +add_edge!(g, Str => Str; weight = -0.5*C_BG_Th) +add_edge!(g, Str => GPE; weight = C_BG_Th) +add_edge!(g, GPE => Str; weight = -0.5*C_BG_Th) +add_edge!(g, GPE => Th; weight = C_Cor_BG_Th) +add_edge!(g, STN => Str; weight = -0.5*C_BG_Th) +add_edge!(g, STN => GPE; weight = C_BG_Th) +add_edge!(g, GPI => STN; weight = -0.5*C_BG_Th) +add_edge!(g, GPI => GPI; weight = C_BG_Th_Cor) # ## Creating the Model diff --git a/docs/src/tutorials/resting_state.jl b/docs/src/tutorials/resting_state.jl index d4165b2f..a9881770 100644 --- a/docs/src/tutorials/resting_state.jl +++ b/docs/src/tutorials/resting_state.jl @@ -19,33 +19,41 @@ using Random using CairoMakie using Statistics using HypothesisTests +using Downloads +using SparseArrays -## read connection matrix from file -weights = CSV.read("../data/weights.csv",DataFrame) +## 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 = Array(weights) ## transform the weights into a matrix -N_bloxs = size(wm)[1] ## number of blox components +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) + +# 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 -blocks = [Generic2dOscillator(name=Symbol(region_names[i]),bn=sqrt(5e-4)) for i in 1:N_bloxs] +blox = [Generic2dOscillator(name=Symbol(region_names[i]),bn=sqrt(5e-4)) for i in 1:N_bloxs] ## add neural mass models to Graph and connect using the connection matrix g = MetaDiGraph() -add_blox!.(Ref(g), blocks) +add_blox!.(Ref(g), blox) create_adjacency_edges!(g, wm) @named sys = system_from_graph(g); # 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. -prob = SDEProblem(sys,rand(-2:0.1:4,76*2), (0.0, 6e5), []) -sol = solve(prob, EulerHeun(), dt=0.5, saveat=5) +tspan = (0.0, 6e5) +prob = SDEProblem(sys,rand(-2:0.1:4,76*2), tspan, []) +sol = solve(prob, EulerHeun(), dt=0.5, saveat=5); -# Let us plot the voltage potential of the first couple of components +# Let us now plot the voltage potential of the first couple of components -v1 = voltage_timeseries(blocks[1], sol) -v2 = voltage_timeseries(blocks[2], sol) +v1 = voltage_timeseries(blox[1], sol) +v2 = voltage_timeseries(blox[2], sol) fig = Figure() ax = Axis(fig[1,1]; xlabel = "time (ms)", ylabel = "Potential") @@ -56,17 +64,20 @@ fig # To evaluate the connectivity of our simulated resting state network, we calculate the statistically significant correlations -cs = [] -for i in 1:Int((length(sol.t)-1)/1000)-1 - solv = Array(sol[1:2:end,(i-1)*1000+1:(i*1000)])' - push!(cs,cor(solv)) +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) + +for (i, t) in enumerate(time_steps[1:end-1]) + V = voltage_timeseries(blox, sol; ts = t:(t + step_sz)) + cs[:,:,i] = cor(V) end -css = stack(cs) -p = zeros(76,76) -for i in 1:76 - for j in 1:76 - p[i,j] = pvalue(OneSampleTTest(css[i,j,:])) +p = zeros(N_bloxs, N_bloxs) +for i in 1:N_bloxs + for j in 1:N_bloxs + p[i,j] = pvalue(OneSampleTTest(cs[i,j,:])) end end @@ -75,3 +86,5 @@ heatmap(log10.(p) .* (p .< 0.05)) heatmap(wm) # Fig.: Connection Adjacency Matrix that was used to connect the neural mass models + +# Notice how the correlation heatmap qualitatively matches the sparsity structure that we printed above with `wm_sparse`. \ No newline at end of file diff --git a/src/Neurographs.jl b/src/Neurographs.jl index d53292a4..4cb4ee7e 100644 --- a/src/Neurographs.jl +++ b/src/Neurographs.jl @@ -57,7 +57,7 @@ end get_sys(g::MetaDiGraph) = get_sys.(get_bloxs(g)) get_dynamics_bloxs(blox) = [blox] -get_dynamics_bloxs(blox::Union{CompositeBlox, AbstractComponent}) = get_blox_parts(blox) +get_dynamics_bloxs(blox::Union{CompositeBlox}) = get_parts(blox) flatten_graph(g::MetaDiGraph) = mapreduce(get_dynamics_bloxs, vcat, get_bloxs(g)) diff --git a/src/blox/blox_utilities.jl b/src/blox/blox_utilities.jl index 431d1d12..2a9daca2 100644 --- a/src/blox/blox_utilities.jl +++ b/src/blox/blox_utilities.jl @@ -57,11 +57,16 @@ function get_neurons(vn::AbstractVector{<:AbstractBlox}) mapreduce(x -> get_neurons(x), vcat, vn) end +get_parts(blox::CompositeBlox) = blox.parts +get_parts(blox) = blox + +get_components(blox::Union{CompositeBlox, Vector{<:AbstractBlox}}) = mapreduce(x -> get_components(x), vcat, get_parts(blox)) +get_components(blox) = blox + get_neuron_color(n::AbstractExciNeuronBlox) = "blue" get_neuron_color(n::AbstractInhNeuronBlox) = "red" get_neuron_color(n::AbstractNeuronBlox) = "black" - function get_discrete_parts(b::Union{AbstractComponent, CompositeBlox}) mapreduce(x -> get_discrete_parts(x), vcat, b.parts) end @@ -179,8 +184,6 @@ get_weight_learning_rules(bc::BloxConnector) = bc.learning_rules get_weight_learning_rules(blox::Union{CompositeBlox, AbstractComponent}) = (get_weight_learning_rules ∘ get_connector)(blox) get_weight_learning_rules(blox) = Dict{Num, AbstractLearningRule}() -get_blox_parts(blox::Union{CompositeBlox, AbstractComponent}) = blox.parts - function get_weight(kwargs, name_blox1, name_blox2) get(kwargs, :weight) do @warn "Connection weight from $name_blox1 to $name_blox2 is not specified. Assuming weight=1" @@ -506,7 +509,7 @@ end function state_timeseries(cb::Union{CompositeBlox, AbstractVector{<:AbstractBlox}}, sol::SciMLBase.AbstractSolution, state::String; ts=nothing) - neurons = get_neurons(cb) + neurons = get_components(cb) state_names = map(neuron -> Symbol(namespaced_nameof(neuron), "₊", state), neurons) if isnothing(ts) diff --git a/src/blox/connections.jl b/src/blox/connections.jl index d7dc5307..583f6ddc 100644 --- a/src/blox/connections.jl +++ b/src/blox/connections.jl @@ -280,8 +280,8 @@ function (bc::BloxConnector)( bloxin::CanonicalMicroCircuitBlox; kwargs... ) - sysparts_out = get_blox_parts(bloxout) - sysparts_in = get_blox_parts(bloxin) + sysparts_out = get_parts(bloxout) + sysparts_in = get_parts(bloxin) wm = get_weightmatrix(kwargs, namespaced_nameof(bloxin), namespaced_nameof(bloxout)) @@ -297,7 +297,7 @@ function (bc::BloxConnector)( kwargs... ) - sysparts_in = get_blox_parts(bloxin) + sysparts_in = get_parts(bloxin) bc(bloxout, sysparts_in[1]; kwargs...) end @@ -307,7 +307,7 @@ function (bc::BloxConnector)( bloxin::ObserverBlox; kwargs... ) - sysparts_out = get_blox_parts(bloxout) + sysparts_out = get_parts(bloxout) bc(sysparts_out[2], bloxin; kwargs...) end diff --git a/src/blox/discrete.jl b/src/blox/discrete.jl index aca753f4..1e3aca72 100644 --- a/src/blox/discrete.jl +++ b/src/blox/discrete.jl @@ -79,7 +79,8 @@ struct SNc <: AbstractModulator function SNc(; name, namespace=nothing, κ_DA=1, N_time_blocks=5, DA_reward=10, λ_DA=0.33, t_event=90.0) sts = @variables R(t) R_(t) - ps = @parameters κ=κ_DA λ=λ_DA jcn=0 [input=true] jcn_=0.0 #HACK: jcn_ stores the value of jcn at time t_event that can be accessed after the simulation + ps = @parameters κ=κ_DA λ=λ_DA jcn=0 [input=true] jcn_=0 #HACK: jcn_ stores the value of jcn at time t_event that can be accessed after the simulation + eqs = [ R ~ min(κ, κ/(λ*jcn + sqrt(eps()))), R_ ~ min(κ, κ/(λ*jcn_ + sqrt(eps())))