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