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

@@ -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"
 

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

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
 

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`.

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))
 

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)

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

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())))