Multimodal comparison of two cortical brain areas

The tutorial shows how maps and mutimodal regional measurements are obtained using siibra-python for language area IFG 44 and primary visual region hOc1, as defined in the Julich-Brain cytoarchitectonic atlas.

import matplotlib.pyplot as plt
from nilearn import plotting
import seaborn as sns
import siibra

assert siibra.__version__ >= "1.0.1"

sns.set_style("dark")

Instantiate parcellation and reference space from the human brain atlas

jubrain = siibra.parcellations.JULICH_BRAIN_CYTOARCHITECTONIC_ATLAS_V3_0_3
pmaps = jubrain.get_map(space="mni152", maptype="statistical")
print(jubrain.publications[0]['citation'])

Amunts, K., Mohlberg, H., Bludau, S., Zilles, K. (2020). Julich-Brain – A 3D probabilistic atlas of human brain’s cytoarchitecture. Science 369, 988-992

Obtain definitions and probabilistic maps in MNI asymmetric space of area IFG 44 and primary visual region hOc1 as defined in the Julich-Brain cytoarchitectonic atlas

specs = ["ifg 44 left", "hoc1 left"]
regions = [jubrain.get_region(spec) for spec in specs]
for r in regions:
    plotting.plot_glass_brain(
        pmaps.fetch(region=r),
        cmap="viridis",
        draw_cross=False,
        colorbar=False,
        annotate=False,
        symmetric_cbar=True,
        title=r.name,
        threshold=0
    )

• 
• 

For each of the two brain areas, query for layer-specific distributions of cell bodies.

fig, axs = plt.subplots(1, len(regions), sharey=True, figsize=(8, 2.7))
for i, region in enumerate(regions):
    layerwise_cellbody_densities = siibra.features.get(region, "layerwise cell density")
    layerwise_cellbody_densities[0].plot(ax=axs[i], textwrap=25)
    print(layerwise_cellbody_densities[0].urls)
    axs[i].set_ylim(25, 115)

['https://doi.org/10.25493/THE9-VQU']
['https://doi.org/10.25493/EVPC-MVH']

Next, retrieve average densities of a selection of monogenetic neurotransmitter receptors.

receptors = ["M1", "M2", "M3", "5-HT1A", "5-HT2", "D1"]
fig, axs = plt.subplots(1, len(regions), sharey=True, figsize=(8, 3.5))
for i, region in enumerate(regions):
    receptor_fingerprints = siibra.features.get(region, "receptor density fingerprint")
    receptor_fingerprints[0].plot(receptors=receptors, ax=axs[i])
    print(receptor_fingerprints[0].urls)
    axs[i].set_ylim(0, 1000)
    axs[i].title.set_size(12)
    transmitters = [
        siibra.vocabularies.RECEPTOR_SYMBOLS[r]['neurotransmitter']['name']
        for r in receptors
    ]
    axs[i].set_xticklabels([f"{r}\n({n})" for r, n in zip(receptors, transmitters)])

['https://doi.org/10.25493/P82M-PVM', 'https://doi.org/10.25493/P82M-PVM']
['https://doi.org/10.25493/P8SD-JMH', 'https://doi.org/10.25493/P8SD-JMH']

Now, for further insight, query for expressions of a selection of genes coding for these receptors. :ignore:

genes = ["CHRM1", "CHRM2", "CHRM3", "HTR1A", "HTR2A", "DRD1"]
fig, axs = plt.subplots(1, len(regions), sharey=True, figsize=(7, 3))
for i, region in enumerate(regions):
    gene_expressions = siibra.features.get(region, "gene expressions", gene=genes)
    assert len(gene_expressions) == 1
    gene_expressions[0].plot(ax=axs[i])
    axs[i].title.set_size(11)
    print(gene_expressions[0].urls)

Traceback (most recent call last):
  File "/home/docs/checkouts/readthedocs.org/user_builds/siibra-python/checkouts/doc_test_ng_resolve/examples/tutorials/2025-paper-fig4.py", line 92, in <module>
    gene_expressions = siibra.features.get(region, "gene expressions", gene=genes)
                       ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/docs/checkouts/readthedocs.org/user_builds/siibra-python/checkouts/doc_test_ng_resolve/siibra/features/feature.py", line 551, in _match
    return cls._match(concept, ftype_candidates, **kwargs)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/docs/checkouts/readthedocs.org/user_builds/siibra-python/checkouts/doc_test_ng_resolve/siibra/features/feature.py", line 533, in _match
    sum((
  File "/home/docs/checkouts/readthedocs.org/user_builds/siibra-python/checkouts/doc_test_ng_resolve/siibra/features/feature.py", line 534, in <genexpr>
    cls._match(concept, t, **kwargs) for t in ftypes
    ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/docs/checkouts/readthedocs.org/user_builds/siibra-python/checkouts/doc_test_ng_resolve/siibra/features/feature.py", line 593, in _match
    live_instances = feature_type._livequery(concept, **kwargs)
                     ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/docs/checkouts/readthedocs.org/user_builds/siibra-python/checkouts/doc_test_ng_resolve/siibra/features/feature.py", line 494, in _livequery
    features = q.query(concept)
               ^^^^^^^^^^^^^^^^
  File "/home/docs/checkouts/readthedocs.org/user_builds/siibra-python/checkouts/doc_test_ng_resolve/siibra/livequeries/allen.py", line 144, in query
    raise InvalidAllenAPIResponseException(
siibra.livequeries.allen.InvalidAllenAPIResponseException: The service "web API of the Allen Brain Atlas for the human microarray expression" is not available at the moment, therefore siibra is not able to fetch gene expression features. This is a known issue which we are investigating: https://github.com/FZJ-INM1-BDA/siibra-python/issues/636.

For each of the two brain areas, collect functional connectivity profiles referring to temporal correlation of fMRI timeseries of several hundred subjects from the Human Connectome Project. We show the strongest connections per brain area for the average connectivity patterns

fts = siibra.features.get(jubrain, "functional connectivity")
for cf in fts:
    if cf.cohort != "HCP":
        continue
    if cf.paradigm == "Resting state (EmpCorrFC concatenated)":
        conn = cf
        break
print(conn.urls)


# plot both average connectivity profiles to target regions
def shorten_name(n):
    return (
        n.replace("Area ", "")
        .replace(" (GapMap)", "")
        .replace("left", "L")
        .replace("right", "R")
    )


plotkwargs = {
    "kind": "bar",
    "width": 0.85,
    "logscale": True,
    "xlabel": "",
    "ylabel": "temporal correlation",
    "rot": 90,
    "ylim": (0.3, 1.1),
    "capsize": 2,
}
fig, axs = plt.subplots(1, len(regions), sharey=True, figsize=(7, 4.5))
for i, region in enumerate(regions):
    plotkwargs["ax"] = axs[i]
    conn.plot(region, max_rows=17, **plotkwargs)
    axs[i].xaxis.set_ticklabels([shorten_name(t.get_text()) for t in axs[i].xaxis.get_majorticklabels()])
    axs[i].set_title(region.name.replace("Area ", ""), fontsize=8)
    axs[i].grid(True, 'minor')
plt.suptitle("Functional Connectivity")

For both brain areas, sample representative cortical image patches at 1µm resolution that were automatically extracted from scans of BigBrain sections. The image patches display clearly the differences in laminar structure of the two regions

selected_sections = [4950, 1345]
fig, axs = plt.subplots(1, len(regions), sharey=True)
for r, sn, ax in zip(regions, selected_sections, axs):
    # find 1 micron sections intersecting the region
    pmap = pmaps.get_volume(r)
    all_patches = siibra.features.get(pmap, "BigBrain1MicronPatch")
    patches = [p for p in all_patches if p.bigbrain_section == sn]

    # select the first patch and access the underlying image data
    patchimg = patches[0].fetch()
    patchdata = patchimg.get_fdata().squeeze()

    # plot the pure image array
    ax.imshow(patchdata, cmap="gray")
    ax.axis("off")
    ax.set_title(r.name)

Estimated memory usage: 1747 MB