From ebb342179467e626050370e1cfae296faa659081 Mon Sep 17 00:00:00 2001 From: David Hofmann <1681922+david-hofmann@users.noreply.github.com> Date: Wed, 13 Nov 2024 20:13:45 +0100 Subject: [PATCH] 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 --- Project.toml | 6 +----- docs/src/tutorials/spectralDCM.jl | 19 ++++++++++++------ src/Neurographs.jl | 2 +- src/blox/canonicalmicrocircuit.jl | 23 +++++++++++----------- src/blox/connections.jl | 32 +++++++++---------------------- test/Project.toml | 20 +++++++++++++++++++ test/datafitting.jl | 19 +++++++----------- 7 files changed, 62 insertions(+), 59 deletions(-) create mode 100644 test/Project.toml diff --git a/Project.toml b/Project.toml index 50643b87..932522f5 100644 --- 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 index 4a361317..edfa48f8 100644 --- 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 index 55611546..807354a1 100644 --- 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 index ba16a104..dd72d05e 100644 --- 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 index b39916be..a3c7ad12 100644 --- 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 new file mode 100644 index 00000000..7e794922 --- /dev/null +++ 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 index 6518bf23..304a5294 100644 --- 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