The netsci
package comes with an out-of-box dataset describing the neuronal network among ~2,000 neurons from a rat brain. In this notebook, we shortly explain the nature of this data and demonstrate a basic network analysis applied upon it.
Essentially, after completing this tutorial, you will know:
A deeper description of this network study, together with its neuroscientific context and implication, can be found in this recent paper.
Let’s get started.
You are more than your genes. You are your connectome. --Sebastian Seung
The human brain is composed of about 100 billion specialized nerve cells, neurons. Neurons can connect to up to 10,000 other target neurons each, altogether forming the vast and complex biological neural networks. The inter-connections between neurons are formed by synapses - structural "junctions" that allow transferring electric impulses, action potentials, from one neuron the another. The precise wiring diagram of all neuronal connections, or the connectome, shapes the network-wide electric activity and underlies the different brain functions like information processing, memory, and, ultimately, behavior. However, how exactly does the network structure affect brain activity and function, is a long-standing and still open question in neuroscience for more than a century.
from www.khanacademy.org
The biological connectome can be further modeled as a mathematical directed graph whose nodes correspond to the neurons, and its edges are the synaptic connections between those neurons. This network-based perspective of the brain allows utilizing tools from graph theory and modern network science to uncover key network structures that support brain functions. Such methods have already identified several network features in cortex architecture, such as the rare but highly-connected hub neurons, cliques of all-to-all connected neurons, and overall small-world topology of the cortical microcircuit.
Arguably, one of the most basic, yet still not fully explained, network structures observed in neuronal circuits are the 3-neuron subgraphs (triplets). When counting the frequency of all possible 3-nodes connectivity patterns in the network, the distribution appears highly unexpected. Specifically, few specific configurations (out of all 16 possible) stand-out and are significantly over-expressed (motifs) when compared with reference randomized networks. More surprisingly, different microcircuits across different brain regions commonly display similar over- and under-expression of the same motifs. However, despite their cross-region universality, both the origin of these motifs and their functional implication have remained elusive.
So, what is the origin of neuronal network motifs?
In the absence of a concrete theory about the principles underlying these motifs, their emergence may be hypothetically associated with active plasticity and learning processes. But another, much more parsimonious, possibility comes to mind. Most types of cortical neurons display highly asymmetric geometry with dendritic and axonal trees typically extending in different directions. As a consequence, the probability of forming a connection in one direction (e.g., down) may be higher than in the other direction (e.g., up). This symmetry breaking may distinctively promote some motifs while depressing others. Thus, it could be the case that the geometry per se “enforces” the complex profile of brain microcircuit motifs.
The importance of this hypothesis is that if it is indeed so, then we can now see learning processes as operating on top of an innate, already structured, cortical skeleton rather than on a tabula rasa network connectivity. This will strongly constrain the degree by which plasticity could further shape neural connectivity, and possible reduce the room for learning.
We will test this hypothesis below. Toward this end, we will utilize the publicly available dataset of the Blue Brain model. Specifically, we analyze here a subcircuit of it (pyramidal neurons for layer 5), that is now accessible via netsci
API's.
!pip install plotly
!pip install holoviews
!pip install netsci # required, when notebook is executed in Google Colab
Requirement already satisfied: plotly in /Users/eyalgal/anaconda3/envs/hv/lib/python3.7/site-packages (4.2.1) Requirement already satisfied: six in /Users/eyalgal/anaconda3/envs/hv/lib/python3.7/site-packages (from plotly) (1.12.0) Requirement already satisfied: retrying>=1.3.3 in /Users/eyalgal/anaconda3/envs/hv/lib/python3.7/site-packages (from plotly) (1.3.3) Requirement already satisfied: holoviews in /Users/eyalgal/anaconda3/envs/hv/lib/python3.7/site-packages (1.12.6) Requirement already satisfied: param<2.0,>=1.8.0 in /Users/eyalgal/anaconda3/envs/hv/lib/python3.7/site-packages (from holoviews) (1.9.2) Requirement already satisfied: numpy>=1.0 in /Users/eyalgal/anaconda3/envs/hv/lib/python3.7/site-packages (from holoviews) (1.17.3) Requirement already satisfied: pyviz_comms>=0.7.2 in /Users/eyalgal/anaconda3/envs/hv/lib/python3.7/site-packages (from holoviews) (0.7.2) Requirement already satisfied: netsci in /Users/eyalgal/Dropbox/PhD/Code/python/netsci (0.0.2.2) Requirement already satisfied: numpy>=1.16.2 in /Users/eyalgal/anaconda3/envs/hv/lib/python3.7/site-packages (from netsci) (1.17.3) Requirement already satisfied: pandas>=0.24.2 in /Users/eyalgal/anaconda3/envs/hv/lib/python3.7/site-packages (from netsci) (0.25.2) Requirement already satisfied: matplotlib>=3.0.3 in /Users/eyalgal/anaconda3/envs/hv/lib/python3.7/site-packages (from netsci) (3.1.1) Requirement already satisfied: seaborn>=0.9.0 in /Users/eyalgal/anaconda3/envs/hv/lib/python3.7/site-packages (from netsci) (0.9.0) Requirement already satisfied: networkx>=2.2 in /Users/eyalgal/anaconda3/envs/hv/lib/python3.7/site-packages (from netsci) (2.4) Requirement already satisfied: python-dateutil>=2.6.1 in /Users/eyalgal/anaconda3/envs/hv/lib/python3.7/site-packages (from pandas>=0.24.2->netsci) (2.8.1) Requirement already satisfied: pytz>=2017.2 in /Users/eyalgal/anaconda3/envs/hv/lib/python3.7/site-packages (from pandas>=0.24.2->netsci) (2019.3) Requirement already satisfied: cycler>=0.10 in /Users/eyalgal/anaconda3/envs/hv/lib/python3.7/site-packages (from matplotlib>=3.0.3->netsci) (0.10.0) Requirement already satisfied: kiwisolver>=1.0.1 in /Users/eyalgal/anaconda3/envs/hv/lib/python3.7/site-packages (from matplotlib>=3.0.3->netsci) (1.1.0) Requirement already satisfied: pyparsing!=2.0.4,!=2.1.2,!=2.1.6,>=2.0.1 in /Users/eyalgal/anaconda3/envs/hv/lib/python3.7/site-packages (from matplotlib>=3.0.3->netsci) (2.4.2) Requirement already satisfied: scipy>=0.14.0 in /Users/eyalgal/anaconda3/envs/hv/lib/python3.7/site-packages (from seaborn>=0.9.0->netsci) (1.3.1) Requirement already satisfied: decorator>=4.3.0 in /Users/eyalgal/anaconda3/envs/hv/lib/python3.7/site-packages (from networkx>=2.2->netsci) (4.4.1) Requirement already satisfied: six>=1.5 in /Users/eyalgal/anaconda3/envs/hv/lib/python3.7/site-packages (from python-dateutil>=2.6.1->pandas>=0.24.2->netsci) (1.12.0) Requirement already satisfied: setuptools in /Users/eyalgal/anaconda3/envs/hv/lib/python3.7/site-packages (from kiwisolver>=1.0.1->matplotlib>=3.0.3->netsci) (41.6.0.post20191030)
import numpy as np
import pandas as pd
import holoviews as hv
from holoviews import dim, opts
import matplotlib.pyplot as plt
import seaborn as sns; # sns.set_context("paper")
from netsci.datasets import load_connectome
import netsci.metrics.motifs as nsm
hv.extension('plotly')
plt.rcParams['figure.dpi'] = 100
The connectome dataset is composed of two parts:
gid
-s and 3d position $ (x,y,z) $.connectome = load_connectome()
display({k: type(v) for k, v in connectome.items()})
{'title': str, 'nodes': pandas.core.frame.DataFrame, 'edges': pandas.core.frame.DataFrame}
The magnitude of the network is
print(connectome['title'])
neurons, synapses = connectome['nodes'], connectome['edges']
n = len(neurons)
print(f'\t{n:,} neurons\n',
f'\t{len(synapses):,} connections (p = {len(synapses)/(n*(n-1)):.2%})\n'
f'\t{synapses["contacts"].sum():,} contacts')
Connectome (L5-TTPC) 2,003 neurons 102,732 connections (p = 2.56%) 650,949 contacts
display(neurons)
gid | x | y | z | |
---|---|---|---|---|
0 | 75211 | 341.846181 | 773.977844 | 636.437126 |
1 | 75214 | 305.869316 | 1200.585727 | 608.787329 |
2 | 75217 | 312.411025 | 1184.177732 | 647.120290 |
3 | 75218 | 275.261806 | 1110.341752 | 630.010095 |
4 | 75221 | 239.869853 | 880.629815 | 609.945456 |
... | ... | ... | ... | ... |
1998 | 81314 | 472.198046 | 1087.211926 | 409.227347 |
1999 | 81317 | 308.608095 | 947.743964 | 624.308077 |
2000 | 81318 | 325.510390 | 1029.783942 | 574.410863 |
2001 | 81319 | 383.361952 | 1161.047906 | 655.450589 |
2002 | 81322 | 331.464765 | 1062.599933 | 651.448519 |
2003 rows × 4 columns
The y-axis denotes the cortical depth. The positive direction is defined to be the pia direction in the column.
# hv.Scatter3D((neurons.x, neurons.y, neurons.z)).opts(
# opts.Scatter3D(azimuth=40, elevation=20, color='z', s=50, cmap='fire'))
hv.Scatter3D((neurons.x, neurons.y, neurons.z)) \
.opts(cmap='fire', color='y', size=3, alpha=.75, width=800, height=800)
display(synapses)
from | to | contacts | |
---|---|---|---|
0 | 75211 | 75221 | 5 |
1 | 75211 | 75276 | 5 |
2 | 75211 | 75473 | 5 |
3 | 75211 | 76012 | 7 |
4 | 75211 | 76095 | 3 |
... | ... | ... | ... |
102727 | 81322 | 80991 | 5 |
102728 | 81322 | 81051 | 5 |
102729 | 81322 | 81111 | 4 |
102730 | 81322 | 81116 | 9 |
102731 | 81322 | 81147 | 4 |
102732 rows × 3 columns
Although representing the network as an edge list as above might space-efficient, sometimes it is more convenient to analyze the adjacency matrix directly. You could, obviously, transform the edge list to the matrix representation yourself (a quick exercise in table operations acrobatics, or straightforward by nested loops). Still, we can simply get it using the adjacency=True
flag as below:
connectome = load_connectome(adjacency=True)
neurons, synapses, A, W = connectome['nodes'], connectome['edges'], connectome['A'], connectome['W']
print(f'{A.shape[0]}x{A.shape[1]} adjacency matrix\n',
f'\t{len(A):,} neurons\n',
f'\t{A.sum():,} connections\n'
f'\t{W.sum():,} contacts')
2003x2003 adjacency matrix 2,003 neurons 102,732 connections 650,949 contacts
Note that for convenience, in addition to the binary adjacency matrix A
, a weighted matrix W
, whose values describe the number of synapses per connection, is also returned. For the rest of our analysis here, we consider only the connections (without the contacts).
Let's examine the how each motif is embedded in space.
tid_order = nsm.triad_order_nn4576
frequency, participations = nsm.motifs(A, participation=True)
frequency, participations = np.take(frequency, tid_order), np.take(participations, tid_order)
print(frequency)
[4187617 2738360 3162677 281224 208612 174242 29961 4328 7555 7808 6694 689 15]
gids = neurons['gid'].values
motifs = pd.concat([pd.DataFrame(gids.take(p), columns=['R','G','B']).assign(motif=i+1)
for i, p in enumerate(participations) if p], ignore_index=True)
display(motifs)
R | G | B | motif | |
---|---|---|---|---|
0 | 75211 | 75221 | 75261 | 1 |
1 | 75211 | 75221 | 75271 | 1 |
2 | 75211 | 75221 | 75277 | 1 |
3 | 75211 | 75221 | 75396 | 1 |
4 | 75211 | 75221 | 75691 | 1 |
... | ... | ... | ... | ... |
10809777 | 81129 | 80370 | 76956 | 13 |
10809778 | 78718 | 77948 | 77520 | 13 |
10809779 | 80751 | 79899 | 77546 | 13 |
10809780 | 80386 | 80335 | 77704 | 13 |
10809781 | 78422 | 77773 | 77708 | 13 |
10809782 rows × 4 columns
assert np.array_equal(frequency, motifs.groupby('motif').size())
motifs.groupby('motif').size().to_frame('Frequency').T.style.format('{:,}')
motif | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | 11 | 12 | 13 |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|
Frequency | 4,187,617 | 2,738,360 | 3,162,677 | 281,224 | 208,612 | 174,242 | 29,961 | 4,328 | 7,555 | 7,808 | 6,694 | 689 | 15 |
Lastly, we examine the positions of the neurons forming each of the motifs. Since the pyramidal neurons are vertical, we'll focus on the y-coordinate (cortical depth).
embedding = motifs \
.merge(neurons.set_index('gid').add_prefix('R'), how='left', left_on='R', right_index=True) \
.merge(neurons.set_index('gid').add_prefix('G'), how='left', left_on='G', right_index=True) \
.merge(neurons.set_index('gid').add_prefix('B'), how='left', left_on='B', right_index=True)
embedding
R | G | B | motif | Rx | Ry | Rz | Gx | Gy | Gz | Bx | By | Bz | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|
0 | 75211 | 75221 | 75261 | 1 | 341.846181 | 773.977844 | 636.437126 | 239.869853 | 880.629815 | 609.945456 | 203.276337 | 968.822791 | 481.909118 |
1 | 75211 | 75221 | 75271 | 1 | 341.846181 | 773.977844 | 636.437126 | 239.869853 | 880.629815 | 609.945456 | 320.183623 | 1180.844858 | 570.809167 |
2 | 75211 | 75221 | 75277 | 1 | 341.846181 | 773.977844 | 636.437126 | 239.869853 | 880.629815 | 609.945456 | 293.199873 | 1090.600882 | 580.749479 |
3 | 75211 | 75221 | 75396 | 1 | 341.846181 | 773.977844 | 636.437126 | 239.869853 | 880.629815 | 609.945456 | 393.450029 | 1067.014389 | 642.183122 |
4 | 75211 | 75221 | 75691 | 1 | 341.846181 | 773.977844 | 636.437126 | 239.869853 | 880.629815 | 609.945456 | 497.137587 | 1024.840713 | 725.135217 |
... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
10809777 | 81129 | 80370 | 76956 | 13 | 192.844296 | 948.897651 | 737.201974 | 178.475477 | 1061.766714 | 692.749867 | 248.743735 | 1104.252847 | 703.846591 |
10809778 | 78718 | 77948 | 77520 | 13 | 472.534544 | 787.782031 | 597.521670 | 385.660359 | 797.075622 | 501.540998 | 359.745437 | 878.731038 | 649.375373 |
10809779 | 80751 | 79899 | 77546 | 13 | 214.454828 | 907.813569 | 437.465360 | 226.882608 | 883.506021 | 430.744574 | 251.492190 | 1077.677984 | 499.352161 |
10809780 | 80386 | 80335 | 77704 | 13 | 335.823705 | 1079.200210 | 626.482609 | 340.106908 | 1018.695726 | 604.070499 | 337.106262 | 998.329943 | 660.917740 |
10809781 | 78422 | 77773 | 77708 | 13 | 277.609104 | 1140.425748 | 652.055906 | 271.051619 | 1144.976403 | 535.146850 | 355.084509 | 1117.287910 | 510.178296 |
10809782 rows × 13 columns
motifs_order = [3, 4, 7, 10, 11]
axes = embedding[embedding['motif'].isin(motifs_order)].groupby('motif') \
.boxplot(column=['Ry', 'Gy', 'By'], layout=(1, len(motifs_order)), figsize=(12,3),
showfliers=False, patch_artist=True, boxprops=dict())
colors = ['pink', 'lightgreen', 'lightblue']
for ax in axes.values:
plt.setp(ax.lines, color='k')
[patch.set(edgecolor='k', facecolor=color) for patch, color in zip(ax.artists, colors)]
axes.values[0].set_ylabel('Height (µm)')
Text(0, 0.5, 'Height (µm)')
The analysis uncovers a close match between the physical positions of the neurons and the synaptic connectivity they form. For instance, for motif #4 the “source” neuron (the presynaptic neuron that projects to the other two) is located, on average, above the other two cells whereas the “sink” neuron (the postsynaptic neuron receiving projections from the other two) is, on average, the lowest.
Books
Papers
Web resources