WIP: cleaning and revisiting spDCM tutorial

docs/Project.toml

@@ -8,13 +8,15 @@ Downloads = "f43a241f-c20a-4ad4-852c-f6b1247861c6"
 HypothesisTests = "09f84164-cd44-5f33-b23f-e6b0d136a0d5"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306"
+MAT = "23992714-dd62-5051-b70f-ba57cb901cac"
 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"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
+StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 StochasticDiffEq = "789caeaf-c7a9-5a7d-9973-96adeb23e2a0"
 
 [compat]

docs/src/tutorials/spectralDCM.jl

@@ -31,6 +31,7 @@ using OrderedCollections
 using CairoMakie
 using ModelingToolkit
 using Random
+using MAT
 
 # # Model simulation
 # ## Define the model
@@ -49,20 +50,28 @@ g = MetaDiGraph()
 regions = [];      # list of neural mass blocks to then connect them to each other with an adjacency matrix `A_true`
 # Now add the different blocks to each region and connect the blocks within each region. 
 # For convenience we use a for loop since the type of blocks belonging to a each region repeat over regions but you could also approach building the system the same way as was shown in previous tutorials:
+vars = matread(joinpath(@__DIR__, "../../../test/spm25_demo.mat"));
+vars["u"]
+tspan = (0.0, 1124.0)
+t_burnin = 102   # ignore the first 102 seconds
+dt = 2   # 2 seconds as measurement interval for fMRI
+times = t_burnin:dt:tspan[2]
 for i = 1:nr
     region = LinearNeuralMass(;name=Symbol("r$(i)₊lm"))
     push!(regions, region)          # store neural mass model in list. We need this list below. If you haven't seen the Julia command `push!` before [see here](http://jlhub.com/julia/manual/en/function/push-exclamation).
 
     ## add Ornstein-Uhlenbeck block as noisy input to the current region
-    input = OUBlox(;name=Symbol("r$(i)₊ou"), σ=0.1)
+    # input = OUBlox(;name=Symbol("r$(i)₊ou"), σ=0.25, τ=2.5)
+    input = ARBlox(;name=Symbol("r$(i)₊ar"), t_stim=vars["u"][:, i], timeseries=times)
+
     add_edge!(g, input => region, weight=1/16)   # Note that 1/16 is taken from SPM12, this stabilizes the balloon model simulation. Alternatively the noise of the Ornstein-Uhlenbeck block or the weight of the edge connecting neuronal activity and balloon model could be reduced to guarantee numerical stability.
 
     ## simulate fMRI signal with BalloonModel which includes the BOLD signal on top of the balloon model dynamics
     measurement = BalloonModel(;name=Symbol("r$(i)₊bm"))
     add_edge!(g, region => measurement, weight=1.0)
 end
 # Next we define the between-region connectivity matrix and connect regions; we 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]]
+A_true = [[-0.5 -0.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
@@ -73,9 +82,7 @@ end
 # ## Run the simulation and plot the results
 
 # setup simulation of the model, time in seconds
-tspan = (0.0, 512.0)
 prob = SDEProblem(simmodel, [], tspan)
-dt = 2   # 2 seconds (units are seconds) as measurement interval for fMRI
 sol = solve(prob, ImplicitRKMil(), saveat=dt);
 
 # we now want to extract all the variables in our model which carry the tag "measurement". For this purpose we can use the Neuroblox function `get_idx_tagged_vars`
@@ -93,10 +100,21 @@ lines!(ax, sol, idxs=idx_m)
 f
 
 # We note that the initial spike is not meaningful and a result of the equilibration of the stochastic process thus we remove it.
-dfsol = DataFrame(sol[ceil(Int, 101/dt):end]);
+dfsol = DataFrame(sol[round(Int, t_burnin/dt):end]);
 
 # ## Estimate and plot the cross-spectral densities
-data = Matrix(dfsol[:, idx_m]);
+# add rescaling of data as done in SPM:
+data = Matrix(dfsol[:, idx_m.+1]);
+data += vars["e"][1:size(data)[1], :]
+data .-= mean(data,dims=1)
+scale = 1/sqrt(sum((data[:]).^2)/length(data))/4;
+data *= scale
+ddata = DataFrame(data, :auto)
+
+# ddata = DataFrame(vars["data"], :auto)    # turn data into DataFrame, name column names after building the model.
+_, obsvars = get_eqidx_tagged_vars(simmodel, "measurement")  # get index of equation of bold state
+rename!(ddata, Symbol.(obsvars))
+
 # We compute the cross-spectral density by fitting a linear model of order `p` and then compute the csd analytically from the parameters of the multivariate autoregressive model
 p = 8
 mar = mar_ml(data, p)   # maximum likelihood estimation of the MAR coefficients and noise covariance matrix
@@ -143,6 +161,10 @@ end
 # Here we define the prior expectation values of the effective connectivity matrix we wish to infer:
 A_prior = 0.01*randn(nr, nr)
 A_prior -= diagm(diag(A_prior))    # remove the diagonal
+# These two parameters are not present in the ground truth thus set them to zero and set their tuning parameter to false: 
+A_prior[3] = 0.0
+A_prior[7] = 0.0
+
 # Since we want to optimize these weights we turn them into symbolic parameters:
 # Add the symbolic weights to the edges and connect regions.
 A = []
@@ -190,7 +212,7 @@ csdsetup = (mar_order = p, freq = freq, dt = dt);
 # and in particular we need to get the measurement variables in the same ordering as the model equations are defined.
 _, s_bold = get_eqidx_tagged_vars(fitmodel, "measurement");    # get bold signal variables
 # Prepare the DCM. This function will setup the computation of the Dynamic Causal Model. The last parameter specifies that wer are using fMRI time series as opposed to LFPs.
-(state, setup) = setup_sDCM(dfsol[:, String.(Symbol.(s_bold))], fitmodel, perturbedfp, csdsetup, priors, hyperpriors, indices, pmean, "fMRI");
+(state, setup) = setup_sDCM(ddata, fitmodel, perturbedfp, csdsetup, priors, hyperpriors, indices, pmean, "fMRI");
 
 # We are now ready to run the optimization procedure! :)
 # That is we loop over run_sDCM_iteration! which will alter `state` after each optimization iteration. It essentially computes the Variational Laplace estimation of expectation and variance of the tunable parameters. 

src/Neuroblox.jl

@@ -237,7 +237,7 @@ export LearningBlox
 export CosineSource, CosineBlox, NoisyCosineBlox, PhaseBlox, ImageStimulus, ConstantInput, ExternalInput, SpikeSource, PoissonSpikeTrain, generate_spike_times
 export DBS, ProtocolDBS, detect_transitions, compute_transition_times, compute_transition_values, get_protocol_duration
 export BandPassFilterBlox
-export OUBlox, OUCouplingBlox
+export OUBlox, OUCouplingBlox, ARBlox
 export phase_inter, phase_sin_blox, phase_cos_blox
 export SynapticConnections, create_rl_loop
 export add_blox!, get_system

src/blox/stochastic.jl

@@ -28,6 +28,26 @@ mutable struct OUBlox <: NeuralMassBlox
     end
 end
 
+mutable struct ARBlox <: StimulusBlox
+    namespace
+    system
+    function ARBlox(;name, namespace=nothing, t_stim, timeseries)
+        sts = @variables x(t)=0.0 [output=true, irreducible=true]
+
+        N = length(timeseries)
+        stim_eqs = Vector{Equation}(undef, N)
+        for i = 1:N
+            stim_eqs[i] = x ~ timeseries[i]
+        end
+        eq = x ~ 0.0
+        cb_stim = t_stim .=> stim_eqs
+        sys = ODESystem(eq, t, sts, []; name, discrete_events = cb_stim)
+
+        new(namespace, sys)
+    end
+end
+
+
 # """
 # Ornstein-Uhlenbeck Coupling Blox
 # This blox takes an input and multiplies that input with

src/datafitting/spDCM_VL.jl

@@ -104,9 +104,6 @@ function LinearAlgebra.eigen(M::Matrix{Dual{T, P, np}}) where {T, P, np}
 end
 
 function transferfunction(freq, derivatives, params, indices)
-    # nr = length(indices[:u])
-    # pars = params[indices[:dspars]]
-    # ∂f = derivatives([pars[1:nr^2], pars[nr^2+1:end]...])
     ∂f = derivatives(params[indices[:dspars]])
     ∂f∂x = ∂f[indices[:sts], indices[:sts]]
     ∂f∂u = ∂f[indices[:sts], indices[:u]]
@@ -119,7 +116,7 @@ function transferfunction(freq, derivatives, params, indices)
     ∂g∂v = ∂g∂x*V
     ∂v∂u = V\∂f∂u              # u is external variable which we don't use right now. With external variable this would read V/dfdu
 
-    nfreq = size(freq, 1)            # number of frequencies
+    nfreq = size(freq, 1)      # number of frequencies
     ng = size(∂g∂x, 1)         # number of outputs
     nu = size(∂v∂u, 2)         # number of inputs
     nk = size(V, 2)            # number of modes
@@ -314,15 +311,16 @@ function integration_step(dfdx, f, v, solenoid=false)
         dfdx = dfdx - Q*dfdx;        
     end
 
-    # (expm(dfdx*t) - I)*inv(dfdx)*f ~~~ could also be done with expv (expv(t, dFdθθ, dFdθθ \ dFdθ) - dFdθθ \ dFdθ) but doesn't work with Dual.
-    # could also be done with exponential! but isn't numerically stable
+    # NB: (exp(dfdx*t) - I)*inv(dfdx)*f, could also be done with expv (expv(t, dFdθθ, dFdθθ \ dFdθ) - dFdθθ \ dFdθ) but doesn't work with Dual.
+    # Could also be done with `exponential!` but isn't numerically stable.
+    # Thus, just use `exp`.
     n = length(f)
     t = exp(v - spm_logdet(dfdx)/n)
 
     if t > exp(16)
         dx = - dfdx \ f   # -inv(dfdx)*f
     else
-        dx = (exp(t * dfdx) - I) * inv(dfdx)*f # (expm(dfdx*t) - I)*inv(dfdx)*f
+        dx = (exp(t * dfdx) - I) * inv(dfdx) * f # (expm(dfdx*t) - I)*inv(dfdx)*f
     end
 
     return dx
@@ -462,15 +460,7 @@ function setup_sDCM(data, model, initcond, csdsetup, priors, hyperpriors, indice
     return (vlstate, vlsetup)
 end
 
-with_stack(f, n) = fetch(schedule(Task(f, n)));
-
 function run_sDCM_iteration!(state::VLState, setup::VLSetup)
-    with_stack(5_000_000) do
-        _run_sDCM_iteration!(state, setup)
-    end
-end
-
-function _run_sDCM_iteration!(state::VLState, setup::VLSetup)
     (;μθ_po, λ, v, ϵ_θ, dFdθ, dFdθθ) = state
 
     f = setup.model_at_x0

src/measurementmodels/fmri.jl

@@ -74,7 +74,7 @@ struct BalloonModel <: ObserverBlox
         p = paramscoping(lnκ=lnκ, lnτ=lnτ, lnϵ=lnϵ)     # finally compile all parameters
         lnκ, lnτ, lnϵ = p                               # assign the modified parameters
 
-        sts = @variables s(t)=1.0 lnu(t)=0.0 lnν(t)=0.0 lnq(t)=0.0 bold(t)=0.0 [irreducible=true, output=true, description="measurement"] jcn(t) [input=true]
+        sts = @variables s(t)=0.0 lnu(t)=0.0 lnν(t)=0.0 lnq(t)=0.0 bold(t)=0.0 [irreducible=true, output=true, description="measurement"] jcn(t) [input=true]
 
         eqs = [
             D(s)   ~ jcn - H[1]*exp(lnκ)*s - H[2]*(exp(lnu) - 1),
@@ -129,7 +129,7 @@ function boldsignal(;name, lnϵ=0.0)
     k1  = 4.3*nu0*E0*TE
 
     params = @parameters lnϵ=lnϵ
-    vars = @variables bold(t) lnν(t) lnq(t)   # TODO: got to be really careful with the current implementation! A simple permutation of this breaks the algorithm!
+    vars = @variables bold(t) lnν(t) lnq(t)
 
     eqs = [
         bold ~ V0*(k1 - k1*exp(lnq) + exp(lnϵ)*r0*E0*TE - exp(lnϵ)*r0*E0*TE*exp(lnq)/exp(lnν) + 1-exp(lnϵ) - (1-exp(lnϵ))*exp(lnν))