rDCM: simulation and HRF creation

[zenkavi]

Hi,

Thank you again for all your work and sharing it with us. I had some specific questions regarding simulation and HRF creation in your toolbox:

• In tapas_rdcm_generate.m​ (L62) why do you append the input three times before generating the neural states?

• Why do you then not use the first output of tapas_dcm_euler_gen.m​ for the BOLD signal but convolve it later in the frequency domain only to take (/downsample to) the middle chunk?

• What does DCM.U.dt represent in your DCM structures? I'm asking this to understand how it works in generating the HRF. In the DCM loaded into the tapas_rdcm_tutorial.m​ this is specified as 1/32. But e.g. in DEM_demo_induced_fmri.m​ that you used in your HBM paper it is set to the TR, which is 2. When I simulate data using DEM_demo_induced_fmri.m​ it breaks tapas_rdcm_estimate.m​ because the HRF isn't computed correctly. How did you get around this for your simulations?

Thanks!

[StefanFraessle]

Hi Zeynap,

Thank you very much for your kind words and also for your questions regarding the rDCM toolbox.

Your first two questions can be answered in one go. The reason why the input appears three times is because we do not want to induce any boundary (edge) effects at the beginning and the end of the time series when doing the FFT and, in particular, the convolution. Therefore, we essentially add the input at the beginning and at the end of the actually input we are interested in. This is also the reason why we take "the middle chunk" when finally getting the predicted BOLD signal time series y. I'm sure there are other ways to do this, but it is a relatively simple and straightforward solution.

Regarding the third question: In principle, DCM.U.dt refers to the microtime resolution of the inputs. This - to some degree - is a relict from the classical variants in DCM for fMRI implemented in SPM. More specifically, in SPM, the inputs are sampled at a higher frequency than the fMRI data - namely, at 16x the temporal resolution of the data. In other words, if your TR (i.e., the sampling rate of your fMRI scans) was 1 second, DCM.U.dt = 1/16. The reason for this is that fMRI data is usually acquired at relatively slow rates and you can obtain better predictors (i.e., regressors in GLMs, or inputs in DCMs) if you sample the predictors at a higher rate. This is also why DCM.U.dt = 1/32 in tapas_rdcm_tutorial.m as this resembled a case of synthetic data sampled at a fast TR = 0.5s.

From this explanation, you might already see that the higher sampling is not really essential to the process of generating data but simply serves to make the predictors "more accurate". Why exactly the inputs are then not sampled at a higher rate in the DEM_demo_induced_fmri.m script, I cannot tell you for sure. My intuition would be that, in this particular case, it is not really necessary because the inputs (i.e., the AR(1)-noise process) do not have any fast transients (in contrast to, e.g., typical task-like block regressors).

Now, in terms of model inversion with tapas_rdcm_estimate.m: My suspicion is that there is something wrong with how you set up the model (in particular the DCM.U). In principle, rDCM switches off the inputs when modeling the resting state. However, if you simply input the DCM structures from the DEM_demo_induced_fmri.m that you have used to generate the synthetic data (i.e., with the AR(1) processes as driving inputs), rDCM will assume that this is a task-based model and therefore throw an error.

In other words, what you need to do is to specify a new DCM structure with the generated synthetic data but without the AR(1) inputs. You can do this by simply removing the field U from the DCM structure that you have used to simulate the resting-state fMRI data: DCM = rmfield(DCM,'U'). The function tapas_rdcm_estimate.m will then automatically specify the resting-state model correctly and the error should not occur anymore.

I hope this is clear. Especially the resting-state part is a bit confusing simply due to the fact that you are essentially generating synthetic fMRI data from a different model (classical DCM for fMRI using an AR(1) process as input) than you are using to invert the data (rDCM).

Let me know if you have any further questions.

All the best, Stefan

[zenkavi]

Thank you so much for the detailed response! It is incredibly helpful.

Following your suggestion to remove DCM.U before feeding the structure into tapas_rdcm_estimate.m has almost worked perfectly except for two issues. I think I've debugged the first but need to work more on the second:

• tapas_rdcm_estimate.m specifies DCM.U.u as a column vector for endogenous DCM's. This seems to cause a problem for tapas_dcm_euler_gen.m that is eventually called when finalizing the output because the column vector of DCM.U.u doesn't work with Ep.C'. Changing

  DCM.U.u     = zeros(size(DCM.Y.y,1)*16, 1);

  to

  DCM.U.u     = zeros(size(DCM.Y.y,1)*16, size(DCM.Y.y,2));

  does fix this but it might not be the best solution for large networks though so I haven't sent a PR. Does this make sense or might I have misspecified something else in my DCM?

• I'm also getting some MEX issues and matlab crashing when I try tapas_rdcm_compute_statistics.m but I haven't debugged this yet.

[zenkavi]

Also in tapas_rdcm_create_regressors.m is it supposed to be an empty matrix

DCM.U.X0      = zeros(size(DCM.U.u,1),0);

or a column of zeros?

DCM.U.X0      = zeros(size(DCM.U.u,1),1);

[StefanFraessle]

Hi Zeynap,

My apologies - I forgot to include this in my previous response: When simply taking the DCMs from the simulation step into tapas_rdcm_estimate.m, you need to first also remove the fields Ep and M from the DCM structure using DCM = rmfield(DCM,'Ep'); and DCM = rmfield(DCM,'M'); These are fields that are written when simulating the synthetic data but which will mess up your model inversion since they should not be specified before model inversion.

Alternatively, you could set up a fresh DCM structure for model inversion to prevent these issues. You could use the function tapas_rdcm_model_specification.m to do that. In your case, this could be achieved simply by calling the function like this: DCM = tapas_rdcm_model_specification(Y,[],args) where

• Y is the full Y-structure (i.e., containing Y.y, Y.dt, and Y.name) that you obtain from the simulation of the synthetic data
• args has a single field args.a which encodes the structure of your A-matrix (i.e., essentially DCM.a of the model from which you have simulated data).

Does that fix the issue? If not, I would be wondering which rDCM version you are using. It might well be that you are still using an older version of the rDCM toolbox which was not yet designed to function with resting-state fMRI data. The most recent version of rDCM should handle the specification of the DCM structure correctly for rs-fMRI data so that no errors are produced and no fixes from your side are necessary.

If this does not solve the issue, please feel free to send me your DCM file that you input into tapas_rdcm_estimate.m - it's so much easier to debug problems when you have the files available to interrogate :)

And, for resting-state data, DCM.U.X0 = zeros(size(DCM.U.u,1),0); is correct.

Regarding the MEX issue and Matlab crashing, this is a bit difficult to judge for me from your description. I would probably need more details on that. Could this be a memory issue?

Hope this helps.

With very best wishes, Stefan

[zenkavi]

Hi Stefan,

Thank you again for the detailed response.

My DCM structure didn't have M or Ep fields so this wasn't the issue. Trying to use DCM = tapas_rdcm_model_specification(Y,[],args) also didn't work because tapas_rdcm_generate.m needs DCM.Tp which for type = 's' is called before tapas_rdcm_compute_signals.m sets DCM.Tp = output.Ep.

Even when I specified DCM.Tp for the output of tapas_rdcm_model_specification.m things still crashed due to the MEX issue (explained below what triggers it).

The only way I've managed to run tapas_rdcm_estimate.m without errors (though that doesn't mean the correct inversion is being made since I get a very different `

Ep.AthanTp.A`) is by making three changes:

1. Comment out lines 42-45 in tapas_rdcm_compute_signals.m. If I don't then DCMs = tapas_rdcm_generate(DCM, options, Inf); in line 41 creates MEX crash with the following error.

Abnormal termination:
Segmentation violation

2. Changing the dimensions of DCM.U.u in tapas_rdcm_estimate.m line 74 to:

DCM.U.u     = zeros(size(DCM.Y.y,1)*16, size(DCM.Y.y,2));

If I don't make this change I get dimension error triggered by line 54 of tapas_dcm_euler_gen.m.

3. Changing DCM.U.X0 = zeros(size(DCM.U.u,1),0); to DCM.U.X0 = zeros(size(DCM.U.u,1),1); on line 77 of tapas_rdcm_create_regressors.m. If I don't make this change then I get the following error:

Unable to perform assignment because the size of the left side is 1-by-1 and the size of the right side is 0-by-1.

Error in tapas_rdcm_compute_signals (line 91)
        DCM_rDCM.Y.y(t,r) = DCM_rDCM.Y.y(t,r) + DCM_rDCM.Tp.baseline(r,:) * DCM_rDCM.U.X0((1+r_dt*(t-1)),:)';

My version in tapas_rdcm_estimate.m is set as output.ver = '2020_v01.2';. Should I be using another version?

I'll email the DCM structure I've been working with.

Thank you again so much for your help.

Best,

Zeynep

[StefanFraessle]

Hi Zeynep,

Thank you for your detailed description and for sending the DCM structure - this has certainly helped to understand the problem.

There is a couple of things to note. First, you indeed stumbled across a small bug in the code :) Thank you very much for digging so deeply into this. There seems to be a small issue for the resting-state implementation of rDCM when choosing type = s. I simply had not implemented this properly and the script still tries to do all the things that it would do for task-based fMRI. This is causing trouble.

I checked my scripts from the HBM paper and I had used type = r because, for resting-state data, the two options are identical (I will explain in a second). This is why I had not noticed the bug before. In other words, if I run tapas_rdcm_estimate.m with the DCM structure you sent me and type = r, everything runs smoothly. This is what I would suggest you do for now - no changes to the code are necessary then. I will provide a fix for this bug as soon as I have a bit more time.

For completeness, let me explain, what type = s and type = r does. There are some differences in terms of outputting the results. But the main difference for model inversion is whether you estimate the "baseline" parameter (i.e., the parameter that shifts the predicted BOLD signal time series up or down) or not. For empirical data, you want to estimate this parameter since fMRI data can be shifted away from zero. However, for simulated data, rDCM fixes the baseline parameter to zero because this is how we generate the synthetic data. This makes life a bit easier for rDCM in simulations. This is also why the two options are identical for rs-fMRI data: There is no baseline for the resting-state model as we set DCM.U.X0 = zeros(size(DCM.U.u,1),0);. So, you can happily proceed with your project by simply using type = r.

Now, looking at your DCM structure, there is one thing that you should definitely change. Your A-matrix does not contain any diagonal entries which means that you do not have any inhibitory self-connections in your model. This is not correct and will most likely lead to incorrect results. Hence, I suggest that you re-run the generation of the synthetic data as well as the inversion of the DCMs when enabling all of the inhibitory self-connections (i.e., putting "1"s on the diagonal of your DCM.a).

Regrading rDCM version: From what you describe, you should have the most recent version - but if you want to be 100% sure, you can grab the latest version of TAPAS from here.

I hope this helps. Let me know if you still run into issues.

All the best Stefan

[zenkavi]

Hi Stefan,

Thank you again so much for taking the time to help me out! Setting type = 'r' does indeed help things run smoothly.

Regarding inhibitory self-connections and more generally how to set DCM.a: Would not setting the diagonals of DCM.a as 1 lead to misleading results because the data generating process (AR(1)) does have an important self-connection effect? Do you have any intuition as to why they are left to be 0 in the default DEM_demo_induced_fMRI.m? Have you changed this in your analyses for your HBM paper?

I'm also a bit confused about the role of DCM.a especially when pruning the connectivity estimates. I need to re-read your 2018 paper and go over tapas_rdcm_sparse.m to understand the details of the pruning. From a first pass it seems like there is a strong assumption that each node can have only one connectivity parameter based on L66 of tapas_rdcm_sparse.m

idx_x = [ones(1,size(DCM.a,1)), zeros(1,size(DCM.b,1)*size(DCM.b,3)), ones(1,size(DCM.c,2))] == 1;

Is that correct? Then does the values of DCM.a even matter when using this for estimation?

Best,

Zeynep

[StefanFraessle]

Hi Zeynep,

Glad everything is running smoothly now.

The inhibitory self-connections always have to be present and are enforced in any variant of DCM. If these connections were not present, the dynamical system (i.e., the DCM) would show "runaway activity" where the system would never settle back to baseline after an input is provided and would tend towards infinity. This is not a realistic model of the brain.

Looking at the DEM_demo_induced_fMRI.m, you are correct that in DCM.a (l. 116 in the function), the self-connections are absent. However, I strongly suspect this is a little bug in the script that does not have any consequences (see below). This is because for the true parameters, the self-connections in the A matrix are present (line 49 and lines 58-60 in the function). Note that the values on the diagonal in pP.A are zero. However, keep in mind that there is a reparametrization of the diagonal parameters in SPM which ensures that the self-connections always stay negative. More precisely, if pP.A(1,1) = 0, this means that the actual self-connection is -0.5*exp(pP.A(1,1)) = -0.5. Hence, my suspicion is that when Adeel defines DCM.a = logical(pP.A) in line 116, this is just a small bug.

And it will not have an impact on model inversion because SPM automatically enforces diagonal elements - in other words, if the diagonal of DCM.a is not present, SPM will set those within spm_dcm_estimate.m (or some other function, I'm not 100% sure). But you can check that the self-connections are actually present after model inversion if you check the entries in CSD.Ep.A (which are the posterior means). You see that here the diagonal is not 0 anymore :)

Unfortunately, for rDCM, I have not included code yet that enforces the self-connections to be present (probably something I should do in a future release). Hence, in rDCM, having zeros on the diagonal will get you into trouble.

Regarding rDCM with sparsity constraints: No, this is not quite correct. The line that you mention above says that all connections to the current target region as well as all the driving inputs to the current target region are present. More specifically, for each target region, idx_x is a vector with R ones (where R is the number of regions), then zeros for the non-existing modulatory influences (as rDCM is only a linear model), and then K ones again (where K is the number of driving inputs). This is the starting scenario of sparse rDCM. The model then starts to prune this full network, based on the data that were acquired.

It's certainly a good idea to check out the 2018 NeuroImage paper if you're interested in more details about rDCM with sparsity constraints.

I hope this makes sense and was helpful.

All the best, Stefan

[zenkavi]

Hi Stefan,

Yep, makes sense.

Once again thank you very much for your detailed responses. I really appreciate them. It's been incredibly helpful since I'm pretty new to all this.

I'd be happy to help with extensions/bug fixes soon or later in the future if it lightens your load. Feel free to ping me anytime.

Best,

Zeynep