Structural network mapping¶

In this tutorial, you will explore the structural network mapping functionalities of Lacuna using the CLI.

What you'll learn:

• Fetch a structural connectome
• Obtain voxel-level disconnectivity maps
• Compute a disconnectivity matrix
• Compute ROI-wise disconnection values

Setup¶

Get the tutorial data.

In [ ]:

remove_output

!lacuna tutorial /tmp/tutorial_bids --force

!lacuna tutorial /tmp/tutorial_bids --force

Check whether mrtrix3 is properly installed as it is required for the analysis.

In [ ]:

remove_output

!mrinfo --version

!mrinfo --version

Fetch structural connectome¶

Structural lesion network mapping requires a structural connectome in the form of a normative tractogram, i.e., white matter fiber bundle reconstruction derived from diffusion-weighted MRI.

When selecting a connectome, the main difference lies in the trade-off between sensitivity and anatomical specificity:

If possible, we recommend using dTOR985 as its high coverage improves sensitivity. It is generated via an untargeted, whole-brain tracking approach. With ~12 million streamlines, this connectome is computationally demanding. The fetching process involves a conversion to .tck format which needs 32gb RAM. To download it you need a Figshare personal token. You can obtain this token by logging into Figshare, navigating to Account, then Integrations, and finally selecting Personal Tokens.

In this tutorial, we use the more lightweight HCP1065 tractogram which is a combination of individual curated tract bundles. It uses trajectory-based recognition against an expert-vetted atlas and topology-informed pruning to remove false connections. This makes it suitable for analyses where minimizing spurious (false-positive) streamlines is a priority. However, its lower coverage lowers sensitivity.

In [ ]:

remove_output

!lacuna fetch hcp1065 \
    --output-dir /tmp/hcp1065_data

!lacuna fetch hcp1065 \ --output-dir /tmp/hcp1065_data

Analysis¶

By filtering the normative tractogram based on the mask image, structural lesion network mapping identifies white matter fiber bundles connected to a specific lesion. By that, we can infer which pathways have been disconnected or disrupted.

What the structural network mapping implementation internally does:

1. Ensure spatial alignment
2. Compute the voxel-level tract density image (TDI) of full tractogram (amount of streamlines per voxel)
3. Filter the normative tractogram by the lesion mask
4. Compute the TDI of the lesion filtered tractogram and divide the lesion TDI by the original TDI to obtain a disconnection map (% of streamlines connected to the lesion).

Now run the analysis.

In [ ]:

remove_output

!lacuna run snm \
    /tmp/tutorial_bids/ \
    /tmp/outputs_snm_run1/ \
    --connectome-path /tmp/hcp1065_data/hcp1065.tck \
    --participant-label 01 \
    --mask-space MNI152NLin6Asym

!lacuna run snm \ /tmp/tutorial_bids/ \ /tmp/outputs_snm_run1/ \ --connectome-path /tmp/hcp1065_data/hcp1065.tck \ --participant-label 01 \ --mask-space MNI152NLin6Asym

The standard run as shown above produces voxel-level disconnectivity maps.

In [7]:

!ls /tmp/outputs_snm_run1/sub-01/ses-01/anat/

!ls /tmp/outputs_snm_run1/sub-01/ses-01/anat/

sub-01_ses-01_desc-provenance.json
sub-01_ses-01_label-acuteinfarct_method-snm_desc-summarystatistics_stats.json
sub-01_ses-01_label-acuteinfarct_method-snm_space-MNI152NLin6Asym_desc-disconnectionpct.json
sub-01_ses-01_label-acuteinfarct_method-snm_space-MNI152NLin6Asym_desc-disconnectionpct.nii.gz
sub-01_ses-01_label-acuteinfarct_space-MNI152NLin2009cAsym_mask.nii.gz

Visualize the output.

In [10]:

import nibabel as nib
from nilearn import plotting

disconnection_path = '/tmp/outputs_snm_run1/sub-01/ses-01/anat/sub-01_ses-01_label-acuteinfarct_method-snm_space-MNI152NLin6Asym_desc-disconnectionpct.nii.gz'
disconnection_img = nib.load(disconnection_path)

disconnection_map = plotting.plot_stat_map(disconnection_img, 
                       radiological=True, 
                       title="Disconnection map from SNM \n(% of streamlines connecting to the lesion)", 
                       draw_cross=False, 
                       colorbar=True,
                       dim=-0.2,
                       cmap='Reds')

import nibabel as nib from nilearn import plotting disconnection_path = '/tmp/outputs_snm_run1/sub-01/ses-01/anat/sub-01_ses-01_label-acuteinfarct_method-snm_space-MNI152NLin6Asym_desc-disconnectionpct.nii.gz' disconnection_img = nib.load(disconnection_path) disconnection_map = plotting.plot_stat_map(disconnection_img, radiological=True, title="Disconnection map from SNM \n(% of streamlines connecting to the lesion)", draw_cross=False, colorbar=True, dim=-0.2, cmap='Reds')

Obtain parcel-level disconnectivity scores and disconnectivity matrix¶

Beyond voxel-level maps, Lacuna covers two further formats of disconnectivity information: region of interest-level disconnectivity scores and a disconnectivity matrix.

These require the --compute-roi-disconnection flag as well as definition of a parcellation atlas via --parcel-atlases.

Run the analysis again with these flags.

In [11]:

!lacuna run snm \
    /tmp/tutorial_bids/ \
    /tmp/outputs_snm_run2/ \
    --connectome-path /tmp/hcp1065_data/hcp1065.tck \
    --participant-label 01 \
    --mask-space MNI152NLin6Asym \
    --compute-disconnectivity-matrix \
    --compute-roi-disconnection \
    --parcel-atlases schaefer2018parcels100networks7

!lacuna run snm \ /tmp/tutorial_bids/ \ /tmp/outputs_snm_run2/ \ --connectome-path /tmp/hcp1065_data/hcp1065.tck \ --participant-label 01 \ --mask-space MNI152NLin6Asym \ --compute-disconnectivity-matrix \ --compute-roi-disconnection \ --parcel-atlases schaefer2018parcels100networks7

Matrix outputs comprise four complementary connectivity representations:

1. Disconnectivity matrix – percentage of streamlines disconnected by the lesion
2. Mask connectivity matrix – number of streamlines traversing the lesion mask
3. Intact connectivity matrix – number of remaining streamlines after excluding those passing through the lesion (i.e., total minus mask connectivity)
4. Full connectome matrix – total number of streamlines irrespective of lesion involvement

In [12]:

!ls /tmp/outputs_snm_run2/sub-01/ses-01/anat/*matrix*

!ls /tmp/outputs_snm_run2/sub-01/ses-01/anat/*matrix*

/tmp/outputs_snm_run2/sub-01/ses-01/anat/sub-01_ses-01_label-acuteinfarct_method-snm_atlas-schaefer2018parcels100networks7_desc-disconnectionpct_connmatrix.json
/tmp/outputs_snm_run2/sub-01/ses-01/anat/sub-01_ses-01_label-acuteinfarct_method-snm_atlas-schaefer2018parcels100networks7_desc-disconnectionpct_connmatrix.tsv
/tmp/outputs_snm_run2/sub-01/ses-01/anat/sub-01_ses-01_label-acuteinfarct_method-snm_atlas-schaefer2018parcels100networks7_desc-fullconnectivitymatrix_connmatrix.json
/tmp/outputs_snm_run2/sub-01/ses-01/anat/sub-01_ses-01_label-acuteinfarct_method-snm_atlas-schaefer2018parcels100networks7_desc-fullconnectivitymatrix_connmatrix.tsv
/tmp/outputs_snm_run2/sub-01/ses-01/anat/sub-01_ses-01_label-acuteinfarct_method-snm_atlas-schaefer2018parcels100networks7_desc-intactconnectivitymatrix_connmatrix.json
/tmp/outputs_snm_run2/sub-01/ses-01/anat/sub-01_ses-01_label-acuteinfarct_method-snm_atlas-schaefer2018parcels100networks7_desc-intactconnectivitymatrix_connmatrix.tsv
/tmp/outputs_snm_run2/sub-01/ses-01/anat/sub-01_ses-01_label-acuteinfarct_method-snm_atlas-schaefer2018parcels100networks7_desc-maskconnectivitymatrix_connmatrix.json
/tmp/outputs_snm_run2/sub-01/ses-01/anat/sub-01_ses-01_label-acuteinfarct_method-snm_atlas-schaefer2018parcels100networks7_desc-maskconnectivitymatrix_connmatrix.tsv
/tmp/outputs_snm_run2/sub-01/ses-01/anat/sub-01_ses-01_label-acuteinfarct_method-snm_atlas-schaefer2018parcels100networks7_desc-matrixstatistics_stats.json

Visualize the disconnectivity of the visual network. Because HCP1065 is mainly consisting of long fiber bundles the sparsity of the matrix is expected.

In [13]:

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

file_path = "/tmp/outputs_snm_run2/sub-01/ses-01/anat/sub-01_ses-01_label-acuteinfarct_method-snm_atlas-schaefer2018parcels100networks7_desc-disconnectionpct_connmatrix.tsv"
mat = pd.read_csv(file_path, sep="\t", header=0, index_col=0)

visual_network_regions = [parc for parc in mat.index if "Vis" in parc]
mat = mat.loc[visual_network_regions, visual_network_regions]
mat = mat.values

plt.figure(figsize=(10, 8))
im = plt.imshow(mat, aspect='auto')
plt.colorbar(im, label="Disconnectivity (%)")

plt.title("Disconnectivity matrix of visual network (Schaefer100x7)")
plt.xlabel("Region")
plt.ylabel("Region")

im.set_clim(0, 100)

plt.tight_layout()
plt.show()

import pandas as pd import numpy as np import matplotlib.pyplot as plt file_path = "/tmp/outputs_snm_run2/sub-01/ses-01/anat/sub-01_ses-01_label-acuteinfarct_method-snm_atlas-schaefer2018parcels100networks7_desc-disconnectionpct_connmatrix.tsv" mat = pd.read_csv(file_path, sep="\t", header=0, index_col=0) visual_network_regions = [parc for parc in mat.index if "Vis" in parc] mat = mat.loc[visual_network_regions, visual_network_regions] mat = mat.values plt.figure(figsize=(10, 8)) im = plt.imshow(mat, aspect='auto') plt.colorbar(im, label="Disconnectivity (%)") plt.title("Disconnectivity matrix of visual network (Schaefer100x7)") plt.xlabel("Region") plt.ylabel("Region") im.set_clim(0, 100) plt.tight_layout() plt.show()

Visualize the parcel-level disconnection (%) on the cortical surface with yabplot

In [16]:

import warnings
import pyvista as pv
import yabplot as yab

warnings.simplefilter("ignore", pv.PyVistaDeprecationWarning)

pv.set_jupyter_backend('static')
pv.global_theme.notebook = True
pv.start_xvfb()

file_path = "/tmp/outputs_snm_run2/sub-01/ses-01/anat/sub-01_ses-01_label-acuteinfarct_method-snm_atlas-schaefer2018parcels100networks7_desc-disconnectionpct_parcelstats.tsv"

df = pd.read_csv(file_path, sep="\t", index_col=0)
values = df.values.flatten()

yab.plot_cortical(
    data=values,
    atlas='schaefer_100',
    views=['left_lateral', 'left_medial', 'right_lateral', 'right_medial'],
    figsize=(800, 300),
    vminmax=(0,10),
    cmap="Reds",
    display_type="static",
)

import warnings import pyvista as pv import yabplot as yab warnings.simplefilter("ignore", pv.PyVistaDeprecationWarning) pv.set_jupyter_backend('static') pv.global_theme.notebook = True pv.start_xvfb() file_path = "/tmp/outputs_snm_run2/sub-01/ses-01/anat/sub-01_ses-01_label-acuteinfarct_method-snm_atlas-schaefer2018parcels100networks7_desc-disconnectionpct_parcelstats.tsv" df = pd.read_csv(file_path, sep="\t", index_col=0) values = df.values.flatten() yab.plot_cortical( data=values, atlas='schaefer_100', views=['left_lateral', 'left_medial', 'right_lateral', 'right_medial'], figsize=(800, 300), vminmax=(0,10), cmap="Reds", display_type="static", )

Run structural network mapping on multiple subjects¶

Lacuna supports processing multiple subjects within a single run. If the --participant-label flag is omitted, the analysis is automatically executed for all subjects detected in the BIDS dataset. With the --verbose we can see which masks are detected.

Execution strategies are optimized internally for each analysis type. For structural network mapping, subjects are processed sequentially. Parallelizing across subjects is intentionally avoided, as concurrent processes would repeatedly access the same tractogram file, leading to race conditions and degraded performance.

Instead, parallelism is implemented at the single-subject level via multithreading, controlled by the --nprocs parameter.

In [ ]:

remove_output

!lacuna run snm \
    /tmp/tutorial_bids/ \
    /tmp/outputs_snm_run3/ \
    --connectome-path /tmp/hcp1065_data/hcp1065.tck \
    --mask-space MNI152NLin6Asym \
    --compute-roi-disconnection \
    --parcel-atlases schaefer2018parcels100networks7 \
    --nprocs 12 \
    --verbose

!lacuna run snm \ /tmp/tutorial_bids/ \ /tmp/outputs_snm_run3/ \ --connectome-path /tmp/hcp1065_data/hcp1065.tck \ --mask-space MNI152NLin6Asym \ --compute-roi-disconnection \ --parcel-atlases schaefer2018parcels100networks7 \ --nprocs 12 \ --verbose

List the parcel-level outputs for multiple subjects

In [18]:

!ls /tmp/outputs_snm_run3/sub-*/*/*/sub-*_ses-01_label-acuteinfarct_atlas-schaefer2018parcels100networks7_method-snm_desc-roidisconnection_parcelstats.tsv

!ls /tmp/outputs_snm_run3/sub-*/*/*/sub-*_ses-01_label-acuteinfarct_atlas-schaefer2018parcels100networks7_method-snm_desc-roidisconnection_parcelstats.tsv

zsh:1: no matches found: /tmp/outputs_snm_run3/sub-*/*/*/sub-*_ses-01_label-acuteinfarct_atlas-schaefer2018parcels100networks7_method-snm_desc-roidisconnection_parcelstats.tsv

Collect parcel-level outputs into a single tsv

In [65]:

!lacuna collect \
    /tmp/outputs_snm_run3/ \
    --pattern "*roidisconnection*" \
    --output-dir /tmp/group_outputs/ \
    --verbose

!lacuna collect \ /tmp/outputs_snm_run3/ \ --pattern "*roidisconnection*" \ --output-dir /tmp/group_outputs/ \ --verbose

2026-03-19 17:57:42 - lacuna.cli.main - INFO - Running collect (group-level aggregation)
2026-03-19 17:57:42 - lacuna.cli.main - INFO - Scanning derivatives directory: /tmp/outputs_snm_run3
2026-03-19 17:57:42 - lacuna.cli.main - INFO - Output directory: /tmp/group_outputs
2026-03-19 17:57:42 - lacuna.cli.main - INFO - Pattern: *roidisconnection*_parcelstats.tsv
2026-03-19 17:57:42 - lacuna.cli.main - INFO - Found 3 file(s) across 3 subject(s)
2026-03-19 17:57:42 - lacuna.cli.main - INFO - Created 1 group-level TSV file(s):
2026-03-19 17:57:42 - lacuna.cli.main - INFO -   - group_atlas-schaefer2018parcels100networks7_method-snm_desc-roidisconnection_parcelstats.tsv

In [67]:

!ls /tmp/group_outputs/*

!ls /tmp/group_outputs/*

/tmp/group_outputs/group_atlas-schaefer2018parcels100networks7_method-snm_desc-roidisconnection_parcelstats.json
/tmp/group_outputs/group_atlas-schaefer2018parcels100networks7_method-snm_desc-roidisconnection_parcelstats.tsv

In [69]:

parcel_df = pd.read_csv("/tmp/group_outputs/group_atlas-schaefer2018parcels100networks7_method-snm_desc-roidisconnection_parcelstats.tsv", sep="\t", index_col=0)
parcel_df

parcel_df = pd.read_csv("/tmp/group_outputs/group_atlas-schaefer2018parcels100networks7_method-snm_desc-roidisconnection_parcelstats.tsv", sep="\t", index_col=0) parcel_df

Out[69]:

	session_id	label	7Networks_LH_Vis_1	7Networks_LH_Vis_2	7Networks_LH_Vis_3	7Networks_LH_Vis_4	7Networks_LH_Vis_5	7Networks_LH_Vis_6	7Networks_LH_Vis_7	7Networks_LH_Vis_8	...	7Networks_RH_Default_Temp_1	7Networks_RH_Default_Temp_2	7Networks_RH_Default_Temp_3	7Networks_RH_Default_PFCv_1	7Networks_RH_Default_PFCv_2	7Networks_RH_Default_PFCdPFCm_1	7Networks_RH_Default_PFCdPFCm_2	7Networks_RH_Default_PFCdPFCm_3	7Networks_RH_Default_pCunPCC_1	7Networks_RH_Default_pCunPCC_2
participant_id
1	1	acuteinfarct	0.568720	0.697245	5.758313	1.385571	4.475859	51.615445	0.0	0.969642	...	0.0	7.530272	0.000000	5.670816	0.371594	0.483325	1.913937	0.0	11.744966	0.029343
2	1	acuteinfarct	0.000000	1.175747	0.000000	0.000000	1.423070	0.000000	0.0	0.000000	...	0.0	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.0	0.000000	0.000000
3	1	acuteinfarct	0.189573	0.000000	0.000000	0.000000	0.000000	0.000000	0.0	0.000000	...	0.0	0.012614	2.232346	0.000000	0.000000	0.000000	0.000000	0.0	0.000000	0.000000

3 rows × 102 columns