Inferring effective connectivity with rDCM

In this example we show how to invert a linear DCM with 50 regions using both rDCM with a fixed network architecture and sparse rDCM. In addition to the RegressionDynamicCausalModeling.jl package, Plots.jl and Random.jl are required.

rDCM with fixed network architecture

using RegressionDynamicCausalModeling
using Plots
using Random

dcm = load_example_DCM() # this loads the example DCM

Linear DCM
a:     50x50 matrix
c:     50x25 matrix
scans: 2714
nr:    50
---------------------------------------------
U (Input)
   u:  43424x25 matrix
   dt: 0.03125s
   names: u_1,...,u_25
---------------------------------------------
Y (BOLD signal)
   y:  2714x50 matrix
   dt: 0.5s
   names: y_1,...,y_50
---------------------------------------------
Ep (True parameters)
    A: 50 x 50 matrix
    C: 50 x 25 matrix
    transit: -0.06 ... -0.05
    decay:   0.10 ... 0.10
    epsilon: -0.05

We can visualize the network architecture of this DCM:

heatmap(
    dcm.a;
    yflip=true,
    title="Network architecture (A matrix)",
    titlefontsize=8,
    legend=:none,
    xlabel="region from",
    ylabel="region to",
    aspect_ratio=1,
    size=(350, 350),
)

We can use the DCM to generate synthetic data

y_noise, _, _, _ = generate_BOLD(dcm; SNR=3, rng=MersenneTwister(42));
dcm.Y.y = y_noise # save BOLD signal in DCM
p1 = plot(dcm.Y.y[:, 1]; label="true signal", title="BOLD signal of region 1")

Specify options/hyperparameters for inversion

opt = Options(RigidInversionParams(); synthetic=true, rng=MersenneTwister(42))

rdcm = RigidRdcm(dcm) # convert the DCM to a rDCM model

Estimate connectivity

output = invert(rdcm, opt)

rigid rDCM output
    F:   -574988.69
    F_r: -14474.90 ... -8427.91
    iterations until convergence per region: 4 ... 5
    Posteriors:
        α: 1359.00,...,1359.00
        β: 3308875.41,...,36960.51
        μ: 50 x 75 matrix
        Σ: 50 element vector of matrices

Based on the estimated parameters we can predict the BOLD signal and compare it to the ground truth

y_pred = predict(rdcm, output)
plot(p1, y_pred[:, 1]; label="predicted signal")

Visualize the posterior mean of the A matrix

heatmap(
    output.m_all[:, 1:50];
    yflip=true,
    title="Posterior mean A matrix",
    titlefontsize=8,
    xlabel="region from",
    ylabel="region to",
    aspect_ratio=1,
    size=(350, 350),
)

sparse rDCM

We load again the example DCM and generate BOLD signal with an SNR of 3.

dcm = load_example_DCM()
y_noise, _, _, _ = generate_BOLD(dcm; SNR=3, rng=MersenneTwister(42));
dcm.Y.y = y_noise;

Convert the linear DCM to a sparse rDCM model with a Bernoulli prior of 0.15 (a priori belief about network sparsity). This means we assume that the a priori probability of a connection being present is low.

rdcm = SparseRdcm(dcm; p0=0.15)

Specify options/hyperparameters for inversion

opt = Options(
    SparseInversionParams(; reruns=100, restrictInputs=true);
    synthetic=true,
    rng=MersenneTwister(42),
)

Invert the model

output = invert(rdcm, opt)

sparse rDCM output
    F:   -572101.37
    F_r: -14436.11 ... -8377.39
    iterations until convergence per region: 9 ... 8
    Posteriors:
        α: 1359.00,...,1359.00
        β: 3167112.30,...,34529.65
        μ: 50 x 75 matrix
        Σ: 50 element vector of matrices
        Z: 50 x 75 matrix

We can inspect the Bernoulli posterior for each connection

heatmap(
    output.z_all[:, 1:50];
    yflip=true,
    title="Posterior over binary indicator variables",
    titlefontsize=8,
    xlabel="region from",
    ylabel="region to",
    aspect_ratio=1,
    size=(350, 350),
)

This page was generated using Literate.jl.