updated CMC code

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

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)

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
 

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.

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

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"

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