Dipole scanning - array exceeds maximum array size preference

Discussions
,
1

Dear BST community,

I am trying to perform LCMV beamforming on some already preprocessed continuous resting state EEG data (512 Hz) with the ICBM152 template. My goal is to extract the source reconstructed time series using scouts from the DK atlas.

Accordingly, due to using the MNI template, as stated here https://neuroimage.usc.edu/brainstorm/Tutorials/SourceEstimation, I used unconstrained solutions:

Unconstrained solutions are particularly appropriate when using the MNI template

I ultimately want one signal per scout. Now, I am a bit confused as to what the next step would be.

My first thought was to use Extract Values, with the Scout function set to PCA, but this still provides 3 signals per scout. I want the signals to oscillate around 0, thus cannot use the "norm for unconstrained sources" from my understanding.

Going over the tutorials once again, it looks like Dipole Scanning (https://neuroimage.usc.edu/brainstorm/Tutorials/TutDipScan) is what I want.

I tried it, but the issue my files are too large. I got the following warning:

Requested 45006x245761 (82.4GB) array exceeds maximum array size preference (64.0GB). This might cause MATLAB to become unresponsive.

All my files are around this size. I am wondering what should I do in this case. Would downsampling help?

Another option is "Unconstrained to flat map" using PCA, but this also runs into the same array size issues.

Thank you,
Paul

2

My first thought was to use Extract Values, with the Scout function set to PCA, but this still provides 3 signals per scout.

The option "Scout function: PCA" groups the various signals per region of interest: 1 signal per ROI for constrained orientations, 3 signals per ROI for unconstrained sources.
https://neuroimage.usc.edu/brainstorm/Tutorials/Scouts#Scout_function

Working with unconstrained source maps makes everything much more complicated, and not all the processing options are available. Using constrained orientations helps decreasing the complexity of the problem, but this simplification is difficult to justify in the case of a template brain.

Another option is "Unconstrained to flat map" using PCA, but this also runs into the same array size issues.

If you want to convert from 3 signals to 1 per time point, you may use the process "Unconstrained to flat map", indeed. But it requires reconstructing the full source matrix (45006x245761) instead of working in the optimized format (kernel=45000xNchannels + data=Nchannels x 250000), you won't be able to run this on such large file, as you already noticed.

Additionally, there is no better justification for using this kind of simplification instead of using an arbitrary source orientation. This approach has not been published anywhere and is not part of our recommended workflows.

Going over the tutorials once again, it looks like Dipole Scanning (https://neuroimage.usc.edu/brainstorm/Tutorials/TutDipScan) is what I want. I tried it, but the issue my files are too large. I got the following warning: Requested 45006x245761 (82.4GB) array exceeds maximum array size preference (64.0GB).

The dipole scanning is designed to work on a few time points, in the context of an ERP/ERF study, not to process long time series.

All my files are around this size. I am wondering what should I do in this case. Would downsampling help?

Processing full resolution cortex surfaces over long recordings is not a solution, you need to simplify your analysis somewhere: reducing the number of time points, reducing the number of sources, working with regions of interest, working with summary metrics at the sensor level.

We optimized the computation of the PSD from resting state recordings at the source level:
https://neuroimage.usc.edu/brainstorm/Tutorials/RestingOmega
Optimization for other metrics are possible, but may require that you code them yourself.

Remember that all the information is available in your EEG recordings, the source estimation does not create any extra information, it just helps localizing it spatially. Before trying to process everything in source space, make sure you can already observe what you expect at the sensor level.

3

Thank you Francois.

Sorry for what is a naive question (this is my first time working with unconstrained source maps), but if my goal is to ultimately have a single time series from each of the 68 scouts of the DK atlas, do you see any way that could help reduce the data size somewhere in the pipeline?

I also downsampled my data to 256 Hz, and the "Unconstrained to flat map" worked (in a somewhat reasonable time window), but as you mentioned, there is no better justification for this approach.

Thank you again.

4

Sorry for what is a naive question (this is my first time working with unconstrained source maps), but if my goal is to ultimately have a single time series from each of the 68 scouts of the DK atlas, do you see any way that could help reduce the data size somewhere in the pipeline?

Downsampling (time dimension), reducing the number of sources (= number of vertices).

To keep the full dimensions: computing the inverse kernel only, then write a Matlab script to write a loop to compute your measure of interest on only one scout at a time without having to generate the full source matrix. This is not trivial.

You could also compute the measure on the three dimensions independently, and then group them in some way at the end...

5

Hello,
I am a newbie here, so I apologise in advance if my question has been answered previously. I think I have a similar problem for my connectivity analysis. I am trying to understand how the choice of the scout function and the aggregation procedure would affect MEG connectivity results. I managed to run the ''scoutime' 1 option (before) for PCA and Mean scout functions, over epoched and preprocessed resting state data, but matlab seems not to be able to handle the after scoutime option. I know I could reduce the N of scouts but my ultimate goal is to do whole brain network analysis.
I guess I'd need to simplify my data, by some sort of averging over scouts time series(?), but I wanted to hear your opinion first. Please let me know if I am missing something obvious.

Here is the piece of my code:

**% Process: Envelope Correlation NxN [2020]**

sFiles_Atla1_Alpha_AverageAfter = bst_process('CallProcess', 'process_henv1n', sFiles, [], ...
'timewindow', [], ...
'scouts', {'Schaefer_100_17net', {'Background+FreeSurfer_Defined_Medial_Wall L', 'Background+FreeSurfer_Defined_Medial_Wall R', 'ContA_IPS_1 L', 'ContA_IPS_1 R', 'ContA_PFCl_1 L', 'ContA_PFCl_1 R', 'ContA_PFCl_2 L', 'ContA_PFCl_2 R', 'ContB_IPL_1 R', 'ContB_PFCld_1 R', 'ContB_PFClv_1 L', 'ContB_PFClv_1 R', 'ContB_Temp_1 R', 'ContC_Cingp_1 L', 'ContC_Cingp_1 R', 'ContC_pCun_1 L', 'ContC_pCun_1 R', 'ContC_pCun_2 L', 'DefaultA_IPL_1 R', 'DefaultA_PFCd_1 L', 'DefaultA_PFCd_1 R', 'DefaultA_PFCm_1 L', 'DefaultA_PFCm_1 R', 'DefaultA_pCunPCC_1 L', 'DefaultA_pCunPCC_1 R', 'DefaultB_IPL_1 L', 'DefaultB_PFCd_1 L', 'DefaultB_PFCd_1 R', 'DefaultB_PFCl_1 L', 'DefaultB_PFCv_1 L', 'DefaultB_PFCv_1 R', 'DefaultB_PFCv_2 L', 'DefaultB_PFCv_2 R', 'DefaultB_Temp_1 L', 'DefaultB_Temp_2 L', 'DefaultC_PHC_1 L', 'DefaultC_PHC_1 R', 'DefaultC_Rsp_1 L', 'DefaultC_Rsp_1 R', 'DorsAttnA_ParOcc_1 L', 'DorsAttnA_ParOcc_1 R', 'DorsAttnA_SPL_1 L', 'DorsAttnA_SPL_1 R', 'DorsAttnA_TempOcc_1 L', 'DorsAttnA_TempOcc_1 R', 'DorsAttnB_FEF_1 L', 'DorsAttnB_FEF_1 R', 'DorsAttnB_PostC_1 L', 'DorsAttnB_PostC_1 R', 'DorsAttnB_PostC_2 L', 'DorsAttnB_PostC_2 R', 'DorsAttnB_PostC_3 L', 'LimbicA_TempPole_1 L', 'LimbicA_TempPole_1 R', 'LimbicA_TempPole_2 L', 'LimbicB_OFC_1 L', 'LimbicB_OFC_1 R', 'SalVentAttnA_FrMed_1 L', 'SalVentAttnA_FrMed_1 R', 'SalVentAttnA_Ins_1 L', 'SalVentAttnA_Ins_1 R', 'SalVentAttnA_Ins_2 L', 'SalVentAttnA_ParMed_1 L', 'SalVentAttnA_ParMed_1 R', 'SalVentAttnA_ParOper_1 L', 'SalVentAttnA_ParOper_1 R', 'SalVentAttnB_IPL_1 R', 'SalVentAttnB_PFCl_1 L', 'SalVentAttnB_PFCl_1 R', 'SalVentAttnB_PFCmp_1 L', 'SalVentAttnB_PFCmp_1 R', 'SomMotA_1 L', 'SomMotA_1 R', 'SomMotA_2 L', 'SomMotA_2 R', 'SomMotA_3 R', 'SomMotA_4 R', 'SomMotB_Aud_1 L', 'SomMotB_Aud_1 R', 'SomMotB_Cent_1 L', 'SomMotB_Cent_1 R', 'SomMotB_S2_1 L', 'SomMotB_S2_1 R', 'SomMotB_S2_2 L', 'SomMotB_S2_2 R', 'TempPar_1 L', 'TempPar_1 R', 'TempPar_2 R', 'TempPar_3 R', 'VisCent_ExStr_1 L', 'VisCent_ExStr_1 R', 'VisCent_ExStr_2 L', 'VisCent_ExStr_2 R', 'VisCent_ExStr_3 L', 'VisCent_ExStr_3 R', 'VisCent_Striate_1 L', 'VisPeri_ExStrInf_1 L', 'VisPeri_ExStrInf_1 R', 'VisPeri_ExStrSup_1 L', 'VisPeri_ExStrSup_1 R', 'VisPeri_StriCal_1 L', 'VisPeri_StriCal_1 R'}}, ...
'scoutfunc', 1, ... % Mean
'scouttime', 2, ... % After
'removeevoked', 0, ...
'tfmeasure', 'hilbert', ... % Hilbert transform
'edit', struct(...
    'Comment', 'Test', ...
    'TimeBands', [], ...
    'Freqs', {{'alpha', '8, 12', 'mean'}}, ...
    'ClusterFuncTime', 'none', ...
    'Measure', 'none', ...
    'Output', 'all', ...
    'RemoveEvoked', 0, ...
    'SaveKernel', 0), ...

'tfsplit', 1, ...
'cohmeasure', 'oenv', ... % Envelope correlation (orthogonalized)
'statdyn', 'static', ... % Static
'win_length', 1.5, ...
'win_overlap', 50, ...
'parallel', 0, ...
'outputmode', 2);

and the matlab error: "Requested 900x225060004 (3018.3GB) array exceeds maximum array size preference (125.6GB). This might cause matalab to become unresposive."

Many thanks in advance.

D.

6

This is discussed in the new tutorials:

1 Like
7

Thank you so much Francois. I've read the tutorials and understand the difference between before/after. However, I'd still like to be able to compute the oenv with the option scoutfunc 1; scoutime 2. Do you reckon this is possible at all?

8

The error message tells you'd need 3018.3GB of memory.
Therefore: no, you can't

1 Like