3.2. New Voxel-scale Connectivity Model [Knox2018]¶
The voxel-scale model from [Knox2018] is the first full brain voxel-scale of the mouse connectome. The model performs Nadaraya-Watson regression to infer the connectivity between each of the voxels in the brain into 12 major brain divisions to each of the voxels in the whole brain. The source space is split between these major brain divisions as to prevent influence from injections performed into adjacent brain divisions.
3.2.1. Assumptions¶
- Spatial smoothness within divisions: the connectivity is assumed to vary smoothly as a function of distance within each of the major brain divisions.
- No influence between divisions: the connectivity is allowed to be discontinuous at division boundaries. These major brain divisions are in fact physically separated by white matter, supporting this assumption.
3.2.2. VoxelConnectivityArray
Class¶
Since the full voxel x voxel connectivity matrix is ~200,000 x ~400,000 elements, it will mostlikely not fit in your memory. Luckily, the connectivity matrix has low rank structure, and we can take advantage of this by only computing the connectivty matrix on the fly, in the area we want to perform computation.
3.2.2.1. Loading the array¶
The easiest way to load the VoxelConnectivityArray
is through the
VoxelModelCache
object:
>>> from mcmodels.core import VoxelModelCache
>>> cache = VoxelModelCache(manifest_file='connectivity/voxel_model_manifest.json')
>>> voxel_array, source_mask, target_mask = cache.get_voxel_connectivity_array()
>>> voxel_array
VoxelConnectivityArray(dtype=float32, shape=(226346, 448962))
This downloads and caches the underlying data locally. Additionally, this loads the factorization of the connectivity matrix into memory (~1G total).
3.2.2.2. Retrieving values from the array¶
No part of the connectivity matrix is computed until the user asks for a value or set of values:
>>> # some given source/target voxels
>>> source, target = 20353, 68902
>>>
>>> # we index the VoxelConnectivyArray object just like it is a numpy ndarray
>>> connection_strength = voxel_array[source, target]
>>>
>>> # a row would be the bi-lateral connection strength from a given voxel
>>> row = voxel_array[source]
>>>
>>> # a column would be the connection strength to a given voxel
>>> # from each voxel in the right hemisphere
>>> column = voxel_array[:, target]
>>>
>>> # indexing the VoxelConnectivityArray object returns numpy ndarray
>>> type(row)
np.ndarray
If one wishes to operate on the full matrix (not recommended unless you have >1TB RAM!), computing the full matrix is similar to loading an hdf5 file:
>>> import sys
>>> full_array = vox_array[:]
>>> sys.getsizeof(full_array)
BIG!!!!!!
3.2.2.3. VoxelConnectivityArray methods¶
numpy.ndarray like methods
VoxelConnectivityArray
also has a few methods implemented from numpy.ndarray
.
These include:
and are called just like their numpy.ndarray
counter parts:
>>> voxel_array.T
VoxelConnectivityArray(dtype=float32, shape=(448962, 226346))
iterator methods
In addition to being able to index VoxelConnectivityArray
as a numpy.ndarray
,
VoxelConnectivityArray
implements several iterating methods:
iterrows
: yields single rowsitercolumns
: yields single columnsiterrows_blocked
: yields blocks of rows given the number of blocks.itercolumns_blocked
: yields blocks of columns rows given the number of blocks.
3.2.3. RegionalizedModel
class¶
Our voxel-scale model can be regionalized as well by integrating the connectivity matrix over some parcellation.
3.2.3.1. Metrics¶
Given a parcellation, integrating the connectivity over source and target regions gives us connection strength. Since the relative sizes of the regions may be vastly different, we can normalize this metric by dividing the connection strength by the size of ether the target region (connection density), the size of the source region (normalized connection strength) or by both the source and target (normalized connection density).
We can again use the VoxelModelCache
object to download and cache
regionalized voxel models:
>>> # this returns a pandas dataframe
>>> normalized_connection_density = cache.get_normalized_connection_density()
Alternatively, we could construct the RegionalizedModel
object using
our VoxelConnectivityArray
and our source/target Mask
objects:
>>> from mcmodels.models.voxel import RegionalizedModel
>>> # get set of summary structures
>>> structure_tree = cache.get_structure_tree()
>>> summary_structures = structure_tree.get_structures_by_set_id([167587189])[0]
>>> # the new ccf does not have sturcture 934 as a structure id
>>> structure_ids = [s['id'] for s in summary_structures if s['id'] != 934]
>>> # get keys
>>> source_key = source_mask.get_key(structure_ids=summary_structures)
>>> target_key = target_mask.get_key(structure_ids=summary_structures)
>>> regionalized_model = RegionalizedModel.from_voxel_array(
... voxel_array, source_key, target_key)
>>> normalized_connection_density = regionalized_model.normalized_connection_density
References
[Knox2018] | (1, 2) Knox et al. ‘High resolution data-driven model of the mouse connectome.’ bioRxiv 293019; doi: https://doi.org/10.1101/293019 |