Forward/Inverse modeling sanity checks

💬 Support & Discussions
1

Objective

I’m playing around with forward/inverse modeling in MNE and am trying to check the assumption:

(1) kernel @ gain == I

Where kernel is the prepared inverse operator kernel, gain is the leadfield matrix, and I is an identity matrix. By extension, I want to check that:

(2) x == kernel @ gain @ x

Where x is the original time series data.

(2) is easier to check due to how MNE stores the leadfield matrix, but I need to obtain kernel and gain that satisfy (1) as well.

Assumptions

  • fwd is an MNE forward solution, in my case computed from a volume source space (BEM) with free orientations, 406 dipoles.
  • inv is the MNE inverse solution, computed using fwd and noise covariance cov from the raw data x.
  • x is stored as an MNE Raw object with nchan == 148

Implementation (2)

To compute the RHS of (2) above, I do the following:

src = mne.minimum_norm.apply_inverse_raw(x, inv, lambda2=1.0, method='dSPM', pick_ori='vector')
RHS = mne.apply_forward_raw(fwd, src, x.info)

At this point, I have to rescale RHS because I used dSPM, need to convert units from nAm to Am:

RHS = RHS.apply_function( lambda x: x * 1e-9 )

Comparing the results, it appears that RHS (orange) is pretty close to the original x, just de-trended:

image
image1161×585 94.3 KB

I imagine de-trending is part of running dSPM, haven’t checked the source to confirm.

Implementation (1)

Which brings me to (1). First, is the following the proper way to reconstruct the gain (leadfield) matrix?

eigen_leads = inv['eigen_leads']['data'].T
eigen_fields = inv['eigen_fields']['data'].T
sing = inv['sing']

gain = eigen_fields @ np.diag(sing) @ eigen_leads

If there were originally 148 channels and I used a volume source space with 406 dipoles,

gain.shape == (148, 1218)

Next, is this the proper way to obtain the kernel?

inv_prep = mne.minimum_norm.prepare_inverse_operator(inv, nave=1, 
                                                     lambda2=1.0,
                                                     method='dSPM')

kernel = mne.minimum_norm.inverse._assemble_kernel(inv_prep, label=None, 
                                                   method='dSPM',
                                                   pick_ori=None)[0]

Using the same parameters,

kernel.shape == (1218, 148)

At this point, I’m stuck:

kernel @ gain

Isn’t even close to being an identity matrix, even though the dimensions (148, 148) are correct.

2

As the problem is ill-posed this assumption will not hold:

to convince you rank(gain) is less than n_channels so rank(kernel @ gain) is less than n_channels while rank(I) == n_sources

Alex

1 Like
3

Thanks; I’ll have to read up on this some more.

Is the code I wrote the correct way to get the kernel and gain matrices?

4

yes it looks correct. To check just do a dot product with the data and check it matches the stc you would obtain with the public API

Alex

5

The units of dSPM are not nAm, they are a pseudo F-statistic. If you want physical units you need to use 'MNE' or 'eLORETA' methods, both of which will give SI units of Am. 'sLORETA' and 'dSPM' scale each source by a factor related to the noise level of each source (based on the noise cov).

1 Like