Synaptic clustering (Notebook)

In this tutorial, you will learn how to visualize synapses on neuronal morphology and perform clustering analysis using a custom distance metric (specifically, the path_distance between synaptic locations).

Setup

!pip install dendrotweaks

If you are using Google Colab, you might also need to restart the session as the installation downgraded some packages (numpy). You can do it manually or programmatically as shown below:

# import os
# os.kill(os.getpid(), 9)

Let’s begin by importing the standard libraries and the dendrotweaks library.

import numpy as np
import matplotlib.pyplot as plt
import os

import dendrotweaks as dd
dd.__version__

'0.5.1'

dd.apply_dark_theme()

os.makedirs('examples', exist_ok=True)
if not os.listdir('examples/'):
    print("Downloading example data...")
    dd.download_example_data('examples')

Load morphology

path_to_model = os.path.join('examples', 'Poirazi_2003')
print(f'Path to model: {path_to_model}')

Path to model: examples/Poirazi_2003

model = dd.Model(path_to_model)

model.list_morphologies()

['original', 'main']

model.load_morphology(file_name='main')

Sorted PointTree(root=Point(idx=0), num_nodes=5077).
Extended 178 nodes.
Sorted PointTree(root=Point(idx=0), num_nodes=5255).

model.domains

{'axon': <Domain(axon, 2, #F0E442, 13 sections)>,
 'basal': <Domain(basal, 31, #31A354, 48 sections)>,
 'oblique': <Domain(oblique, 43, #8C564B, 72 sections)>,
 'soma': <Domain(soma, 1, #E69F00, 1 sections)>,
 'trunk': <Domain(trunk, 41, #56B4E9, 25 sections)>,
 'tuft': <Domain(tuft, 42, #A55194, 22 sections)>}

fig, ax = plt.subplots(figsize=(5, 5))
model.sec_tree.plot(ax, 
                    show_points=False,
                    show_lines=True,
                    annotate=False)

We also need to add the default mechanisms (including AMPA.mod) to be able to place synapses on our morphology.

model.add_default_mechanisms(recompile=False)

Synapses

Next, we will randomly allocate 100 AMPA synapses on the sections of the tuft domain.

model.remove_all_stimuli()
model.remove_all_recordings()

soma = model.sec_tree[0]
model.add_recording(soma, loc=0.5)

segments = model.get_segments(['tuft'])
model.add_population(
    'excitatory',
    segments, 
    N=100, 
    syn_type='AMPA'
)

Typically, in model.populations["excitatory"].synapses, we store a mapping from segments to synapses that belong to those segments.

However, we can access all the synapses in this population as a list using the flat_synapses property of the Population class.

synapses = model.populations['excitatory'].flat_synapses

Two important attributes of a Synapse that we will further need are sec and loc, i.e., the section to which the synapse belongs and the relative location of the synapse along that section, respectively (0 - start, 1 - end). These two attributes unambiguously define the position of a synapse.

synapses[0].sec, synapses[0].loc

(NeuronSection(idx=159), 0.026000000000000002)

Distance to soma

We can use syn.sec and syn.loc, for example, to find the path distance from a synapse to the soma.

distances_to_soma = [
    syn.sec.path_distance(relative_position=syn.loc)
     for syn in synapses
    ]

fig, ax = plt.subplots(1, 1, figsize=(8, 2))
ax.hist(distances_to_soma, bins=20, color='C0', edgecolor='k', alpha=0.7)
ax.set_xlabel('Distance to soma (µm)')
ax.set_ylabel('Number of synapses')
ax.set_title('Distribution of synapse locations along the dendritic tree')
ax.set_xlim(0, max(distances_to_soma) * 1.1)
plt.show()

Synapse coordinates

Consider now that we want to find the 3D coordinates for a synapse that correspond to its more abstract coordinates defined by sec and loc attributes. We will utilize a general-purpose get_location_coordinates method of the Section class.

For example, let’s calculate the x, y, z coordinates in the middle of a section.

x, y, z = model.sec_tree[50].get_location_coordinates(loc=0.5)

To get the 3D coordinates of each synapse in our previously obtained list of synapses, we will use the syn.xyz() method that internally relies on the get_location_coordinates method.

syn_coords = [syn.xyz() for syn in synapses]

We can now use these coordinates to visualize the synapses on the section tree plot.

fig, ax = plt.subplots(1, 1, figsize=(8, 6))
model.sec_tree.plot(ax, show_points=False, show_lines=True, annotate=False)
xs = [coord[0] for coord in syn_coords]
ys = [coord[1] for coord in syn_coords]
ax.scatter(xs, ys, color='r', edgecolor='k', alpha=0.7, label='Synapse locations', s=10)

<matplotlib.collections.PathCollection at 0x70629a5a6bd0>

Cluster analysis

Next, we will see how to use an external clustering algorithm, such as DBSCAN (Density-Based Spatial Clustering for Applications with Noise), to split synapses into clusters.

A clustering algorithm essentially relies on calculating distances between data points. Consider the following toy example with some data points on a 2D plane. A distance used for assigning cluster labels here is simply the Euclidean distance between two points.

#!pip install scikit-learn

from sklearn.datasets import make_blobs
from sklearn.cluster import KMeans
X, _ = make_blobs(n_samples=100, centers=3, n_features=2, random_state=42)
kmeans = KMeans(n_clusters=3, random_state=0).fit(X)
point1 = X[22]
point2 = X[13]

fig, ax = plt.subplots(figsize=(10, 5))
ax.scatter(X[:, 0], X[:, 1], c=kmeans.labels_, cmap='tab10', edgecolor='k', alpha=0.7)
ax.plot([point1[0], point2[0]], [point1[1], point2[1]], 'm--', label='Distance between points')
ax.legend()
ax.set_title('Example 2D Dataset for Clustering')
plt.show()

However, in our case, these distances are not Euclidean distances on a 2D plane but path distances between two points on the dendritic tree. We will precompute these distances and then pass them to the DBSCAN algorithm.

from dendrotweaks.analysis import calculate_pairwise_synaptic_distances

pairwise_distances_matrix = calculate_pairwise_synaptic_distances(
    model.populations['excitatory'].flat_synapses
)

fig, ax = plt.subplots(1, 1, figsize=(6, 5))
cax = ax.matshow(pairwise_distances_matrix, cmap='viridis')
fig.colorbar(cax)
ax.set_title('Pairwise Path Distances Between Synapses')
plt.show()

# Validate properties of the distance matrix
print(f"Matrix symmetric? {np.allclose(pairwise_distances_matrix, pairwise_distances_matrix.T)}")
print(f"Diagonal all zeros? {np.allclose(np.diag(pairwise_distances_matrix), 0)}")

Matrix symmetric? True
Diagonal all zeros? True

fig, ax = plt.subplots(1, 1, figsize=(6, 5))
cutoff = 10
distances_cutoff = np.where(pairwise_distances_matrix > cutoff, cutoff, pairwise_distances_matrix)
cax = ax.matshow(distances_cutoff, cmap='viridis', vmax=cutoff)
fig.colorbar(cax)
ax.set_title(f'Pairwise Path Distances Between Synapses (cutoff at {cutoff} µm)')
plt.show()

Now we are ready to run the DBSCAN algorithm

from sklearn.cluster import DBSCAN

db = DBSCAN(eps=10, min_samples=3, metric='precomputed')
labels = db.fit_predict(pairwise_distances_matrix)

We can visualize the number of synapses in each cluster. Note that the -1 cluster (the outliers) is not shown

fig, ax = plt.subplots(1, 1, figsize=(8, 4))
filtered_labels = labels[labels != -1]  # Exclude noise points
unique_labels, counts = np.unique(filtered_labels, return_counts=True)

colors = [plt.cm.viridis_r(label / max(unique_labels)) for label in unique_labels]
ax.bar(unique_labels, counts, color=colors, edgecolor='k', alpha=0.7)
ax.set_xlabel('Cluster label')
ax.set_ylabel('Number of synapses')
ax.set_title('DBSCAN Clustering of Synapses Based on Path Distance')
plt.show()

We can select the synapses that belong to a single cluster and visualize a smaller pairwise distance matrix for that cluster

def get_synapses_in_cluster(cluster_label, synapses, labels):
    return [syn for syn, label in zip(synapses, labels) if label == cluster_label]

get_synapses_in_cluster(0, synapses, labels)

[<Synapse(NeuronSection(idx=159)(0.026))>,
 <Synapse(NeuronSection(idx=159)(0.297))>,
 <Synapse(NeuronSection(idx=159)(0.587))>,
 <Synapse(NeuronSection(idx=163)(0.009))>,
 <Synapse(NeuronSection(idx=163)(0.077))>,
 <Synapse(NeuronSection(idx=170)(0.008))>,
 <Synapse(NeuronSection(idx=170)(0.036))>,
 <Synapse(NeuronSection(idx=170)(0.150))>,
 <Synapse(NeuronSection(idx=170)(0.303))>,
 <Synapse(NeuronSection(idx=170)(0.304))>,
 <Synapse(NeuronSection(idx=170)(0.363))>]

cluster_0 = get_synapses_in_cluster(0, synapses, labels)
cluster_distances = calculate_pairwise_synaptic_distances(cluster_0)
fig, ax = plt.subplots(1, 1, figsize=(6, 5))
cax = ax.matshow(cluster_distances, cmap='viridis')
fig.colorbar(cax)
ax.set_title('Pairwise Distances Within Cluster 0')
plt.show()

We can finally visualize the synapses on the tuft dendrites and color them by cluster. First, we will visualize them all toghether and then each cluster separately

fig, ax = plt.subplots(1, 1, figsize=(8, 6))

model.sec_tree.plot(ax, show_points=False, show_lines=True, annotate=False)

xs = [coord[0] for coord in syn_coords]
ys = [coord[1] for coord in syn_coords]

for label in set(labels):
    if label == -1:
        ax.scatter([xs[i] for i in range(len(labels)) if labels[i] == label],
                   [ys[i] for i in range(len(labels)) if labels[i] == label],
                   color='w', edgecolor='k', alpha=0.7, label='Noise', s=10, zorder=0)
    else:
        color = plt.cm.viridis_r(label / max(labels))
        
        label_name = f'Cluster {label}'
        ax.scatter([xs[i] for i in range(len(labels)) if labels[i] == label],
                   [ys[i] for i in range(len(labels)) if labels[i] == label],
                   color=color, edgecolor='k', alpha=0.7, label=label_name, s=10)
ax.set_title('Synapse Clusters Based on Path Distance')
ax.legend()
ax.set_xlim(-200, 200)
ax.set_ylim(350, 550)

(350.0, 550.0)

# Plot unclustered synapses
noise_indices = [i for i in range(len(labels)) if labels[i] == -1]
if noise_indices:
    fig, ax = plt.subplots(1, 1, figsize=(8, 6))
    model.sec_tree.plot(ax, show_points=False, show_lines=True, annotate=False)
    ax.scatter([xs[i] for i in noise_indices],
               [ys[i] for i in noise_indices],
               color='w', edgecolor='k', alpha=0.7, label='Unclustered', s=10, zorder=0)
    ax.set_xlim(-200, 200)
    ax.set_ylim(350, 550)
    ax.set_title('Unclustered Synapses')

# Plot a few example clusters
unique_clusters = [label for label in set(labels) if label != -1]
selected_clusters = np.random.choice(unique_clusters, size=min(5, len(unique_clusters)), replace=False)

for label in sorted(selected_clusters):
    fig, ax = plt.subplots(1, 1, figsize=(8, 6))
    model.sec_tree.plot(ax, show_points=False, show_lines=True, annotate=False)
    
    cluster_indices = [i for i in range(len(labels)) if labels[i] == label]
    cluster_distances = pairwise_distances_matrix[np.ix_(cluster_indices, cluster_indices)]
    np.fill_diagonal(cluster_distances, np.inf)  # Ignore self-distances
    min_distances_to_neighbor = np.min(cluster_distances, axis=1)
    max_nearest_neighbor_dist = np.max(min_distances_to_neighbor)
    
    color = plt.cm.viridis_r(label / max(labels))
    
    ax.scatter([xs[i] for i in cluster_indices],
               [ys[i] for i in cluster_indices],
               color=color, edgecolor='k', alpha=0.7, s=10)
    ax.set_xlim(-200, 200)
    ax.set_ylim(350, 550)
    ax.set_title(f'Cluster {label} of {len(cluster_indices)} synapses\n'
                 f'Max nearest-neighbor distance: {max_nearest_neighbor_dist:.2f} µm')