brain-modelling-group/neural-flows

Estimation, analysis and decomposition of brainwave spatiotemporal dynamics

Neural Flows: estimation of velocities, detection and classification singularities in 3D+t brain data

If you use this toolbox, please cite:

Brief Description

This is neural-flows, a library to estimate wave velocities and singularities in 3D+t in neuroimaging data.
The figure below gives you an idea of a typical workflow you can run with this toolbox:
Each link of the path is a component of the core/ namely:

• data interpolation,
• flow estimation,
• singularity detection and classification, and
• streamline tracing.

The toolbox has a number of ancillary modules, including:

• io/ that handles read and write operations of .json and .mat files.
• visualisation/ contains lots of plotting functions including some that produce movies like the one you can see at the end of this document.
• external/, which contains code implemented by others, the most prominent being C-NEM invoked by so-called meshless methods, and the Hungarian method used track singularities; and finally,
• utils/ which gathers small standalone functions that support the main modules.

Getting started

Dependencies:

MATLAB

• tested on R2018b & R2020a
• in principle should work on any OS that MATLAB supports; tested on Windows 10, Debian 10, Ubuntu 20.04, Ubuntu 16.l4, OpenSuse 15.1

CNEM, https://m2p.cnrs.fr/sphinxdocs/cnem/index.html

• only relevant if you want to use the phase-gradient method for calculating flows and meshless method for tracing streamlines.
• meshless methods speed up calculations significantly but installing cnem can be tricky.
• used version (v03, 2014-05-04) available here: https://ff-m2p.cnrs.fr/frs/?group_id=14
• aka cnemlib: https://au.mathworks.com/matlabcentral/linkexchange/links/3875
• if you have trouble, contact the author of the library (mail can be found on the cnrs.fr website linked above) who has been super helpful when we had issues.

Installation:

If you download the zip folder, first unzip and then ...

• Windows 10: see the instructions here.
• Linux: on the terminal change directory to neura-flows folder, then start matlab. That's it. The toolbox paths should be appended automatically.

If you clone the repo, same rules as above apply.

What can the toolbox do?

A bunch of things, mostly the ones described in the diagram above.

To get started follow the next steps (timestamp: Tue 29 Jun 2021 14:49:38 AEST):

Out of the box examples

0. (Optional) Inside directory neural-flows/, make a new directory called scratch/ (by default we will store results there, you can change it later ...). If this directory doesn't exist the toolbox will print a warning message and then proceed to create it.

1. Call the following function from within the top level directory neural-flows/

fmri_waves()

or

rotating_wave()

The first example takes about 7 minutes to run everything and then it pops up figures with flow mode decomposition and singularity tracking! The second one takes about 9 minutes. The execution times are estimated using a high-end workstation (circa 2018). On a laptop it can take 2.5 longer to run (o.O it's heavy stuff).

Flow mode decomposition

Singularity statistics

Movies

Movies are not produced by default as they require more computational resources. However, we provide the visualisation tools to make movies like this one:

What if I want to use my own data? Input data format

The toolbox reads data from a .mat file. It will expect the following variables to be stored.

- `data`: a 2D array of size [number of timepoints x number of ROIs/nodes/sources] with the neuroimaging data for unstructured datasets (i.e., source reconstructed EEG/MEG data, whole-brain models/brain networks, parcellated fmri data) 
- `locs`: a 2D array of size [number of ROIs/nodes/sources x 3] with the x,y,z coordinates of ROIs/nodes/sources, expressed in mm or m.
- `ht`  : sampling period expressed in ms or s.

neural-flows's Issues

Write startup function

Silence CNEM console output

CNEM prints lots of messages to matlab's workspace and to the terminal.
Check if we can silence them.

Required:

• rebuild cnem and mex files

Document all functions consistently

Make sure all functions follow the template for documentation input/output arguments, requirements

Enable singularity tracking as part of core functionality

1D singularity tracking is functional but not integrated into main_core()

Consolidate generation of streamlines seeding locations

Currently there are three/four functions that generating random seeding locations. Main feature is the density of seeding locs. Needs to handle 2D/3D space and structured/unstructured domains.

Reconsider how to handle CNEM output files

CNEM writes binary files at every (time) iteration. This may increase execution times, in particular in server-class machines with HDDs.

TODO:

• suppress output.
• benchmark with/without writing to disk.

Parallelise streamlines calculation on structured grids

Currently working functionality but sequentially.
Need to implement parallel version similar to the way interpolation of data works.

Implement Helmholtz decomposition using curl & divergence

This decomposition may be useful for classification of patterns:

Implementation is required for:

• 3D structured
• 2D structured and 2D unstructured

Implement HS for unstructured grids in 2D unstructured domains

Is your feature request related to a problem? Please describe.
This is a new feature

Describe the solution you'd like
An implementation of HS algorithm to operate on 2D unstructured domains embedded in 2D

Describe alternatives you've considered
Nothing available yet

Debug tests

Some tests are not passing because of changes in the input .json file.

Implement in Python

Yeah ... I said that

Enable struct use instead of file objects

neural-flows uses a type of matfile that allows for partial i/o operations. With current SSD w/r speeds it shouldn't be a problem. This choice was made to keep memory usage to minimum. The disadvantage is that this feature may bring an HPC down if the storage directory is on a network drive rather than scratch space.

Describe the solution you'd like
For small datasets, use a struct() with equivalent hierarchy and organisation as the current file_object. Keep everything in memory and write at the end.

Describe alternatives you've considered
Memory mapped files, but they don't work very well with parallelised loops.
Pick any two: minimise bookkeping complexity, enable parallel computations, minimise memory usage.

meshless streamline tracing does not trace properly

Describe the bug
Streamlines tracing using CNEM produce nonparallel streamlines for a perfect travelling wave

To Reproduce
Steps to reproduce the behavior:

1. Run travelling_reflected_wave('u-a-hs-mb') and then travelling_reflected_wave('u-a-hs-ml')

2. Streamlines should be parallel and look like a collection of vertical lines on the x-y plane. This is achieved with mesh-based methods but not for meshless streamlines tracing.

Expected behavior
both methods should produce the same results

Screenshots

Desktop (please complete the following information):

• OS: Ubuntu
• Matlab version [e.g. 2018b]
• Device specs

(RAM, processor)

Additional context
Add any other context about the problem here.

Generalise routing analysis

Streamline tracing for phase-based cnem flows fails

Describe the bug
Streamline tracing fails for rotating_wave('upc') - phase-based flows using CNEM

To Reproduce

1. From neural-flows/examples, run:

rotating_wave('upc') 

Expected behavior
Execute without errors

Screenshots

interpoltaion ok...

Error using  * 
Incorrect dimensions for matrix multiplication. Check that the number of columns in the first matrix matches the number of rows in the second matrix. To perform
elementwise multiplication, use '.*'.

Error in m_cnem3d_interpol/interpolate (line 134)
            Var_Int=obj.Mat_INT*[Var;Var_Add_Node];

Error in streams3d_tracing_cnem/tracing_cnem_loop (line 123)
        flowfield_interp = cnem_interpol_obj.interpolate(flow_field); % has same length as number of live points

Error in streams3d_tracing_cnem/tracing_cnem_phase_step (line 172)
        output_cell = tracing_cnem_loop(locs, boundary_faces, flow_field, seeding_locs, time_step, max_stream_length);

Error in streams3d_tracing_cnem/tracing_cnem_time_loop (line 55)
       streamlines.paths = tracing_cnem_fun_step(tt);

Error in streams3d_tracing_cnem (line 52)
tracing_cnem_time_loop()

Error in streams3d_trace (line 55)
            streams3d_tracing_cnem(obj_streams, obj_flows, params)

Error in main_core (line 62)
    streams3d_trace(tmp_params)

Error in rotating_wave (line 36)
output_params = main_core(input_params);
 

Desktop (please complete the following information):

• OS: Ubuntu 20.04
• Matlab version: 2018b
• Device specs: tesseract

Generalise consistent interface for plotting functions

• Decide on consistent naming convention
• Split functions into high- & low-level plotting functions
• Define interface for the two groups

High-level plotting functions:

• Interface: input: params | output: none or figure_handle

Low-level plotting functions:

• Interface: input data and graphics properties | output: axes_handles

Write readme with installation/dependencies and basic example

Disable parpool IdleTimeout until results are collated into a single file

If IdleTimeout is reached, the parapool shuts down and the Compiste files with partial results become accessible, an we see the following error

Error using Composite/subsref (line 22)
The Composite is invalid; it may have been saved and loaded, which is not allowed, or the parallel pool has
been shut down.

Error in spmd_parfor_with_matfiles (line 78)
       worker_fname = WorkerFname{this_temp_file};

Error in streams3d_tracing_mlab_parallel (line 41)
    [obj_interp] = spmd_parfor_with_matfiles(tpts, parfun, obj_streams, streamtracing_3d_storage_expression);

Error in streams3d_trace (line 58)
                streams3d_tracing_mlab_parallel(obj_flows, obj_streams, params)

Error in main_core (line 62)
    streams3d_trace(tmp_params)

Error in fmri_waves (line 50)
output_params = main_core(input_params);