amnsbr / cubnm

A toolbox for biophysical network modeling on GPUs

cuBNM toolbox is designed for efficient biophysical network modeling of the brain on GPUs and CPUs.

Overview

cuBNM toolbox simulates neuronal activity of network nodes (neural mass models) which are connected through the structural connectome using GPUs/CPUs. Currently three models (rWW, rWWEx and Kuramoto) are implemented, but the modular design of the code makes it possible to include additional models in future. The simulated activity of model neurons is fed into the Balloon-Windkessel model to calculate simulated BOLD signal. Functional connectivity (FC) and functional connectivity dynamics (FCD) from the simulated BOLD signal are calculated efficiently on GPUs/CPUs and compared to FC and FCD matrices derived from empirical BOLD signals to assess similarity (goodness-of-fit) of the simulated to empirical BOLD signal.

The toolbox supports parameter optimization algorithms including grid search and evolutionary optimizers (via pymoo), such as the covariance matrix adaptation-evolution strategy (CMA-ES). Parallelization within the grid or the iterations of evolutionary optimization is done at the level of simulations (across the GPU ‘blocks’), and nodes (across each block’s ‘threads’). The models can incorporate global or regional free parameters that are fit to empirical data using the provided optimization algorithms. Regional parameters can be homogeneous or vary across nodes based on a parameterized combination of fixed maps or independent free parameters for each node or group of nodes.

GPU usage is the primary focus of the toolbox but it also supports running the simulations on single or multiple cores of CPU. CPUs will be used if no GPUs are detected or if requested by the user.

At its core the toolbox runs highly parallelized simulations using C++/CUDA, while the user interface is written in Python and allows for user control over simulation configurations:

Documentation

Please find the documentations on installation, usage examples and API at https://cubnm.readthedocs.io.

cubnm's Issues

[DOC] Add tutorial on running a single simulation

Including:

• Selecting and visualizing SC
• Setting model parameters
• Running the simulation on GPU and CPU
• Visualizing simulated BOLD and model state variables time series
• Visualizing simulated FC and FCD

[ENH] Skip looping through non-existent connections

Currently in the integration loop tmpglobalinput to each node j is calculated by looping through all other nodes k and adding up SC[j,k] * S_E[k] even if SC[j,k] is zero. Skipping non-existent (zero) connections may increase efficiency, particularly in thresholded SCs.

[ENH] Option to use floats vs doubles

In the core code currently there is an unused preprocessor flag to "USE_FLOATS" which is a remnant of the legacy code (https://github.com/amnsbr/bnm_cuda). After making sure it works with the updated code of cuBNM, we might be able to use this as an installation option which can be controlled by setting an environment variable before pip install or by using pip install cubnm[float].

However, ideally we should be able to make this a runtime option so that it doesn't require a separate installation, but I can imagine that this can be very challenging.

It is also important to think if this option should really exist, and if it is really okay to use float precision in some simulations.

[BUG] Simulations not reproducible across some platforms

Using the development Docker container (both in Docker and Singularity) the generated simulations are different from the ones that are generated on other machines including Juseless, Colab and Kaggle, as well as the stable release container. This is probably not because of differences in compilation, since installing the exact same wheel on e.g. Juseless vs. the Docker container creates different simulations. The most likely culprit for the difference is the difference of random number sequences caused by the different libraries in each platform. From what I've learned the standard C++ random number generators are not guaranteed to be reproducible across platforms. Instead random number generators from e.g. boost library might be [ref]. Therefore using boost's random number generator may be a solution.

[DOC] Guide on how to add new models

Currently adding new models will require modifying the source code. Adding a new model involves adding different source and header files on C++/CUDA side as well as adding its Python interface. There should be a guide on how to do this.

[BUG] Core code exiting

Core code does exit(1) in some cases, and this will cause the Python program to exit as well. This happens for example when the simulation is too short and "Number of BOLD volumes (after removing initial volumes) is too low for FC calculations". Then especially if it is run within a Jupyter notebook, the error message printed from the core code is not shown to the user, and it will be unclear why the program has exited. Such cases should be handled more gracefully. Ideally Python should catch the errors and print them without exiting.

[ENH] User-defined integration dt

Currently neuronal model dt is fixed to 0.1 msec and BOLD dt is fixed to 1 msec. Users should be able to define custom dt's. They are currently hardcoded in the DerivedModel::Constants structs. In addition the integration loop expects BOLD dt to be 10x neuronal dt because there is an inner loop of 10 neuronal model integration steps within each step of BOLD model integration. Therefore, the steps to do this are:

• Make dt and bold_dt a model member variable.
• On Python side accept them as parameters, while ensuring they are divisible.
• Calculate n_inner as bold_dt / dt and adapt the integration loop accordingly.

[DOC] Add tutorial on running CMAES optimization of a heterogeneous model

Expand on currently existing tutorial at https://www.kaggle.com/code/aminsaberi/cubnm-0-0-2-demo-cmaes-heterogeneous.

• Include it in readthedocs
• Add more description
• Add an example of node-wise heterogeneity
• Compare resulting GOF with homogeneous model

[ENH] Add Jansen-Rit model

[BUG] Fixed issue in CPU-GPU identity tests revealed Kuramoto CPU-GPU test fails

In a8d2a9c and edf6e25 I replaced model-specific tests with a generic ./tests/sim/test.py file which runs the tests on all the current and future models without having to create a new test file for each new model. When I ran this generic test on rWW, rWWEx and Kuramoto, all 40 tests except one failed: test_identical_cpu_gpu[Kuramoto-do_delay:1]. However, previously all 10 tests of ./tests/sim/test_Kuramoto.py were successful. Then I realized that there was actually a silly bug in the previous test_identical_cpu_gpu functions: same SimGroup instance was used to run simulations on both CPU and GPU, then they were added to a dictionary sgs[force_cpu] where force_cpu was True or False, without making a copy of them. Therefore eventually when sgs[True] and sgs[False] were being compared, identical instances with identical data were being compared, including the GPU simulation data. This issue was fixed in the generic test.

[BUG] Linking to pre-existing GSL libraries compiled without `-fPIC`

When building from source setup.py searches for libgsl.a and libgslcblas.a in some known paths which will be linked in compilation of bnm.cu and run_simulations.cpp. If there are preexisting GSL libraries but they have not been compiled with -fPIC compilation fails with this error message:

  /bin/ld: /usr/lib/x86_64-linux-gnu/libgslcblas.a(xerbla.o): warning: relocation against `stderr@@GLIBC_2.2.5' in read-only section `.text'
  /bin/ld: /usr/lib/x86_64-linux-gnu/libgslcblas.a(xerbla.o): relocation R_X86_64_PC32 against symbol `stderr@@GLIBC_2.2.5' can not be used when making a shared object; recompile with -fPIC

To fix this the setup needs to identify if GSL libraries are compiled without -fPIC and if this is the case avoid using them.

[ENH] Return noise array to user

This would be useful for e.g. checking noise statistics (e.g. how much auto-correlated they are, or how much they are correlated across nodes) and checking how much it influences model state variables.

[ENH] Create Docker container

This should hopefully make the toolbox available across more platforms, including Windows.

[ENH] Enable band/high/low-pass filtering of the simulated BOLD

This will be tricky to implement on GPU but can be done on C++/Python side, i.e. simulated BOLD is calculated on GPU, returned to CPU [C++/Python] to be filtered, and then FC and FCD are again calculated on GPU. To do this the mean_bold and ssd_bold calculations should not be done within the bnm kernel (which is the current behavior), but rather should be done in a separate kernel after filtering.

[ENH] Add Kuramoto model

Add Kuramoto phase oscillator model as a simple initial example of oscillator models. This will also require adding support for calculating global input based on phase.

[DOC] Add tutorial on running CMAES optimization of a homogeneous model

Expand on currently existing tutorial at https://www.kaggle.com/code/aminsaberi/cubnm-0-0-2-demo-cmaes-homogeneous.

• Include it in readthedocs
• Add more description

[BUG] Import error with installing from wheels: `undefined symbol: shm_open`

On some Linux machines the wheels can be installed successfully from pip (therefore are compatible with manylinux_2_28) but on import cuBNM raise the following error:

<path-to-virtual-env>/lib/python3.10/site-packages/cuBNM/core.cpython-310-x86_64-linux-gnu.so: undefined symbol: shm_open

In current implementation any ImportError including undefined symbol in the error message are followed by a recommendation to build the package from source, but the root of the problem needs to be fixed.

Also the original ImportError should be raised instead of being printed.

[ENH] Command line interface

Currently the toolbox can only be used as a module that can be imported. But some users probably would benefit more from a command line interface. Also this will be a good entry point when we containarize the toolbox. There could be different commands for running one/several specific simulations with defined parameters (e.g. cubnm run), for running a grid search (e.g. cubnm grid) or for running an evolutionary optimization (e.g. cubnm optimize --optimizer CMAES).

[ENH] Optimizer batch mode

Currently only one optimizer run can be executed at any given time. Sometimes (especially on GPUs) it would be more efficient if multiple optimizer runs are executed in parallel on the same GPU to use the full capacity of the device. These multiple runs can be different runs (with different seeds) of the same input data (i.e., same subject), or on different input data.

[ENH] Better control on heterogeneity maps

This will involve the following changes:

• Allow user to use any custom numpy array as the heterogeneity maps (rather than accepting a text file as an input)
• Make the selection of preloaded maps in datasets more intuitive. For example, rather than '6maps', accept a list of map names, e.g. ['myelinemap', 'fcg1'].

[ENH] Enable modification of some Balloon-Windkessel parameters

Some of the Balloon-Windkessel model parameters (k1, k2 and k3) are dependent on scanner field strength (see Henizle 2016 and e.g. Demirtaş 2019 using this approach). The user should be able to choose the source of these model parameters, and in case of heinzle2016 choose the parameters associated with the scanner field strength of the empirical data.

[DOC] Add tutorial on running a grid search

Expand on currently existing tutorial at https://www.kaggle.com/code/aminsaberi/cubnm-0-0-2-demo-grid.

• Include it in readthedocs
• Add more description
• Add grids of model state variable averages in addition to GOF

[ENH] Add Wilson-Cowan model

[BUG] Compilation warnings

This is not per se a bug, but ideally there shouldn't be many (or any) compilation warnings when the CUDA or C++ code is being compiled, and currently there are many. They should be investigated and resolved one by one.