updated CMC code (#474)

* switched from symbolic arrays to list of symbolic parameters. Needed to be able to set tunable flag for individual parameters of the array. Added untune! utility function to set tunable flag to false on list elements. * added new approach to set tunable flag for the connection matrix also to the tutorial * use generate_weight_param also in StimulusBlox and ObserverBlox connectors and set standard connection weight to non-tunable * fixed and cleaned code * minor changes * update MTK compat * added Project.toml to test * fixed MTK compat once more * added more packages to test environment needed in GraphDynamics tests * remove test deps and unnecessary deps * add Test to tect project * remove Test dep from main Project * add missing test deps * add Statistics dep * removed [extras] from Project.toml * added [extras] again. --------- Co-authored-by: haris organtzidis <[email protected]>
Neuroblox · Nov 13, 2024 · ebb3421 · ebb3421
1 parent 40d4a9e
commit ebb3421
Show file tree

Hide file tree

Showing 7 changed files with 62 additions and 59 deletions.
diff --git a/Project.toml b/Project.toml
@@ -34,7 +34,6 @@ Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 StochasticDiffEq = "789caeaf-c7a9-5a7d-9973-96adeb23e2a0"
 Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
-Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 ToeplitzMatrices = "c751599d-da0a-543b-9d20-d0a503d91d24"
 
 [weakdeps]
@@ -57,7 +56,7 @@ GraphDynamics = "0.2"
 Graphs = "1"
 Interpolations = "0.14, 0.15"
 MetaGraphs = "0.7"
-ModelingToolkit = "9.46.0 - 9.49"
+ModelingToolkit = "9.46.0 - 9.50.0"
 ModelingToolkitStandardLibrary = "2"
 MuladdMacro = "0.2"
 OrderedCollections = "1.6.3"
@@ -84,6 +83,3 @@ SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 StochasticDiffEq = "789caeaf-c7a9-5a7d-9973-96adeb23e2a0"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 Turing = "fce5fe82-541a-59a6-adf8-730c64b5f9a0"
-
-[targets]
-test = ["DelayDiffEq", "LinearAlgebra", "ForwardDiff", "Distributions", "Random", "SparseArrays", "MAT", "OrdinaryDiffEq", "StochasticDiffEq", "SafeTestsets", "Turing", "Test"]
diff --git a/docs/src/tutorials/spectralDCM.jl b/docs/src/tutorials/spectralDCM.jl
@@ -1,8 +1,12 @@
-# # Spectral Dynamic Causal Modeling Tutorial
+# # Solving Inverse Problems with Spectral Dynamic Causal Modeling
 # # Introduction
 #
-# In this tutorial we will introduce how to perform a spectral Dynamic Causal Modeling analysis on simulated data [1,2].
-# To do so we roughly resemble the procedure in the [SPM12](https://www.fil.ion.ucl.ac.uk/spm/software/spm12/) script `DEM_demo_induced_fMRI.m` in [Neuroblox](https://www.neuroblox.org/).
+# Neuroblox provides you with a comprehensive environment for simulations as we have explored previously, but its functionality doesn't stop there.
+# We will now pivot and turn our attention to a different kind of problem: 
+# inferring model parameters, that is solving inverse problems, from time series. 
+# The method of choice is one of the most widely spread in imaging neuroscience, spectral Dynamic Causal Modeling (spDCM)[1,2]. 
+# In this tutorial we will introduce how to perform a spDCM analysis on simulated data.
+# To do so we roughly reproduce the procedure in the [SPM12](https://www.fil.ion.ucl.ac.uk/spm/software/spm12/) script `DEM_demo_induced_fMRI.m` in [Neuroblox](https://www.neuroblox.org/).
 # This work was also presented in Hofmann et al.[2]
 #
 # In this tutorial we will define a circuit of three linear neuronal mass models, all driven by an Ornstein-Uhlenbeck process.
@@ -56,6 +60,8 @@ end
 # Next we define the between-region connectivity matrix and make sure that it is diagonally dominant to guarantee numerical stability (see Gershgorin theorem).
 A_true = 0.1*randn(nr, nr)
 A_true -= diagm(map(a -> sum(abs, a), eachrow(A_true)))    # ensure diagonal dominance of matrix
+# Instead of a random matrix use the same matrix as is defined in [3]
+A_true = [[-0.5 -2 0]; [0.4 -0.5 -0.3]; [0 0.2 -0.5]]
 for idx in CartesianIndices(A_true)
     add_edge!(g, regions[idx[1]] => regions[idx[2]]; :weight => A_true[idx[1], idx[2]])
 end
@@ -127,7 +133,7 @@ for i = 1:nr
 end
 
 A_prior = 0.01*randn(nr, nr)
-A_prior -= diagm(diag(A_prior))    # ensure diagonal dominance of matrix
+A_prior -= diagm(diag(A_prior))    # remove the diagonal
 # Since we want to optimize these weights we turn them into symbolic parameters:
 # Add the symbolic weights to the edges and connect regions.
 A = []
@@ -137,7 +143,7 @@ for (i, a) in enumerate(vec(A_prior))
 end
 # With the function `untune!`` we can list indices of parameters whose tunable flag should be set to false.
 # For instance the first element in the second row:
-untune!(A, [4])
+untune!(A, [])
 for (i, idx) in enumerate(CartesianIndices(A_prior))
     if idx[1] == idx[2]
         add_edge!(g, regions[idx[1]] => regions[idx[2]]; :weight => -exp(A[i])/2)  # -exp(A[i])/2: treatement of diagonal elements in SPM12 to make diagonal dominance (see Gershgorin Theorem) more likely but it is not guaranteed
@@ -201,4 +207,5 @@ ecbarplot(state, setup, A_true)
 
 # ## References
 # [1] [Novelli, Leonardo, Karl Friston, and Adeel Razi. “Spectral Dynamic Causal Modeling: A Didactic Introduction and Its Relationship with Functional Connectivity.” Network Neuroscience 8, no. 1 (April 1, 2024): 178–202.](https://doi.org/10.1162/netn_a_00348) \
-# [2] [Hofmann, David, Anthony G. Chesebro, Chris Rackauckas, Lilianne R. Mujica-Parodi, Karl J. Friston, Alan Edelman, and Helmut H. Strey. “Leveraging Julia’s Automated Differentiation and Symbolic Computation to Increase Spectral DCM Flexibility and Speed.” bioRxiv: The Preprint Server for Biology, 2023.](https://doi.org/10.1101/2023.10.27.564407)
+# [2] [Hofmann, David, Anthony G. Chesebro, Chris Rackauckas, Lilianne R. Mujica-Parodi, Karl J. Friston, Alan Edelman, and Helmut H. Strey. “Leveraging Julia’s Automated Differentiation and Symbolic Computation to Increase Spectral DCM Flexibility and Speed.” bioRxiv: The Preprint Server for Biology, 2023.](https://doi.org/10.1101/2023.10.27.564407) \
+# [3] [Friston, Karl J., Joshua Kahan, Bharat Biswal, and Adeel Razi. “A DCM for Resting State fMRI.” NeuroImage 94 (July 2014): 396–407.](https://linkinghub.elsevier.com/retrieve/pii/S1053811913012135)
diff --git a/src/Neurographs.jl b/src/Neurographs.jl
@@ -25,7 +25,7 @@ function add_edge!(g::MetaDiGraph, p::Pair; kwargs...)
         add_blox!(g, dest)
         dest_idx = nv(g)
     end
-    
+
     add_edge!(g, src_idx, dest_idx, Dict(kwargs))
 end
 

diff --git a/src/blox/canonicalmicrocircuit.jl b/src/blox/canonicalmicrocircuit.jl
@@ -35,18 +35,17 @@ mutable struct CanonicalMicroCircuitBlox <: CompositeBlox
 
         g = MetaDiGraph()
         sblox_parts = vcat(ss, sp, ii, dp)
-        add_blox!.(Ref(g), sblox_parts)
-
-        add_edge!(g, 1, 1, :weight, -800.0)
-        add_edge!(g, 2, 1, :weight, -800.0)
-        add_edge!(g, 3, 1, :weight, -1600.0)
-        add_edge!(g, 1, 2, :weight,  800.0)
-        add_edge!(g, 2, 2, :weight, -800.0)
-        add_edge!(g, 1, 3, :weight,  800.0)
-        add_edge!(g, 3, 3, :weight, -800.0)
-        add_edge!(g, 4, 3, :weight,  400.0)
-        add_edge!(g, 3, 4, :weight, -400.0)
-        add_edge!(g, 4, 4, :weight, -200.0)
+
+        add_edge!(g, ss => ss; :weight => -800.0)
+        add_edge!(g, sp => ss; :weight => -800.0)
+        add_edge!(g, ii => ss; :weight => -1600.0)
+        add_edge!(g, ss => sp; :weight =>  800.0)
+        add_edge!(g, sp => sp; :weight => -800.0)
+        add_edge!(g, ss => ii; :weight =>  800.0)
+        add_edge!(g, ii => ii; :weight => -800.0)
+        add_edge!(g, dp => ii; :weight =>  400.0)
+        add_edge!(g, ii => dp; :weight => -400.0)
+        add_edge!(g, dp => dp; :weight => -200.0)
 
         # Construct a BloxConnector object from the graph
         # containing all connection equations from lower levels and this level.

diff --git a/src/blox/connections.jl b/src/blox/connections.jl
@@ -70,7 +70,7 @@ function generate_weight_param(blox_out, blox_in; kwargs...)
     if typeof(weight) == Num   # Symbol
         w = weight
     else
-        w = only(@parameters $(w_name)=weight)
+        w = only(@parameters $(w_name)=weight [tunable=false])
     end    
 
     return w
@@ -351,7 +351,7 @@ function (bc::BloxConnector)(
     push!(bc.weights, r)
 
     eq = sys_in.jcn ~ sigmoid(x, r)*w
-    
+
     accumulate_equation!(bc, eq)
 end
 
@@ -410,28 +410,20 @@ end
 # additional dispatch to connect to hemodynamic observer blox
 function (bc::BloxConnector)(
     bloxout::NeuralMassBlox, 
-    bloxin::ObserverBlox; 
-    weight=1,
-    delay=0,
-    density=0.1
-)
-    # Need t for the delay term
-    @variables t
+    bloxin::ObserverBlox;
+    kwargs...)
 
     sys_out = get_namespaced_sys(bloxout)
     sys_in = get_namespaced_sys(bloxin)
 
     if typeof(bloxout.output) == Num
-        w_name = Symbol("w_$(nameof(sys_out))_$(nameof(sys_in))")
-        if typeof(weight) == Num # Symbol
-            w = weight
-        else
-            w = only(@parameters $(w_name)=weight [tunable=false])
-        end    
+        w = generate_weight_param(bloxout, bloxin; kwargs...)
         push!(bc.weights, w)
         x = namespace_expr(bloxout.output, sys_out, nameof(sys_out))
         eq = sys_in.jcn ~ x*w
     else
+        # Need t for the delay term
+        @variables t
         # Define & accumulate delay parameter
         # Don't accumulate if zero
         τ_name = Symbol("τ_$(nameof(sys_out))_$(nameof(sys_in))")
@@ -453,18 +445,12 @@ end
 function (bc::BloxConnector)(
     bloxout::StimulusBlox,
     bloxin::NeuralMassBlox;
-    weight=1
-)
+    kwargs...)
 
     sys_out = get_namespaced_sys(bloxout)
     sys_in = get_namespaced_sys(bloxin)
 
-    w_name = Symbol("w_$(nameof(sys_out))_$(nameof(sys_in))")
-    if typeof(weight) == Num # Symbol
-        w = weight
-    else
-        w = only(@parameters $(w_name)=weight)
-    end    
+    w = generate_weight_param(bloxout, bloxin; kwargs...)
     push!(bc.weights, w)
 
     x = namespace_expr(bloxout.output, sys_out, nameof(sys_out))

diff --git a/test/Project.toml b/test/Project.toml
@@ -0,0 +1,20 @@
+[deps]
+CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
+DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
+DelayDiffEq = "bcd4f6db-9728-5f36-b5f7-82caef46ccdb"
+Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
+GraphDynamics = "bcd5d0fe-e6b7-4ef1-9848-780c183c7f4c"
+Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+MAT = "23992714-dd62-5051-b70f-ba57cb901cac"
+MetaGraphs = "626554b9-1ddb-594c-aa3c-2596fe9399a5"
+ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
+Neuroblox = "769b91e5-4c60-41ee-bfae-153c84203cb2"
+OrderedCollections = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
+OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
+SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
+StochasticDiffEq = "789caeaf-c7a9-5a7d-9973-96adeb23e2a0"
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
diff --git a/test/datafitting.jl b/test/datafitting.jl
@@ -13,7 +13,6 @@ using MAT
     dt = 2.0                                 # time bin in seconds
     freq = range(min(128, ns*dt)^-1, max(8, 2*dt)^-1, 32)  # define frequencies at which to evaluate the CSD
 
-
     ########## assemble the model ##########
     g = MetaDiGraph()
     regions = Dict()
@@ -23,13 +22,11 @@ using MAT
         add_blox!(g, region)
         regions[ii] = nv(g)    # store index of neural mass model
         taskinput = ExternalInput(;name=Symbol("r$(ii)₊ei"), I=1.0)
-        add_blox!(g, taskinput)
-        add_edge!(g, nv(g), regions[ii], Dict(:weight => C))
+        add_edge!(g, taskinput => region, weight = C)
         # add hemodynamic observer
         observer = BalloonModel(;name=Symbol("r$(ii)₊bm"), lnκ=lnκ, lnϵ=lnϵ)
-        add_blox!(g, observer)
         # connect observer with neuronal signal
-        add_edge!(g, regions[ii], nv(g), Dict(:weight => 1.0))
+        add_edge!(g, region => observer, weight = 1.0)
     end
 
     # add symbolic weights
@@ -130,14 +127,12 @@ end
         add_blox!(g, region)
         regions[ii] = nv(g)    # store index of neural mass model
         input = ExternalInput(;name=Symbol("r$(ii)₊ei"), I=1.0)
-        add_blox!(g, input)
-        add_edge!(g, nv(g), nv(g) - 1, Dict(:weight => C))
+        add_edge!(g, input => region; weight = C)
 
         # add lead field (LFP measurement)
         measurement = LeadField(;name=Symbol("r$(ii)₊lf"))
-        add_blox!(g, measurement)
         # connect measurement with neuronal signal
-        add_edge!(g, nv(g) - 2, nv(g), Dict(:weight => 1.0))
+        add_edge!(g, region => measurement; weight = 1.0)
     end
 
     nl = Int((nrr^2-nrr)/2)   # number of links unidirectional
@@ -191,15 +186,15 @@ end
     end
     indices = Dict(:dspars => collect(1:np))
     # Noise parameter mean
-    modelparam[:lnα] = zeros(Float64, 2, nrr);           # intrinsic fluctuations, ln(α) as in equation 2 of Friston et al. 2014 
+    modelparam[:lnα] = zeros(Float64, 2, nrr);        # intrinsic fluctuations, ln(α) as in equation 2 of Friston et al. 2014 
     n = length(modelparam[:lnα]);
     indices[:lnα] = collect(np+1:np+n);
     np += n;
-    modelparam[:lnβ] = [-16.0, -16.0];           # global observation noise, ln(β) as above
+    modelparam[:lnβ] = [-16.0, -16.0];                # global observation noise, ln(β) as above
     n = length(modelparam[:lnβ]);
     indices[:lnβ] = collect(np+1:np+n);
     np += n;
-    modelparam[:lnγ] = [-16.0, -16.0];   # region specific observation noise
+    modelparam[:lnγ] = [-16.0, -16.0];                # region specific observation noise
     indices[:lnγ] = collect(np+1:np+nrr);
     np += nrr
     indices[:u] = idx_u
-Original file line number
+Diff line change
@@ Expand Up / @@ -25,7 +25,7 @@ function add_edge!(g::MetaDiGraph, p::Pair; kwargs...) @@
             add_blox!(g, dest)
             dest_idx = nv(g)
         end
         add_edge!(g, src_idx, dest_idx, Dict(kwargs))
     end
@@ Expand Down @@