Part III: Computational Neuroscience | Chapter 4

Connectomics

Graph theory, small-world networks, rich-club organization, and the C. elegans connectome

Mapping the Brain's Wiring

Connectomics aims to comprehensively map neural connections at scales ranging from individual synapses (electron microscopy connectomes) to brain regions (macroscale connectomes from diffusion MRI). Graph theory provides the mathematical framework for analyzing these networks, revealing organizational principles such as small-world topology, modular structure, and rich-club connectivity.

The C. elegans connectome — the only complete connectome of a nervous system — has served as a benchmark for network analysis methods. This chapter covers graph-theoretic measures, generative models of brain networks, and the relationship between structural connectivity and functional dynamics.

1. Graph Theory for Brain Networks

A brain network is represented as a graph $G = (V, E)$ where nodes $V$ represent neurons or brain regions and edges $E$ represent connections. The adjacency matrix$A_{ij} = 1$ if nodes $i$ and $j$ are connected (or a weight for weighted graphs).

Derivation 1: Clustering Coefficient and Path Length

The clustering coefficient measures local connectivity. For node $i$ with degree$k_i$ and $t_i$ triangles:

$$C_i = \frac{2 t_i}{k_i(k_i - 1)}, \quad \bar{C} = \frac{1}{N}\sum_i C_i$$

The characteristic path length is the average shortest path between all node pairs:

$$L = \frac{1}{N(N-1)} \sum_{i \neq j} d(i, j)$$

For a random graph (Erdos-Renyi) with connection probability $p$:

$$C_{\text{random}} = p, \quad L_{\text{random}} \approx \frac{\ln N}{\ln(pN)}$$

A small-world network (Watts and Strogatz, 1998) has high clustering ($C \gg C_{\text{random}}$) but short path length ($L \approx L_{\text{random}}$). The small-world index$\sigma = (C/C_{\text{random}}) / (L/L_{\text{random}}) > 1$ quantifies this property. Brain networks consistently show small-world topology, combining local processing efficiency with global integration.

2. Rich-Club Organization

Rich-club organization describes the tendency of high-degree nodes (hubs) to be preferentially interconnected. In brain networks, rich-club hubs include regions like the precuneus, superior frontal cortex, and superior parietal cortex — areas critical for integrating information across functional domains.

Derivation 2: Rich-Club Coefficient

The rich-club coefficient at degree threshold $k$ is the fraction of possible edges among nodes with degree greater than $k$:

$$\phi(k) = \frac{2 E_{>k}}{N_{>k}(N_{>k} - 1)}$$

where $N_{>k}$ is the number of nodes with degree greater than $k$ and $E_{>k}$ is the number of edges among them. To assess significance, we compare against degree-preserving random networks. The normalized rich-club coefficient is:

$$\phi_{\text{norm}}(k) = \frac{\phi(k)}{\langle\phi_{\text{random}}(k)\rangle}$$

Values $\phi_{\text{norm}}(k) > 1$ indicate rich-club organization exceeding what degree distribution alone would predict. Van den Heuvel and Sporns (2011) showed that the human connectome has significant rich-club organization, with hub regions forming a densely connected "backbone" supporting global brain communication.

Derivation 3: Modularity and Community Detection

Brain networks are organized into modules (communities) corresponding to functional systems. The modularity quality function (Newman, 2004) is:

$$Q = \frac{1}{2m}\sum_{ij}\left[A_{ij} - \frac{k_i k_j}{2m}\right]\delta(c_i, c_j)$$

where $m$ is the total number of edges, $k_i$ is the degree of node $i$,$c_i$ is the community assignment, and $\delta$ is the Kronecker delta. The null model $k_i k_j / 2m$ is the expected edge weight in a random graph preserving degrees. Maximizing $Q$ (NP-hard, but approximated by the Louvain algorithm) reveals brain modules corresponding to:

Visual system (occipital regions)
Sensorimotor system (pre/postcentral)
Default mode network (medial frontal, posterior cingulate)
Frontoparietal control network

Typical brain network modularity is $Q \approx 0.4\text{--}0.6$, significantly higher than random networks ($Q \approx 0$).

3. The C. elegans Connectome

The nematode C. elegans has 302 neurons with approximately 7,000 chemical synapses and 900 gap junctions, fully mapped by White et al. (1986). This connectome serves as the gold standard for testing network analysis methods and understanding structure-function relationships in neural circuits.

Derivation 4: Network Motifs and Their Significance

Network motifs are recurrent subgraph patterns that appear significantly more often than in random networks. For a directed graph with $N$ nodes, the number of possible 3-node motifs is 13. The z-score for motif $m$ is:

$$Z_m = \frac{n_m - \langle n_m^{\text{rand}} \rangle}{\sigma_m^{\text{rand}}}$$

where $n_m$ is the observed count and $\langle n_m^{\text{rand}} \rangle$ is the mean across randomized networks. The significance profile (normalized z-scores) is:

$$SP_m = \frac{Z_m}{\sqrt{\sum_{m'} Z_{m'}^2}}$$

In C. elegans, feedforward loops, bi-parallel motifs, and reciprocal connections are over-represented. These motifs serve computational functions: feedforward loops enable signal persistence, reciprocal connections support oscillations, and convergent motifs implement coincidence detection.

Derivation 5: Spectral Analysis of the Graph Laplacian

The graph Laplacian $\mathbf{L} = \mathbf{D} - \mathbf{A}$ (where $\mathbf{D}$ is the degree matrix) captures the diffusive dynamics on a network. The eigendecomposition:

$$\mathbf{L} = \sum_{k=0}^{N-1} \lambda_k \mathbf{u}_k \mathbf{u}_k^T, \quad 0 = \lambda_0 \leq \lambda_1 \leq \ldots \leq \lambda_{N-1}$$

The second-smallest eigenvalue $\lambda_1$ (algebraic connectivity or Fiedler value) measures how well-connected the graph is. The Fiedler vector $\mathbf{u}_1$ provides an optimal bipartition. Diffusion on the graph follows:

$$\mathbf{x}(t) = \sum_k e^{-\lambda_k t} (\mathbf{u}_k^T \mathbf{x}_0) \mathbf{u}_k$$

Slow modes ($\lambda_k$ small) correspond to large-scale functional networks, while fast modes correspond to local activity. This spectral approach predicts functional connectivity (correlated BOLD signals) from structural connectivity, with remarkable accuracy (Honey et al., 2009).

4. Historical Development

• 1986: White et al. publish the complete connectome of C. elegans, the first whole-nervous-system wiring diagram.
• 1998: Watts and Strogatz introduce the small-world network model, demonstrating its prevalence in real networks.
• 2002: Sporns, Tononi, and Kotter coin the term "connectome" for a comprehensive map of neural connections.
• 2004: Newman introduces the modularity function, enabling systematic community detection in complex networks.
• 2009: The Human Connectome Project is launched, mapping structural and functional connectivity in 1,200 subjects.
• 2011: Van den Heuvel and Sporns discover rich-club organization in the human structural connectome.
• 2019: Zheng et al. publish a complete electron microscopy connectome of the Drosophila brain hemisphere (~25,000 neurons).
• 2024: Large-scale EM connectomes of mouse cortex provide circuit-level wiring diagrams of mammalian cortex.

5. Applications

Neurological Biomarkers

Graph-theoretic measures from connectomics serve as biomarkers for Alzheimer's disease (reduced hub connectivity), schizophrenia (altered modularity), and traumatic brain injury (disrupted rich-club). Network-based statistics identify disease-affected subnetworks.

Brain Stimulation Planning

Connectome-informed models predict how stimulation of one region propagates through the network. This guides optimal electrode placement for DBS and TMS targeting in depression and OCD treatment.

AI Architecture Design

Brain-inspired network topologies (small-world, modular) improve artificial neural network efficiency. Connectome-constrained models achieve better generalization with fewer parameters.

Surgical Planning

Connectome analysis predicts cognitive outcomes of brain surgery by identifying which network connections would be disrupted. This enables personalized risk assessment for epilepsy surgery and tumor resection.

6. Computational Exploration

Connectomics: Small-World Networks, Rich-Club, and Graph Spectral Analysis

Python

script.py321 lines

import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt

print("=" * 72)
print("CONNECTOMICS: GRAPH THEORY, SMALL-WORLD, AND RICH-CLUB")
print("=" * 72)

np.random.seed(42)

# --------------------------------------------------
# 1. SMALL-WORLD NETWORK (WATTS-STROGATZ)
# --------------------------------------------------
print()
print("1. WATTS-STROGATZ SMALL-WORLD MODEL")
print("-" * 55)

def watts_strogatz(N, K, p):
    """Generate Watts-Strogatz small-world network."""
    A = np.zeros((N, N))
    # Regular ring lattice
    for i in range(N):
        for j in range(1, K//2 + 1):
            A[i, (i+j) % N] = 1
            A[i, (i-j) % N] = 1
    # Rewire with probability p
    for i in range(N):
        for j in range(1, K//2 + 1):
            if np.random.random() < p:
                A[i, (i+j) % N] = 0
                A[(i+j) % N, i] = 0
                new_target = np.random.randint(N)
                while new_target == i or A[i, new_target] == 1:
                    new_target = np.random.randint(N)
                A[i, new_target] = 1
                A[new_target, i] = 1
    return A

def clustering_coeff(A):
    N = A.shape[0]
    cc = 0
    for i in range(N):
        neighbors = np.where(A[i] > 0)[0]
        ki = len(neighbors)
        if ki < 2:
            continue
        triangles = 0
        for ni in range(len(neighbors)):
            for nj in range(ni+1, len(neighbors)):
                if A[neighbors[ni], neighbors[nj]] > 0:
                    triangles += 1
        cc += 2 * triangles / (ki * (ki - 1))
    return cc / N

def avg_path_length(A):
    N = A.shape[0]
    # BFS for shortest paths
    total = 0
    count = 0
    for source in range(N):
        dist = -np.ones(N)
        dist[source] = 0
        queue = [source]
        while queue:
            node = queue.pop(0)
            for neighbor in np.where(A[int(node)] > 0)[0]:
                if dist[neighbor] == -1:
                    dist[neighbor] = dist[int(node)] + 1
                    queue.append(neighbor)
        for d in dist:
            if d > 0:
                total += d
                count += 1
    return total / max(count, 1)

N_ws = 100
K_ws = 8
rewire_probs = [0, 0.001, 0.01, 0.05, 0.1, 0.3, 0.5, 1.0]
cc_vals = []
apl_vals = []

for p in rewire_probs:
    A = watts_strogatz(N_ws, K_ws, p)
    cc = clustering_coeff(A)
    apl = avg_path_length(A)
    cc_vals.append(cc)
    apl_vals.append(apl)
    print("  p = {:.3f}: C = {:.3f}, L = {:.2f}".format(p, cc, apl))

# Normalize
cc_norm = np.array(cc_vals) / cc_vals[0]
apl_norm = np.array(apl_vals) / apl_vals[0]

# --------------------------------------------------
# 2. RICH-CLUB COEFFICIENT
# --------------------------------------------------
print()
print("2. RICH-CLUB ANALYSIS")
print("-" * 55)

# Create a network with rich-club structure
N_rc = 80
A_rc = np.zeros((N_rc, N_rc))

# Hub nodes (first 10) are densely interconnected
n_hubs = 10
for i in range(n_hubs):
    for j in range(i+1, n_hubs):
        if np.random.random() < 0.8:
            A_rc[i,j] = A_rc[j,i] = 1

# Hub-to-periphery connections
for i in range(n_hubs):
    for j in range(n_hubs, N_rc):
        if np.random.random() < 0.15:
            A_rc[i,j] = A_rc[j,i] = 1

# Periphery-to-periphery (sparse)
for i in range(n_hubs, N_rc):
    for j in range(i+1, N_rc):
        if np.random.random() < 0.03:
            A_rc[i,j] = A_rc[j,i] = 1

degrees_rc = np.sum(A_rc, axis=1)
k_thresholds = range(1, int(np.max(degrees_rc)))
phi_vals = []
phi_rand_vals = []

for k_thresh in k_thresholds:
    high_deg = np.where(degrees_rc > k_thresh)[0]
    n_high = len(high_deg)
    if n_high < 2:
        break
    edges_among = 0
    for i in high_deg:
        for j in high_deg:
            if i < j:
                edges_among += A_rc[i, j]
    possible = n_high * (n_high - 1) / 2
    phi = edges_among / max(possible, 1)
    phi_vals.append(phi)

# Random null (approximate)
    total_edges = np.sum(A_rc) / 2
    p_rand = 2 * total_edges / (N_rc * (N_rc - 1))
    phi_rand_vals.append(p_rand * n_high)

print("  Network: {} nodes, {} hubs".format(N_rc, n_hubs))
print("  Hub mean degree: {:.1f}".format(np.mean(degrees_rc[:n_hubs])))
print("  Periphery mean degree: {:.1f}".format(np.mean(degrees_rc[n_hubs:])))

# --------------------------------------------------
# 3. MODULARITY
# --------------------------------------------------
print()
print("3. MODULARITY ANALYSIS")
print("-" * 55)

# Create modular network
N_mod = 60
n_modules = 3
module_size = N_mod // n_modules
A_mod = np.zeros((N_mod, N_mod))
true_labels = np.zeros(N_mod, dtype=int)

for mod in range(n_modules):
    start = mod * module_size
    end = start + module_size
    true_labels[start:end] = mod
    for i in range(start, end):
        for j in range(i+1, end):
            if np.random.random() < 0.3:  # within-module
                A_mod[i,j] = A_mod[j,i] = 1
        for j in range(end, N_mod):
            if np.random.random() < 0.02:  # between-module
                A_mod[i,j] = A_mod[j,i] = 1

# Compute modularity
m_total = np.sum(A_mod) / 2
degrees_mod = np.sum(A_mod, axis=1)
Q = 0
for i in range(N_mod):
    for j in range(N_mod):
        if true_labels[i] == true_labels[j]:
            Q += A_mod[i,j] - degrees_mod[i] * degrees_mod[j] / (2 * m_total)
Q /= (2 * m_total)

print("  Modules: {}, Module size: {}".format(n_modules, module_size))
print("  Modularity Q = {:.4f}".format(Q))
print("  Within-module density: {:.3f}".format(0.3))
print("  Between-module density: {:.3f}".format(0.02))

# --------------------------------------------------
# 4. GRAPH LAPLACIAN SPECTRAL ANALYSIS
# --------------------------------------------------
print()
print("4. GRAPH LAPLACIAN SPECTRAL ANALYSIS")
print("-" * 55)

D_mod = np.diag(degrees_mod)
L_mod = D_mod - A_mod
eigenvalues, eigenvectors = np.linalg.eigh(L_mod)

print("  First 6 eigenvalues:")
for k in range(6):
    print("    lambda_{} = {:.4f}".format(k, eigenvalues[k]))
print("  Fiedler value (lambda_1): {:.4f}".format(eigenvalues[1]))
print("  Spectral gap: {:.4f}".format(eigenvalues[1]))

# --------------------------------------------------
# 5. C. ELEGANS-LIKE NETWORK ANALYSIS
# --------------------------------------------------
print()
print("5. C. ELEGANS-LIKE NETWORK SIMULATION")
print("-" * 55)

N_ce = 100  # simplified version
# Create network with C. elegans-like properties
A_ce = np.zeros((N_ce, N_ce))
# Sensory neurons (0-29), interneurons (30-69), motor neurons (70-99)
for i in range(30):
    targets = np.random.choice(range(30, 70), size=3, replace=False)
    for t in targets:
        A_ce[i, t] = 1
for i in range(30, 70):
    targets_inter = np.random.choice(range(30, 70), size=2, replace=False)
    targets_motor = np.random.choice(range(70, 100), size=2, replace=False)
    for t in np.concatenate([targets_inter, targets_motor]):
        if t != i:
            A_ce[i, int(t)] = 1
for i in range(70, 100):
    targets = np.random.choice(range(70, 100), size=1)
    for t in targets:
        if t != i:
            A_ce[i, int(t)] = 1

degrees_ce = np.sum(A_ce, axis=1) + np.sum(A_ce, axis=0)
print("  Neurons: {} (30 sensory, 40 inter, 30 motor)".format(N_ce))
print("  Edges: {}".format(int(np.sum(A_ce))))
print("  Mean degree: {:.1f}".format(np.mean(degrees_ce)))
print("  Sensory mean out-degree: {:.1f}".format(np.mean(np.sum(A_ce[:30], axis=1))))
print("  Interneuron mean out-degree: {:.1f}".format(np.mean(np.sum(A_ce[30:70], axis=1))))
print("  Motor mean out-degree: {:.1f}".format(np.mean(np.sum(A_ce[70:], axis=1))))

# --------------------------------------------------
# PLOTTING
# --------------------------------------------------
fig, axes = plt.subplots(2, 3, figsize=(16, 10))
fig.suptitle("Connectomics: Small-World, Rich-Club, Modularity, and Spectral Analysis",
             fontsize=14, fontweight='bold', color='white')
fig.patch.set_facecolor('#0a0a0a')

for ax in axes.flat:
    ax.set_facecolor('#111111')
    ax.tick_params(colors='white', labelsize=8)
    ax.xaxis.label.set_color('white')
    ax.yaxis.label.set_color('white')
    ax.title.set_color('#f0abfc')
    for spine in ax.spines.values():
        spine.set_color('#333333')

# Panel 1: Small-world properties
axes[0,0].semilogx(rewire_probs[1:], cc_norm[1:], 'o-', color='#ec4899', lw=2, markersize=6, label='C/C(0)')
axes[0,0].semilogx(rewire_probs[1:], apl_norm[1:], 's-', color='#60a5fa', lw=2, markersize=6, label='L/L(0)')
axes[0,0].axhspan(0.1, 0.5, alpha=0.1, color='#fbbf24')
axes[0,0].set_xlabel('Rewiring probability p')
axes[0,0].set_ylabel('Normalized value')
axes[0,0].set_title('Watts-Strogatz Small World')
axes[0,0].legend(fontsize=7, facecolor='#1a1a1a', edgecolor='#333', labelcolor='white')

# Panel 2: Degree distribution
axes[0,1].hist(degrees_rc[:n_hubs], bins=10, color='#ec4899', alpha=0.7, label='Hubs', density=True)
axes[0,1].hist(degrees_rc[n_hubs:], bins=10, color='#60a5fa', alpha=0.7, label='Periphery', density=True)
axes[0,1].set_xlabel('Degree')
axes[0,1].set_ylabel('Density')
axes[0,1].set_title('Rich-Club: Degree Distribution')
axes[0,1].legend(fontsize=7, facecolor='#1a1a1a', edgecolor='#333', labelcolor='white')

# Panel 3: Rich-club coefficient
k_range = range(1, len(phi_vals) + 1)
axes[0,2].plot(list(k_range), phi_vals, 'o-', color='#d946ef', lw=2, markersize=5)
axes[0,2].set_xlabel('Degree threshold k')
axes[0,2].set_ylabel('Rich-club coeff phi(k)')
axes[0,2].set_title('Rich-Club Coefficient')

# Panel 4: Adjacency matrix (modular)
reorder = np.argsort(true_labels)
A_reordered = A_mod[np.ix_(reorder, reorder)]
axes[1,0].imshow(A_reordered, cmap='magma', aspect='equal', interpolation='none')
for i in range(1, n_modules):
    axes[1,0].axhline(i * module_size - 0.5, color='white', alpha=0.5, lw=0.5)
    axes[1,0].axvline(i * module_size - 0.5, color='white', alpha=0.5, lw=0.5)
axes[1,0].set_title('Modular Network (Q={:.3f})'.format(Q))
axes[1,0].set_xlabel('Node')
axes[1,0].set_ylabel('Node')

# Panel 5: Laplacian eigenvalue spectrum
axes[1,1].plot(eigenvalues[:20], 'o-', color='#ec4899', lw=2, markersize=5)
axes[1,1].axhline(0, color='white', alpha=0.3)
axes[1,1].set_xlabel('Eigenvalue index')
axes[1,1].set_ylabel('Eigenvalue')
axes[1,1].set_title('Laplacian Spectrum (first 20)')

# Panel 6: Fiedler vector (community structure)
fiedler = eigenvectors[:, 1]
colors_fiedler = ['#ec4899' if true_labels[i] == 0 else '#d946ef' if true_labels[i] == 1 else '#60a5fa'
                  for i in range(N_mod)]
axes[1,2].scatter(range(N_mod), fiedler, c=colors_fiedler, s=20, alpha=0.8)
axes[1,2].axhline(0, color='white', alpha=0.3, linestyle='--')
axes[1,2].set_xlabel('Node index')
axes[1,2].set_ylabel('Fiedler vector value')
axes[1,2].set_title('Fiedler Vector (Graph Bipartition)')

plt.tight_layout()
plt.savefig('output.png', dpi=150, bbox_inches='tight', facecolor='#0a0a0a')
plt.close()
print()
print("[Plot saved: small-world, degree distribution, rich-club,")
print(" modular adjacency, Laplacian spectrum, and Fiedler vector]")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

Chapter Summary

• Small-world topology: brain networks have high clustering $C \gg C_{\text{random}}$ with short path length $L \approx L_{\text{random}}$.
• Rich-club: hub regions are densely interconnected, forming a communication backbone with $\phi_{\text{norm}} > 1$.
• Modularity: brain networks decompose into functional modules (visual, motor, default mode) with $Q \approx 0.4\text{--}0.6$.
• Spectral analysis: the graph Laplacian eigenvalues predict diffusion dynamics; the Fiedler vector reveals community structure.
• C. elegans: the complete connectome (302 neurons) demonstrates feedforward, recurrent, and modular organization principles.

← Brain Imaging Neurological Disorders →