Network assortativity and distance-preserving surrogates

This example demonstrates how to calculate network assortativity and generate distance-preserving surrogate connectomes for null hypothesis testing. We’ll use structural connectivity data and myelin maps from Bazinet et al. (2023) to show how brain network topology relates to spatial patterns of cortical features.

First, let’s fetch the structural connectivity and myelin data from Bazinet et al. (2023). This dataset contains a 400-parcel structural connectome and corresponding T1w/T2w myelin maps:

import importlib
import numpy as np
import pickle
from netneurotools.datasets import fetch_bazinet_assortativity

bazinet_path = fetch_bazinet_assortativity()

with open(f'{bazinet_path}/data/human_SC_s400.pickle', 'rb') as file:
    human_data = pickle.load(file)

myelin = human_data['t1t2']
SC = human_data['adj']

print(f'Structural connectivity shape: {SC.shape}')
print(f'Myelin data shape: {myelin.shape}')

Please cite the following papers if you are using this function:
  [primary]:
    Vincent Bazinet, Justine Y Hansen, Reinder Vos de Wael, Boris C Bernhardt, Martijn P van den Heuvel, and Bratislav Misic. Assortative mixing in micro-architecturally annotated brain connectomes. Nature Communications, 14(1):2850, 2023.
[fetch_single_file] Downloading data from
https://github.com/netneurolab/bazinet_assortativity/archive/baab96d173b8c5a2656
67a4259cf228713dabcc7.tar.gz ...
[fetch_single_file] Downloaded 13639680 of ? bytes.
[fetch_single_file] Downloaded 27533312 of ? bytes.
[fetch_single_file] Downloaded 40640512 of ? bytes.
[fetch_single_file] Downloaded 52961280 of ? bytes.
[fetch_single_file] Downloaded 57933824 of ? bytes.
[fetch_single_file] Downloaded 61865984 of ? bytes.
[fetch_single_file] Downloaded 65798144 of ? bytes.
[fetch_single_file] Downloaded 69730304 of ? bytes.
[fetch_single_file] Downloaded 86597632 of ? bytes.
[fetch_single_file] Downloaded 102244352 of ? bytes.
[fetch_single_file] Downloaded 106700800 of ? bytes.
[fetch_single_file] Downloaded 109314048 of ? bytes.
[fetch_single_file] Downloaded 112459776 of ? bytes.
[fetch_single_file] Downloaded 115613696 of ? bytes.
[fetch_single_file] Downloaded 124256256 of ? bytes.
[fetch_single_file] Downloaded 134225920 of ? bytes.
[fetch_single_file] Downloaded 137101312 of ? bytes.
[fetch_single_file] Downloaded 139984896 of ? bytes.
[fetch_single_file] Downloaded 143663104 of ? bytes.
[fetch_single_file] Downloaded 148381696 of ? bytes.
[fetch_single_file] Downloaded 151257088 of ? bytes.
[fetch_single_file] Downloaded 156246016 of ? bytes.
[fetch_single_file] Downloaded 159391744 of ? bytes.
[fetch_single_file] Downloaded 165937152 of ? bytes.
[fetch_single_file] Downloaded 168828928 of ? bytes.
[fetch_single_file] Downloaded 171712512 of ? bytes.
[fetch_single_file] Downloaded 174620672 of ? bytes.
[fetch_single_file] Downloaded 177717248 of ? bytes.
[fetch_single_file] Downloaded 191627264 of ? bytes.
[fetch_single_file] Downloaded 203948032 of ? bytes.
[fetch_single_file] Downloaded 206577664 of ? bytes.
[fetch_single_file] Downloaded 209461248 of ? bytes.
[fetch_single_file] Downloaded 212869120 of ? bytes.
[fetch_single_file] Downloaded 216866816 of ? bytes.
[fetch_single_file] Downloaded 219684864 of ? bytes.
[fetch_single_file] Downloaded 222568448 of ? bytes.
[fetch_single_file] Downloaded 229384192 of ? bytes.
[fetch_single_file] Downloaded 232267776 of ? bytes.
[fetch_single_file] Downloaded 235413504 of ? bytes.
[fetch_single_file] Downloaded 238559232 of ? bytes.
[fetch_single_file] Downloaded 242745344 of ? bytes.
[fetch_single_file] Downloaded 246161408 of ? bytes.
[fetch_single_file] Downloaded 266346496 of ? bytes.
[fetch_single_file] Downloaded 289406976 of ? bytes.
[fetch_single_file] Downloaded 313532416 of ? bytes.
[fetch_single_file] Downloaded 336068608 of ? bytes.
[fetch_single_file] Downloaded 355196928 of ? bytes.
[fetch_single_file] Downloaded 379592704 of ? bytes.
[fetch_single_file] Downloaded 401088512 of ? bytes.
[fetch_single_file] Downloaded 421527552 of ? bytes.
[fetch_single_file] Downloaded 443514880 of ? bytes.
[fetch_single_file] Downloaded 462692352 of ? bytes.
[fetch_single_file] Downloaded 486014976 of ? bytes.
[fetch_single_file] Downloaded 507256832 of ? bytes.
[fetch_single_file] Downloaded 529514496 of ? bytes.
[fetch_single_file] Downloaded 548749312 of ? bytes.
[fetch_single_file] Downloaded 568049664 of ? bytes.
[fetch_single_file] Downloaded 592977920 of ? bytes.
[fetch_single_file]  ...done. (61 seconds, 1 min)

Downloaded ds-bazinet_assortativity to /home/runner/nnt-data
/home/runner/work/netneurotools/netneurotools/examples/plot_assortativity.py:27: DeprecationWarning: numpy.core.numeric is deprecated and has been renamed to numpy._core.numeric. The numpy._core namespace contains private NumPy internals and its use is discouraged, as NumPy internals can change without warning in any release. In practice, most real-world usage of numpy.core is to access functionality in the public NumPy API. If that is the case, use the public NumPy API. If not, you are using NumPy internals. If you would still like to access an internal attribute, use numpy._core.numeric._frombuffer.
  human_data = pickle.load(file)
Structural connectivity shape: (400, 400)
Myelin data shape: (400,)

Let’s visualize the structural connectivity matrix. This shows the weighted connections between all pairs of brain regions:

import matplotlib.pyplot as plt

fig, ax = plt.subplots(figsize=(6, 6))
im = ax.imshow(SC, cmap='Greys', vmin=SC[SC > 0].min(), vmax=SC.max())
ax.set(xlabel='ROI', ylabel='ROI', title='Structural Connectivity')
fig.colorbar(im, ax=ax)

<matplotlib.colorbar.Colorbar object at 0x7f3afd76a1e0>

To visualize the myelin data on the cortical surface, we need to project the parcellated data to the vertex level. We’ll use the Schaefer 400-parcel atlas to do this mapping:

from netneurotools.datasets import fetch_schaefer2018
from netneurotools.interface import parcels_to_vertices

parc = fetch_schaefer2018('fsaverage')['400Parcels7Networks']
myelin_hemi = (myelin[:200], myelin[200:])
myelin_vertex, _, _ = parcels_to_vertices(myelin_hemi, parc, hemi='both')

Please cite the following papers if you are using this function:
  [primary]:
    Alexander Schaefer, Ru Kong, Evan M Gordon, Timothy O Laumann, Xi-Nian Zuo, Avram J Holmes, Simon B Eickhoff, and BT Thomas Yeo. Local-global parcellation of the human cerebral cortex from intrinsic functional connectivity mri. Cerebral cortex, 28(9):3095–3114, 2018.
[fetch_single_file] Downloading data from
https://files.osf.io/v1/resources/udpv8/providers/osfstorage/673273c54b79ea9a506
2f4b3 ...
[fetch_single_file]  ...done. (1 seconds, 0 min)

Downloaded atl-schaefer2018 to /home/runner/nnt-data

Now we can plot the myelin data on the inflated cortical surface using PyVista:

from netneurotools.plotting import pv_plot_surface

# This is necessary for headless rendering when building the sphinx gallery.
# If you are running this code locally, you might not need this line.
import os
os.environ["VTK_DEFAULT_OPENGL_WINDOW"] = "vtkOSOpenGLRenderWindow"

pv_plot_surface(
    myelin_vertex,
    'fsaverage',
    'inflated',
    hemi='both',
    cmap="Spectral_r",
    clim=[np.nanmin(myelin_vertex), np.nanmax(myelin_vertex)],
    cbar_title='T1w/T2w ratio',
    lighting_style='plastic',
    jupyter_backend="static"
)

<PIL.Image.Image image mode=RGB size=1400x1000 at 0x7F3AFCA23980>

<pyvista.plotting.plotter.Plotter object at 0x7f3afca22150>

Network assortativity measures the tendency for nodes with similar features to be connected. Let’s calculate the assortativity of the myelin distribution with respect to the structural connectivity network:

from netneurotools.metrics import assortativity_und

assort = assortativity_und(myelin, SC, use_numba=False)
print(f'Network assortativity: {assort:.3f}')

Network assortativity: 0.463

We can visualize this by creating a scatterplot showing the myelin values of connected node pairs, with point colors indicating connection strength:

SC_E = SC > 0
Xi = np.repeat(myelin[np.newaxis, :], 400, axis=0)
Xj = np.repeat(myelin[:, np.newaxis], 400, axis=1)
sort_ids = np.argsort(SC[SC_E])

fig, ax = plt.subplots(figsize=(6, 6))
scatter = ax.scatter(
    Xi[SC_E][sort_ids],
    Xj[SC_E][sort_ids],
    c=SC[SC_E][sort_ids],
    cmap='Greys',
    s=3,
    rasterized=True
)
ax.set(xlabel='Myelin (node i)', ylabel='Myelin (node j)')
ax.set_title(f'Assortativity = {assort:.3f}')
fig.colorbar(scatter, ax=ax, label='Connection strength')

<matplotlib.colorbar.Colorbar object at 0x7f3b0c3bc6e0>

To test whether this assortativity is significant, we need null models that preserve key topological properties. We’ll generate distance-preserving surrogate networks that maintain the degree distribution and spatial embedding of connections.

First, let’s calculate the Euclidean distances between parcel centroids:

from scipy.spatial.distance import cdist

get_parcel_centroids = importlib.import_module(
    'neuromaps.nulls.spins'
).get_parcel_centroids
fetch_atlas = importlib.import_module('neuromaps.datasets').fetch_atlas
neuromaps_images = importlib.import_module('neuromaps.images')
relabel_gifti = neuromaps_images.relabel_gifti
annot_to_gifti = neuromaps_images.annot_to_gifti

parc_gii = [relabel_gifti(annot_to_gifti(parc_hemi))[0] for parc_hemi in parc]
fsaverage_atlas = fetch_atlas('fsaverage', '164k')['pial']
parc_centroids, _ = get_parcel_centroids(fsaverage_atlas, parc_gii)
parc_dist = cdist(parc_centroids, parc_centroids)

print(f'Distance matrix shape: {parc_dist.shape}')

[fetch_files] Downloading data from
https://files.osf.io/v1/resources/4mw3a/providers/osfstorage/60b684b2cb2a5e01e96
8c918 ...
[fetch_files]  ...done. (1 seconds, 0 min)

[fetch_files] Extracting data from
/home/runner/neuromaps-data/030343e2c52d5517d71ab63dbec6a5ae/fsaverage164k.tar.g
z...
[fetch_files] .. done.

Distance matrix shape: (400, 400)

This is also a good point to make the workflow boundaries explicit. The tractography-derived connectome, parcel geometry, and centroid estimation all come from upstream resources, while netneurotools provides the analysis interface points used here: annotation-based assortativity, parcel/vertex conversion, and distance-preserving null models.

fig, ax = plt.subplots(figsize=(11, 2.8))
ax.axis('off')

boxes = [
    (0.02, 0.55, 0.20, 0.26, "Upstream\ntractography +\nparcellation", "#e9ecef"),
    (0.27, 0.55, 0.18, 0.26, "Surface geometry\n+ centroids\n(neuromaps)", "#e9ecef"),
    (0.50, 0.55, 0.18, 0.26, "netneurotools\ndatasets / interface", "#d8f3dc"),
    (
        0.73,
        0.55,
        0.22,
        0.26,
        "netneurotools\nassortativity +\nnull rewiring",
        "#d8f3dc",
    ),
    (0.73, 0.14, 0.22, 0.20, "Empirical vs. null\ninterpretation", "#fff3bf"),
]

for x0, y0, width, height, label, facecolor in boxes:
    rect = plt.Rectangle((x0, y0), width, height, facecolor=facecolor,
                         edgecolor='black', linewidth=1.0)
    ax.add_patch(rect)
    ax.text(x0 + width / 2, y0 + height / 2, label,
            ha='center', va='center', fontsize=10)

arrowprops = dict(arrowstyle='->', lw=1.5, color='black')
ax.annotate('', xy=(0.27, 0.68), xytext=(0.22, 0.68), arrowprops=arrowprops)
ax.annotate('', xy=(0.50, 0.68), xytext=(0.45, 0.68), arrowprops=arrowprops)
ax.annotate('', xy=(0.73, 0.68), xytext=(0.68, 0.68), arrowprops=arrowprops)
ax.annotate('', xy=(0.84, 0.34), xytext=(0.84, 0.55), arrowprops=arrowprops)

ax.set_title('Workflow context for this example', fontsize=11)

Text(0.5, 1.0, 'Workflow context for this example')

Now we’ll generate distance-preserving surrogate connectomes. This process rewires the network while preserving both the degree distribution and the relationship between connection probability and distance.

Note: Generating many surrogates can be time-consuming. For this example, we’ll create just a few surrogates:

from netneurotools.networks import match_length_degree_distribution

n_surrogates = 10
surr_all = np.zeros((n_surrogates, 400, 400))

for i in range(n_surrogates):
    surr_all[i] = match_length_degree_distribution(SC, parc_dist)[1]

print(f'Generated {n_surrogates} surrogate networks')

Generated 10 surrogate networks

Let’s visualize a few of these surrogate connectomes to see how they compare to the empirical network:

fig, axes = plt.subplots(1, 3, figsize=(15, 5))
for i, ax in enumerate(axes):
    im = ax.imshow(
        surr_all[i],
        cmap='Greys',
        vmin=surr_all[i][surr_all[i] > 0].min(),
        vmax=surr_all[i].max()
    )
    ax.set_title(f'Surrogate {i + 1}')
    ax.set(xlabel='ROI', ylabel='ROI')
fig.colorbar(im, ax=axes.ravel().tolist(), label='Connection strength')

<matplotlib.colorbar.Colorbar object at 0x7f3b0c41caa0>

Now we’ll calculate the assortativity for each surrogate network to create a null distribution:

assort_surr = np.zeros(n_surrogates)
for i in range(n_surrogates):
    assort_surr[i] = assortativity_und(myelin, surr_all[i])

print(f'Mean surrogate assortativity: {np.mean(assort_surr):.3f}')
print(f'Std surrogate assortativity: {np.std(assort_surr):.3f}')

Mean surrogate assortativity: 0.401
Std surrogate assortativity: 0.010

Finally, we can compare the empirical assortativity to the null distribution using a boxplot. This visualization shows whether the observed assortativity is significantly different from what we’d expect by chance:

fig, ax = plt.subplots(figsize=(4, 6))

flierprops = dict(
    marker='+',
    markerfacecolor='lightgray',
    markeredgecolor='lightgray'
)

bplot = ax.boxplot(
    assort_surr,
    widths=0.3,
    patch_artist=True,
    showfliers=True,
    showcaps=False,
    flierprops=flierprops
)

for box in bplot['boxes']:
    box.set(facecolor='lightgray', edgecolor='black')
for median in bplot['medians']:
    median.set(color='black', linewidth=2)

ax.scatter(1, assort, color='red', s=100, zorder=3, label='Empirical')
ax.set_xticks([1])
ax.set_xticklabels(['Surrogates'])
ax.set_ylabel('Assortativity')
ax.set_title('Empirical vs. Null Distribution')
ax.legend()
ax.set_xlim(0.7, 1.3)

plt.tight_layout()

The red dot shows the empirical assortativity, while the boxplot shows the distribution of assortativity values from the distance-preserving surrogate networks. If the empirical value falls outside the null distribution, this suggests that the observed pattern of myelin assortativity is not explained by the spatial constraints and topology of the connectome alone.

The limitation is equally important: geometry-aware nulls still depend on the quality of the upstream surface / centroid information, and generating large surrogate ensembles can be computationally expensive.