[FieldTrip] Issue plotting Brainnetome aligned EEG source data

[Jack Fogarty]

Hi Jan Mathijs,

Thanks for your consideration. I have some responses to some of your points below:

1. "It is a bit unclear to me what the relevance of the atlas is in all this."

Ultimately, I am looking to explore EEG connectivity between brainnetome ROIs. Specifically, in my full pipeline I am:
a) Creating Brainnetome atlas aligned template grid
b) Using the template grid for source reconstruction across all EEG trials (i.e. epochs)
c) Muliplying the spatial filter with the individual EEG epochs to get trial-level source data
d) Projecting the dipole moments to their max in each trial
e) Using ft_volumelookup to select and then average the dipole mom for each brainnetome ROI to create a 'virtual channel' timecourse
f) Running ft_freqanalysis and ft_connectivity analysis on the virtual channel data

2. "The atlas should not be used a priori to define the ‘inside’ compartment. You can map the dipoles onto the anatomical parcels after you performed the source reconstruction."

Is this still the case if I am using a standard BEM model? I decided to create a template grid pre-aligned to the Brainnetome atlas to make it simple to compare subject source-level data and interpret the outcomes relative to specific ROIs in the brainnetome atlas. Is this a poor approach? 

I did first try the typical 'source reconstruction -> interpolation to atlas -> source parcellate' pipeline but I ran into a lot of problems trying to interpolate my source reconstruction onto the Brainnetome atlas. Further digging led me to the current approach as a possible solution.

3. "Following up on my previous message: I didn’t think of the fact that you are using EEG, with a BEM volume conduction model. In this case the inwardshift is probably not going to work."

You are right there, as far I can figure. I played a lot with various inwardshift, resolution, and other settings when I was preparing the template_grid (sourcemodel), and I also mucked around with ft_sourceplot projection methods. In some cases, certain options did look like they improved the surface plot but I think they were only hiding some of the problem. Inward shift only appeared to cover the surface when I had shifted many gridpoints 'outside' the brain volume; although, regardless, as soon as I align the template grid to the brainnetome atlas the problem returns.

4. "Is there a particular reason why you’d want to visualize your results on the surface?"

Self-gratification, perhaps. Really, I was hoping to plot the outcomes on a surface because it is a nice and simple way to display results, both for checking that things 'look sensible' and for potential publications where we might use this atlas. Is it the data type that makes this tricky? I am wondering if my pipeline is ok, but perhaps I need to give up on the surface plot, or create one specifically to represent data at ROIs of the Brainnetome atlas. 

In case it is helpful for you, or for other interested readers, my full code for steps a-e is below. Thanks again for your thoughts!

Best,
Jack



%% Fieldtrip: eLORETA ROI timecourses
% https://www.fieldtriptoolbox.org/example/create_single-subject_grids_in_individual_head_space_that_are_all_aligned_in_brain_atlas_based_mni_space/
% https://www.fieldtriptoolbox.org/tutorial/salzburg/#reconstruct-the-time-course-of-activity-at-a-particular-brain-location

    ft_defaults
    [ftver, ftpath] = ft_version;

    % Load the Headmodel
    load(fullfile(ftpath, 'template/headmodel/standard_bem.mat')); % Colin27

    % Prepare the template grid 'standard sourcemodel'
    cfg = [];
    cfg.xgrid  = 'auto';
    cfg.ygrid  = 'auto';
    cfg.zgrid  = 'auto';
    cfg.unit   = 'mm';
    cfg.tight  = 'yes';
    cfg.resolution  = 7; % What does this result in regard to the final resolution for reporting?
    cfg.headmodel   = vol;
    template_grid   = ft_prepare_sourcemodel(cfg);
    
    % Check alignment of the template grid and headmodel
    figure;
    hold on
    ft_plot_mesh(template_grid.pos(template_grid.inside,:));
    ft_plot_headmodel(vol,  'facecolor', 'cortex', 'edgecolor', 'none');
    ft_plot_axes(vol);
    alpha 0.5
    camlight
    
    % Read in Brainnetome Atlas; ensure its units are consistent with template
    atlas = ft_read_atlas(fullfile(ftpath, 'template/atlas/brainnetome/BNA_MPM_thr25_1.25mm.nii'));
    
    % Ensure units are consistent between atlas and sourcemodel
    atlas = ft_convert_units(atlas,'mm');
    template_grid = ft_convert_units(template_grid,'mm');
    
    % Check the coordsys field exists in sourcemodel
    if ~isfield(template_grid, 'coordsys');
        template_grid.coordsys = 'mni';
    end

    % Create binary mask for all grid points inside the Brainnetome atlas locations
    cfg = [];
    cfg.atlas      = atlas;
    cfg.roi        = atlas.tissuelabel;  % here you can also specify a single label, i.e. single ROI
    cfg.inputcoord = 'mni';
    mask           = ft_volumelookup(cfg, template_grid);
    template_grid.inside = false(template_grid.dim);
    template_grid.inside(mask == 1) = true;

    % Check grid points inside the Brainnetome atlas
    figure;
    ft_plot_mesh(template_grid.pos(template_grid.inside,:),'vertexcolor','black','vertexalpha',0.5);
    
    clearvars -except atlas template_grid vol
    
    %% Prep EEG source analysis
    
    setenv('PATH', 'C:\Program Files\OpenMEEG\bin');
    setenv('LD_LIBRARY_PATH', 'C:\Program Files\OpenMEEG\lib'); 
    
    % Load preprocessed EEG data
    load('C:\Users\jf752\Desktop\Processing\source_analysis\test_data');
    
    % Realign electrodes
    cfg = [];
    cfg.method = 'interactive';
    cfg.headshape = vol.bnd(1);
    cfg.elec = cln_data.elec;
    elec_aligned = ft_electroderealign(cfg);

    cln_data.elec = elec_aligned;
    
    % Convert time for each trial so that they are equivalent (for rawtimelock)
    for i = 1:length(cln_data.time);
        cln_data.time{1,i} = 0:0.002:1.998;
    end
    
    clear i elec_aligned cfg
    
    % Compute sensor-level time-locked data
    cfg = [];
    cfg.trials = 'all';
    cfg.keeptrials = 'yes';
    cfg.covariance = 'yes';
    
    timlock = ft_timelockanalysis(cfg,cln_data);

    % Prepare the leadfield
    cfg = [];
    cfg.elec      = cln_data.elec;            
    cfg.channel   = cln_data.label;
    cfg.headmodel = vol; % volume conduction model   
    cfg.sourcemodel.pos    = template_grid.pos; % number of voxels

    leadfield = ft_prepare_leadfield(cfg);
    
    clear cfg
    
    %% Run source localisation and calculate single-trial ROI timecourses
    
    % Compute eLORETA time-course based on the total data
    cfg = [];
    cfg.method       = 'eloreta';
    cfg.elec         = timlock.elec;
    cfg.sourcemodel  = leadfield;
    cfg.headmodel    = vol;
    cfg.projectnoise = 'yes';
    cfg.keepmom      = 'yes';
    cfg.keepfilter   = 'yes';
    
    source = ft_sourceanalysis(cfg, timlock);

    % Calculate 3D trial moms based on Michael Glassen's FTMailing Request
    % https://mailman.science.ru.nl/pipermail/fieldtrip/2020-May/039989.html
    for f = 1:length(source.avg.filter);
        for t = 1:length(cln_data.trial);
            if isempty(source.avg.filter{f});
                x = [];
            else
                x = source.avg.filter{f}*cln_data.trial{t};
            end
        trial_source.trial(t).mom{f} = x;
        end
    end
    
    % Mimic Source.avg config to project trial dipole moments in direction of maximal power
    tmp = source;
    
    for t = 1:length(cln_data.trial);
        tmp.avg.mom = trial_source.trial(t).mom;
        
        cfg = [];
        cfg.projectmom = 'yes';
        tmp = ft_sourcedescriptives(cfg, tmp);
        
        trial_source.trial(t).mom = tmp.avg.mom;
        trial_source.trial(t).pow = tmp.avg.pow;
        trial_source.trial(t).ori = tmp.avg.ori;
        trial_source.trial(t).label = tmp.avg.label;
        trial_source.trial(t).filterdimord = tmp.avg.filterdimord;        
    end
        
    clear f t x tmp ans cfg
    
    % Calculate the average signal in each ROI
    for t = 1:length(cln_data.trial);
        for r = 1:length(atlas.tissuelabel);
            cfg = [];
            cfg.atlas      = atlas;
            cfg.roi        = atlas.tissuelabel{1,r};  % here you can also specify a single label, i.e. single ROI
            cfg.inputcoord = 'mni';
            mask           = ft_volumelookup(cfg, template_grid);

            mask_idx = find(mask);

            for m = 1:length(mask_idx);
                if isempty(trial_source.trial(t).mom{mask_idx(m)});
                   tmp(m,:) = nan(1,length(cln_data.trial{1}));
                else
                   tmp(m,:) = trial_source.trial(t).mom{mask_idx(m)};
                end
            end

            roi(r,:) = nanmean(tmp);
            clear tmp mask mask_idx ans m
        end
        roi_trial_dat{t} = roi;
        clear roi r
    end