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 13877248 of ? bytes.
[fetch_single_file] Downloaded 27525120 of ? bytes.
[fetch_single_file] Downloaded 41426944 of ? bytes.
[fetch_single_file] Downloaded 55312384 of ? bytes.
[fetch_single_file] Downloaded 58982400 of ? bytes.
[fetch_single_file] Downloaded 63184896 of ? bytes.
[fetch_single_file] Downloaded 67117056 of ? bytes.
[fetch_single_file] Downloaded 74711040 of ? bytes.
[fetch_single_file] Downloaded 94642176 of ? bytes.
[fetch_single_file] Downloaded 105644032 of ? bytes.
[fetch_single_file] Downloaded 108724224 of ? bytes.
[fetch_single_file] Downloaded 111673344 of ? bytes.
[fetch_single_file] Downloaded 114556928 of ? bytes.
[fetch_single_file] Downloaded 118489088 of ? bytes.
[fetch_single_file] Downloaded 133439488 of ? bytes.
[fetch_single_file] Downloaded 136323072 of ? bytes.
[fetch_single_file] Downloaded 139206656 of ? bytes.
[fetch_single_file] Downloaded 142868480 of ? bytes.
[fetch_single_file] Downloaded 147849216 of ? bytes.
[fetch_single_file] Downloaded 150732800 of ? bytes.
[fetch_single_file] Downloaded 154632192 of ? bytes.
[fetch_single_file] Downloaded 158867456 of ? bytes.
[fetch_single_file] Downloaded 162627584 of ? bytes.
[fetch_single_file] Downloaded 168566784 of ? bytes.
[fetch_single_file] Downloaded 171687936 of ? bytes.
[fetch_single_file] Downloaded 174596096 of ? bytes.
[fetch_single_file] Downloaded 177733632 of ? bytes.
[fetch_single_file] Downloaded 194576384 of ? bytes.
[fetch_single_file] Downloaded 204218368 of ? bytes.
[fetch_single_file] Downloaded 207101952 of ? bytes.
[fetch_single_file] Downloaded 209977344 of ? bytes.
[fetch_single_file] Downloaded 213385216 of ? bytes.
[fetch_single_file] Downloaded 217579520 of ? bytes.
[fetch_single_file] Downloaded 220471296 of ? bytes.
[fetch_single_file] Downloaded 227278848 of ? bytes.
[fetch_single_file] Downloaded 230408192 of ? bytes.
[fetch_single_file] Downloaded 233308160 of ? bytes.
[fetch_single_file] Downloaded 236453888 of ? bytes.
[fetch_single_file] Downloaded 239607808 of ? bytes.
[fetch_single_file] Downloaded 244064256 of ? bytes.
[fetch_single_file] Downloaded 250716160 of ? bytes.
[fetch_single_file] Downloaded 275521536 of ? bytes.
[fetch_single_file] Downloaded 300941312 of ? bytes.
[fetch_single_file] Downloaded 321658880 of ? bytes.
[fetch_single_file] Downloaded 343670784 of ? bytes.
[fetch_single_file] Downloaded 364339200 of ? bytes.
[fetch_single_file] Downloaded 385089536 of ? bytes.
[fetch_single_file] Downloaded 406847488 of ? bytes.
[fetch_single_file] Downloaded 425877504 of ? bytes.
[fetch_single_file] Downloaded 447209472 of ? bytes.
[fetch_single_file] Downloaded 468197376 of ? bytes.
[fetch_single_file] Downloaded 487333888 of ? bytes.
[fetch_single_file] Downloaded 509329408 of ? bytes.
[fetch_single_file] Downloaded 531111936 of ? bytes.
[fetch_single_file] Downloaded 552345600 of ? bytes.
[fetch_single_file] Downloaded 574357504 of ? bytes.
[fetch_single_file] Downloaded 596385792 of ? bytes.
[fetch_single_file]  ...done. (60 seconds, 0 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 0x7f94157a5190>

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. (3 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 0x7F942B378620>

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

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 0x7f94252cf8c0>

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. (2 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 0x7f9425322ba0>

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.