Merge pull request #273 from Neuroblox/fmri_scope

implementing scope for fmri.jl
Neuroblox · Jul 18, 2023 · 9538c63 · 9538c63
2 parents d7de86a + fd98cf0
commit 9538c63
Show file tree

Hide file tree

Showing 3 changed files with 22 additions and 55 deletions.
diff --git a/src/blox/neural_mass.jl b/src/blox/neural_mass.jl
@@ -2,11 +2,9 @@
 D = Differential(t)
 
 mutable struct LinearNeuralMassBlox <: NBComponent
-    # τ::Num
     connector::Num
     odesystem::ODESystem
     function LinearNeuralMassBlox(;name)
-        # params = @parameters τ=τ
         states = @variables x(t) jcn(t)
         eqs = D(x) ~ jcn
         odesys = ODESystem(eqs, t, states, []; name=name)
@@ -16,22 +14,23 @@ end
 
 mutable struct HarmonicOscillatorBlox <: NeuralMassBlox
     # all parameters are Num as to allow symbolic expressions
-    ω::Num
-    ζ::Num
-    k::Num
-    h::Num
+    p_dict::Dict{Symbol,Union{Real,Num}}
     connector::Num
     noDetail::Vector{Num}
     detail::Vector{Num}
     initial::Dict{Num, Tuple{Float64, Float64}}
     odesystem::ODESystem
     function HarmonicOscillatorBlox(;name, ω=25*(2*pi), ζ=1.0, k=625*(2*pi), h=35.0)
-        params = @parameters ω=ω ζ=ζ k=k h=h
+        para_dict = scope_dict(Dict{Symbol,Union{Real,Num}}(:ω => ω,:ζ => ζ,:k => k,:h => h))
+        ω=para_dict[:ω]
+        ζ=para_dict[:ζ]
+        k=para_dict[:k]
+        h=para_dict[:h]
         sts    = @variables x(t)=1.0 y(t)=1.0 jcn(t)=0.0
         eqs    = [D(x) ~ y-(2*ω*ζ*x)+ k*(2/π)*(atan((jcn)/h))
                   D(y) ~ -(ω^2)*x]
-        odesys = ODESystem(eqs, t, sts, params; name=name)
-        new(ω, ζ, k, h, odesys.x,[odesys.x],[odesys.x,odesys.y],
+        odesys = ODESystem(eqs, t, sts, values(para_dict); name=name)
+        new(para_dict, odesys.x,[odesys.x],[odesys.x,odesys.y],
             Dict(odesys.x => (-1.0,1.0), odesys.y => (-1.0,1.0)),
             odesys)
     end
@@ -44,10 +43,7 @@ const harmonic_oscillator = HarmonicOscillatorBlox
 #    return HarmonicOscillatorImage
 
 mutable struct JansenRitCBlox <: NeuralMassBlox
-    τ::Num
-    H::Num
-    λ::Num
-    r::Num
+    p_dict::Dict{Symbol,Union{Real,Num}}
     connector::Num
     noDetail::Vector{Num}
     detail::Vector{Num}
@@ -63,7 +59,7 @@ mutable struct JansenRitCBlox <: NeuralMassBlox
         eqs    = [D(x) ~ y - ((2/τ)*x),
                 D(y) ~ -x/(τ*τ) + (H/τ)*((2*λ)/(1 + exp(-r*(jcn))) - λ)]
         odesys = ODESystem(eqs, t, sts, values(para_dict); name=name)
-        new(τ, H, λ, r, odesys.x,[odesys.x],[odesys.x,odesys.y],
+        new(para_dict, odesys.x,[odesys.x],[odesys.x,odesys.y],
             Dict(odesys.x => (-1.0,1.0), odesys.y => (-1.0,1.0)),
             odesys)
     end
@@ -72,10 +68,7 @@ end
 const jansen_ritC = JansenRitCBlox
 
 mutable struct  JansenRitSCBlox <: NeuralMassBlox
-    τ::Num
-    H::Num
-    λ::Num
-    r::Num
+    p_dict::Dict{Symbol,Union{Real,Num}}
     connector::Num
     noDetail::Vector{Num}
     detail::Vector{Num}
@@ -91,7 +84,7 @@ mutable struct  JansenRitSCBlox <: NeuralMassBlox
         eqs    = [D(x) ~ y - ((2/τ)*x),
                   D(y) ~ -x/(τ*τ) + (H/τ)*((2*λ)/(1 + exp(-r*(jcn))) - λ)]
         odesys = ODESystem(eqs, t, sts, values(para_dict); name=name)
-        new(τ, H, λ, r, odesys.x,[odesys.x],[odesys.x,odesys.y],
+        new(para_dict, odesys.x,[odesys.x],[odesys.x,odesys.y],
             Dict(odesys.x => (-1.0,1.0), odesys.y => (-1.0,1.0)),
             odesys)
     end

diff --git a/src/measurementmodels/fmri.jl b/src/measurementmodels/fmri.jl
@@ -24,8 +24,7 @@ lnτ  : logarithmic prefactor to transit time H[3], set to 0 for standard parame
 returns an ODESystem of the biophysical model for the hemodynamics
 """
 mutable struct Hemodynamics <: NBComponent
-    lnκ::Num
-    lnτ::Num
+    p_dict::Dict{Symbol, Union{Real, Num}}
     connector::Num
     odesystem::ODESystem
     function Hemodynamics(;name, lnκ=0.0, lnτ=0.0)
@@ -37,8 +36,9 @@ mutable struct Hemodynamics <: NBComponent
             H(5) - resting state oxygen extraction                (E0)
         =#
         H = [0.64, 0.32, 2.00, 0.32, 0.4]
-
-        params = @parameters lnκ=lnκ lnτ=lnτ
+        para_dict = scope_dict(Dict{Symbol, Union{Real,Num}}(:lnκ => lnκ, :lnτ => lnτ))
+        lnκ=para_dict[:lnκ]
+        lnτ=para_dict[:lnτ]
         states = @variables s(t) lnf(t) lnν(t) lnq(t) jcn(t)
 
         eqs = [
@@ -47,15 +47,13 @@ mutable struct Hemodynamics <: NBComponent
             D(lnν) ~ (exp(lnf) - exp(lnν)^(H[4]^-1)) / (H[3]*exp(lnτ)*exp(lnν)),
             D(lnq) ~ (exp(lnf)/exp(lnq)*((1 - (1 - H[5])^(exp(lnf)^-1))/H[5]) - exp(lnν)^(H[4]^-1 - 1))/(H[3]*exp(lnτ))
         ]
-        odesys = ODESystem(eqs, t, states, params; name=name)
-        new(lnκ, lnτ, Num(0), odesys)
+        odesys = ODESystem(eqs, t, states, values(para_dict); name=name)
+        new(para_dict, Num(0), odesys)
     end
 end
 
 
 mutable struct LinHemo <: NBComponent
-    lnκ::Num
-    lnτ::Num
     connector::Num
     bloxinput::Num
     odesystem::ODESystem
@@ -69,33 +67,10 @@ mutable struct LinHemo <: NBComponent
         add_vertex!(g, :blox, hemo)
         add_edge!(g, 1, 2, :weight, 1.0)
         linhemo = ODEfromGraph(g; name=name)
-        new(lnκ, lnτ, linhemo.nmm₊x, linhemo.jcn, linhemo)
+        new(linhemo.nmm₊x, linhemo.jcn, linhemo)
     end
 end
 
-# function hemodynamics(;name, lnκ=0.0, lnτ=0.0)
-#     #= hemodynamic parameters
-#         H(1) - signal decay                                   d(ds/dt)/ds)
-#         H(2) - autoregulation                                 d(ds/dt)/df)
-#         H(3) - transit time                                   (t0)
-#         H(4) - exponent for Fout(v)                           (alpha)
-#         H(5) - resting state oxygen extraction                (E0)
-#     =#
-#     H = [0.64, 0.32, 2.00, 0.32, 0.4]
-
-#     params = @parameters lnκ=lnκ lnτ=lnτ
-#     states = @variables s(t) lnf(t) lnν(t) lnq(t) jcn(t)
-
-#     eqs = [
-#         D(s)   ~ jcn - H[1]*exp(lnκ)*s - H[2]*(exp(lnf) - 1),
-#         D(lnf) ~ s / exp(lnf),
-#         D(lnν) ~ (exp(lnf) - exp(lnν)^(H[4]^-1)) / (H[3]*exp(lnτ)*exp(lnν)),
-#         D(lnq) ~ (exp(lnf)/exp(lnq)*((1 - (1 - H[5])^(exp(lnf)^-1))/H[5]) - exp(lnν)^(H[4]^-1 - 1))/(H[3]*exp(lnτ))
-#     ]
-
-#     return ODESystem(eqs, t, states, params; name=name)
-# end
-
 """
 BOLD signal model as described in: 
 

diff --git a/test/datafitting.jl b/test/datafitting.jl
@@ -9,11 +9,10 @@ nd = ncol(data)                         # number of dimensions
 
 ########## assemble the model ##########
 
-# TODO: figure out a good way to deal with shared parameters over regions
-# @parameters lnκ=0.0     # define brain-wide decay parameter for hemodynamics
+@parameters lnκ=0.0     # define brain-wide decay parameter for hemodynamics
 g = MetaDiGraph()
 for ii = 1:nd
-    region = LinHemo(;name=Symbol("r$ii"))
+    region = LinHemo(;name=Symbol("r$ii"), lnκ=lnκ)
     add_blox!(g, region)
 end
 
@@ -83,4 +82,4 @@ results = spectralVI(data, neuronmodel, bold, initcond, csdsetup, params, hyperp
 ### COMPARE RESULTS WITH MATLAB RESULTS ###
 @show results.F, vars["F"]
 @test results.F < vars["F"]*0.99
-@test_broken results.F > vars["F"]*1.01
+@test results.F > vars["F"]*1.01