Bug in DCM-ERP tutorial: can't invert due to missing BEM data

[pranaysy]

Overview of problem

After specifying the model, as outlined in the first four steps of the DCM for evoked responses tutorial, inversion immediately fails with:

SPM26: spm_eeg_lgainmat                    14:06:10 - 28/04/2026
        ----------------------------------------------------------------
Error using ft_read_headmodel (line 46)
file 'C:\Users\at16431\OneDrive - University of Bristol\Ashley_AD_EEG\Anna_AI_taskcodes\Berk_new\New_spm\SPM\SPM\canonical\single_subj_T1_EEG_BEM.mat' does not exist

What's going on

The dataset linked in the docs for DCM for evoked response, points to dcm-erp-tyrer.zip with the following hash:

sha256sum: 56f364b9db49f3a9d1623754fce063424227095b6eaf236d5069615de27d659b

The contents of this archive are a single subject's processed data in SPM format:

mfaeffspmeeg_examplecontrolsubject.mat
mfaeffspmeeg_examplecontrolsubject.dat

Loading these data using SPM and inspecting their existing inversion fields reveals two things: there are five inversions in total with two unique gain matrices, and crucially these gain matrices are not supplied:

D = spm_eeg_load('mfaeffspmeeg_examplecontrolsubject');
numel(D.inv)
% ans =
%      5

% Loop over all inversions
cellfun(@(x)(fprintf('Gain matrix file: %s with comment: %s\n', x.gainmat, x.comment{:})), D.inv, 'UniformOutput', false)
  Gain matrix file: SPMgainmatrix_aeffspmeeg_8028_priming3_proc_1.mat with comment: 
  Gain matrix file: SPMgainmatrix_aeffspmeeg_8028_priming3_proc_1.mat with comment: ERP_CMM_NMDA_DCM_SENSOR_PRIMINGpmeeg_8028_priming3_proc.mat
  Gain matrix file: SPMgainmatrix_mfaeffspmeeg_8028_priming3_proc_3.mat with comment: DCM_ERP_CMM_NMDA_DCM_SENSOR_PRIMINGpmeeg_8028_priming3_proc.mat
  Gain matrix file: SPMgainmatrix_mfaeffspmeeg_8028_priming3_proc_3.mat with comment: DCM_30-Apr-2021
  Gain matrix file: SPMgainmatrix_mfaeffspmeeg_8028_priming3_proc_3.mat with comment: DCM_30-Apr-2021

Due to the missing leadfields, SPM tries to compute one from the forward model. However that fails because the paths to the BEM surfaces are hardcoded in each inversion:

unique(cellfun(@(x)(x.forward.vol), D.inv, 'UniformOutput', false)')
ans =
  2×1 cell array
    {'/users/at16431/Documents/MATLAB/SPM/spm12/spm12/canonical/single_subj_T1_EEG_BEM.mat'                                                           }
    {'C:\Users\at16431\OneDrive - University of Bristol\Ashley_AD_EEG\Anna_AI_taskcodes\Berk_new\New_spm\SPM\SPM\canonical\single_subj_T1_EEG_BEM.mat'}

The inversions appears to be using BEM surfaces derived from canonical meshes provided with SPM, however the paths to those files are local to the system used for respective inversions. And these paths don't automatically update over to a user's SPM installation unfortuntately.

Solution

Three possible paths forward:

1. Supply the gain matrix along with processed data in the archive for maximum reproducibility
2. Use same data and inversions, but provide a helper function to update paths for meshes etc on users' systems to their SPM installations
3. Use same data, discard all existing inversions and create a forward model from scratch using SPM defaults

I vote for 1 > 3 > 2, as we've done the same in a reliable manner previously: https://github.com/pranaysy/MultimodalDCM

[vlitvak]

Good catch! I think 1 might not work unless you manually change the path in inv to be relative. And that might fail as well if the dataset is not in the current folder. So I would vote for 3. We just need to add recreating the head model to the turorial either using the GUI or batch.

[pranaysy]

Thanks @vlitvak, yes I think 3 would be simpler actually. That's what I used to get it to work in the end, using this snippet:

% Load data
D = spm_eeg_load('mfaeffspmeeg_examplecontrolsubject.mat');

% Remove all existing inversions
if isfield(D, 'inv'); D = rmfield(D, 'inv'); end;

% Save and reload
D.save()
D.reload()

% Make a new forward model using template mesh and 'normal' resolution
val = 1;
D.val = val;
D = spm_eeg_inv_mesh_ui(D, val, true, 2);
D = spm_eeg_inv_datareg_ui(D, val);         % Default coreg

% Make BEM surfaces for EEG
D.inv{val}.forward = struct([]);
D.inv{val}.forward(1).voltype = 'EEG BEM';
D = spm_eeg_inv_forward(D);
D.save()

% Check forward solution to confirm all ok
spm_eeg_inv_checkforward(D, val);

(In case it helps others, or perhaps to incorporate in the tutorial)

---

I further note four things:

1. With the above snippet, I am able to reproduce the tutorial results While log-evidence is in favor of model 2 over model 1, just like the tutorial - there is a large numerical difference in the relative F values: ~200 in the tutorial docs, and ~150 in mine. Probably does not matter, and likely due to SPM/MATLAB version differences I suspect.

2. The sensor-level data seem to have a fair amount of activity towards both pst edges: near 0ms/baseline as well as after 400ms. This could add to the oscillation and convergence issues. I think a Hanning taper would be appropriate here. And at the same time in the tutorial docs it would be useful to illustrate what the canonical waveform for these data look like. Fig 1.2 show for all sensors, and a clear waveform, while present, is a bit hard to see.

3. Neither of the two models from the tutorial actually converge in the default 64 iterations of NLSI. Estimation halts before the convergence criteria are met. However, increasing the number of iterations is not possible through the GUI. This might be okay for just demo purposes.

4.. The model's states are in an oscillation for all sources, which may either be due to the pst edge activity or lack of convergence:

I suppose these issues are fine since this is a tutorial, but they can still be addressed to some extent to guide the user better.

[vlitvak]

Please do update the tutorial, and probably it'd be a good idea to update the data, removing Ashley's original invs altogether. Feel free to discuss the other points with Ashley, I can give you her email if you don't have it, or we could find out if she has a GitHub account and add her as a contributor.

[pranaysy]

Thanks, I'll do some more tests because there are ongoing issues with this particular model in SPM that are currently being fixed. I'll update once done.

Ashley's GitHub is here by the way: https://github.com/AshleyTyrer and the DCMs she has provided in her demo repository are the ~same as the SPM tutorial: they have large oscillations. I'll know more once I have fixed issues with the model itself, and then can involve her.

[pranaysy]

Hi @vlitvak
I now have a set of fixes for the tutorial.

• Fix 1 involves the lead field generation from SPM's canonical meshes from scratch, as in my comment above.
• Fix 2 concerns the NMDA model itself. Previously, it was not getting initialised correctly to its steady-state solution as documented in this issue. It is fixed with these commits on the main spm repo: https://github.com/spm/spm/pull/122/commits and the model now starts from a correct baseline before seeing input.
• Fix 3 concerns tapering of the data. There is enough activity at the edges that the model tries to explain them, resulting in oscillations. Models fit to Hanning-tapered data suppresses that edge activity and results in clear evoked activity that's damped properly.

After applying the above three fixes, the BMS results in the tutorials change quite a bit, with evidence favoring the 'full' model. The following are figures showing source activity with each fix, as well as the BMS results.

If things look fine, then the tutorial can be updated but I think only after the fixes to the model are shipped in the next SPM release. Otherwise users won't be able to reproduce the tutorial from the current release.

---

Here's activity in sources after fix 1 (baseline):

---

Here's activity in sources after fix 2, with oscillations that are more damped but still present:

---

Here's activity in sources after fix 3, showing nice evoked responses that return to baseline:

---

Finally, here are BMS results after applying the three fixes:

In the current tutorial, evidence for model 2 is ~200 greater than model 1, but that flips here.

---

[vlitvak]

Thanks, that's great! We'll be recommending people to use the latest development SPM, so it might make sense to update the tutorial and not wait. Just make sure Ashley is in the loop.