Embedding a neuron model in a network model¶

This tutorial guides the user through how to use ISF to embed a neuron model in a dense connectome model, mapping synapses onto a neuron morphology based on empirical data. We provide a network modeling pipeline that has been described in detail in Udvary et al. (2022).

To generate such synapse maps from scratch, the user must download the barrel cortex model from harvard Dataverse:

pixi run download_bc_model

To get started, adapt the desired output directory below:

[1]:

from pathlib import Path
tutorial_output_dir = f"{Path.home()}/isf_tutorial_output"  # <-- Change this to your desired output directory

[2]:

%matplotlib inline
import Interface as I
from getting_started import getting_started_dir
example_data_dir = I.os.path.join(getting_started_dir, 'example_data')
db = I.DataBase(tutorial_output_dir).create_sub_db("network_modeling")

[INFO] ISF: Current version: heads/data+0.gebc00ec75.dirty
[INFO] ISF: Current pid: 100207
[INFO] ISF: Loading mechanisms:
[ATTENTION] ISF: The source folder has uncommited changes!



[INFO] ISF: Loaded modules with __version__ attribute are:
IPython: 8.12.2, Interface: heads/data+0.gebc00ec75.dirty, PIL: 10.4.0, _brotli: 1.0.9, _csv: 1.0, _ctypes: 1.1.0, _curses: b'2.2', _decimal: 1.70, argparse: 1.1, backcall: 0.2.0, blosc: 1.11.1, bluepyopt: 1.9.126, brotli: 1.0.9, certifi: 2024.08.30, cffi: 1.17.0, charset_normalizer: 3.4.0, click: 7.1.2, cloudpickle: 3.1.0, colorama: 0.4.6, comm: 0.2.2, csv: 1.0, ctypes: 1.1.0, cycler: 0.12.1, cytoolz: 0.12.3, dash: 2.18.2, dask: 2.30.0, dateutil: 2.9.0, deap: 1.4, debugpy: 1.8.5, decimal: 1.70, decorator: 5.1.1, defusedxml: 0.7.1, distributed: 2.30.0, distutils: 3.8.20, django: 1.8.19, entrypoints: 0.4, executing: 2.1.0, fasteners: 0.17.3, flask: 1.1.4, fsspec: 2024.10.0, future: 1.0.0, greenlet: 3.1.1, idna: 3.10, ipaddress: 1.0, ipykernel: 6.29.5, ipywidgets: 8.1.5, isf_pandas_msgpack: 0.2.2, itsdangerous: 1.1.0, jedi: 0.19.1, jinja2: 2.11.3, joblib: 1.4.2, json: 2.0.9, jupyter_client: 7.3.4, jupyter_core: 5.7.2, kiwisolver: 1.4.5, logging: 0.5.1.2, markupsafe: 2.0.1, matplotlib: 3.5.1, matplotlib_inline: 0.1.7, msgpack: 1.0.8, neuron: 7.8.2+, numcodecs: 0.12.1, numexpr: 2.8.6, numpy: 1.19.2, packaging: 24.2, pandas: 1.1.3, parameters: 0.2.1, parso: 0.8.4, past: 1.0.0, pexpect: 4.9.0, pickleshare: 0.7.5, platform: 1.0.8, platformdirs: 4.3.6, plotly: 5.24.1, prompt_toolkit: 3.0.48, psutil: 6.0.0, ptyprocess: 0.7.0, pure_eval: 0.2.3, pydevd: 2.9.5, pygments: 2.18.0, pyparsing: 3.1.4, pytz: 2024.2, re: 2.2.1, requests: 2.32.3, scandir: 1.10.0, scipy: 1.5.2, seaborn: 0.12.2, six: 1.16.0, sklearn: 0.23.2, socketserver: 0.4, socks: 1.7.1, sortedcontainers: 2.4.0, stack_data: 0.6.2, statsmodels: 0.13.2, sumatra: 0.7.4, tables: 3.8.0, tblib: 3.0.0, tlz: 0.12.3, toolz: 1.0.0, tqdm: 4.67.1, traitlets: 5.14.3, urllib3: 2.2.3, wcwidth: 0.2.13, werkzeug: 1.0.1, yaml: 5.3.1, zarr: 2.15.0, zlib: 1.0, zmq: 26.2.0, zstandard: 0.19.0

Registering the cell morphology in the desired reference frame¶

The first step to embedding a neuron model into a network model, is to align the reference frames of both.

From external resources (e.g. the NeuroMorph pipeline), you need a morphology .hoc-file. The coordinates in the hoc morphology file need to be anchored at the desired location.

As an example, we provide such file of a Layer-5 Pyramidal Tract Neuron (L5PT), whose coordinates are anchored in the C2 column of a rat barrel cortex (BC)

[3]:

anatomy_dir = I.os.path.join(example_data_dir, 'anatomical_constraints')
I.os.listdir(anatomy_dir)

[3]:

['89_L5_CDK20050712_nr6L5B_dend_PC_neuron_transform_registered_C2.hoc',
 '86_L5_CDK20041214_nr3L5B_dend_PC_neuron_transform_registered_C2center_scaled_diameters.hoc',
 '86_L5_CDK20041214_nr3L5B_dend_PC_neuron_transform_registered_C2center.swc',
 '86_L5_CDK20041214_nr3L5B_dend_PC_neuron_transform_registered_C2center.hoc',
 'example_embedding_86_C2_center',
 'presynaptic_somata_soma_locations.csv',
 'NumberOfConnectedCells.csv']

[4]:

path_to_hoc = I.os.path.join(
    anatomy_dir,
    '86_L5_CDK20041214_nr3L5B_dend_PC_neuron_transform_registered_C2center.hoc')
path_to_scaled_hoc = I.os.path.join(
    anatomy_dir,
    'anatomical_constraints', \
    '86_L5_CDK20041214_nr3L5B_dend_PC_neuron_transform_registered_C2center_scaled_diameters.hoc')

We copy these morphology files to our DataBase.

[5]:

if not 'anatomical_constraints' in db.keys():
    db.create_managed_folder('anatomical_constraints')
    I.shutil.copy(path_to_hoc, db['anatomical_constraints'])
    I.shutil.copy(path_to_scaled_hoc, db['anatomical_constraints'])

Calculating synapse positions on the neuron morphology¶

We use singlecell_input_mapper to create an anatomical model of how that cell is integrated in the brain are (here: the barrel cortex). For more information on how this is done, see Egger et al. 2014. This module creates an anatomical reconstruction of axo-dendritic connections depending on the spatial distribution of cells, bouton density, and anatomical constraints of post-synaptic targets of the postsynaptic cell, such as morphology and synapse density. To do so, it needs the following input:

1. The neuron morphology (the .hoc file that we already copied over)

2. Density of total Post-Synaptic Target sites (PST) across the dendritic tree of the postynaptic cell (used for normalization)

3. Anatomical constraints of PSTs: the amount of synapses per length unit and area unit, depending on pre- and post-synaptic celltype, and the location along the dendritic tree of the postsynaptic cell. These are normalized using the previously mentioned PST densities.

4. Bouton densities across the entire brain area of interest.

5. Spatial distribution of cells, depending on their type.

Under the hood, multiple anatomical realizations will be computed, all of which consistent with the input data. From this distribution of anatomical realizations, the one that is closest to the average is chosen, which can be refered to as a “representative realization”. The result is saved in the same folder as the hoc morphology. This takes about 4 hrs to compute, but you can continue with a precomputed result.

To compute it yourself, run the cell below. To copy a precomputed result, you can skip to the next code cell

Computing it yourself¶

Here, we only use \(2\) samples. This is generally not enough to get a sufficiently large distribution to draw a representative sample from. In our experience, \(50\) samples is about enough to get a convergent distribution of synapse embeddings.

A sample size of \(2\) is however sufficient for the sake of this tutorial.

[ ]:

celltype = 'L5tt'  # Layer 5 thick-tufted
path = db['anatomical_constraints'].join('86_L5_CDK20041214_nr3L5B_dend_PC_neuron_transform_registered_C2center.hoc')
I.map_singlecell_inputs(
    cellName=path,
    cellTypeName=celltype,
    nrOfSamples=2  # <--- Increase this number to get more samples. 50 is sufficient for convergence, but 2 is faster
)

# Adapt the path if you have generated a new anatomical model
path_to_anatomical_model = db['anatomical_constraints'].join('example_embedding_86_C2_center')

Copying a pre-computed result¶

[7]:

# Pre-computed result: --------------------
from distutils.dir_util import copy_tree
path_to_anatomical_model = I.os.path.join(anatomy_dir, "example_embedding_86_C2_center")
silent = copy_tree(path_to_anatomical_model, db['anatomical_constraints'])

Parallelizing the generation of anatomical models¶

We can use a dask Client to create our .syn and .con files for one morphology per process:

[ ]:

client = I.get_client(timeout=10)
client

[24]:

morphology_paths = [db['anatomical_constraints'].join(f) for f in db['anatomical_constraints'].listdir() if f.endswith('.hoc')]
morphology_paths

[24]:

['/home/bgmeulem/isf_tutorial_output/network_modeling/db/anatomical_constraints/86_C2_center.hoc',
 '/home/bgmeulem/isf_tutorial_output/network_modeling/db/anatomical_constraints/86_L5_CDK20041214_nr3L5B_dend_PC_neuron_transform_registered_C2center_scaled_diameters.hoc']

This takes as long as creating an anatomical realization for a single cell, but thanks to parallellization, not much longer than that. At least, until the amount of cells exceed the amount of available threads.

[ ]:

delayed_map_singlecell_inputs = I.dask.delayed(I.map_singlecell_inputs)  # make the single cell mapper a delayed function
delayeds = [delayed_map_singlecell_inputs(p, 'L5tt') for p in morphology_paths] # call it with the morphologies
delayeds = I.dask.delayed(delayeds) # bundle everything in one delayed object
futures = client.compute(delayeds)  # compute the result

You can restart the client to interrupt this computation if don’t want to compute this right now

[ ]:

client.restart()

Inspecting the network model¶

Visualization¶

Let’s start by visualizing the locations of the somata to get a sense of our network model.

[26]:

landmarks = []
cell_types = []

somata = I.pd.read_csv(
    I.os.path.join(anatomy_dir, "presynaptic_somata_soma_locations.csv"),
    skiprows=[0,1,2,3],
    header=None,
    delimiter="\t",
    names=['Type', 'cell ID', 'x', 'y', 'z']
    )
somata[['cell_type', 'column']] = somata['Type'].str.split('_', n=1, expand=True)
somata

[26]:

[27]:

I.plt.style.use('ggplot')

Plot per cell type¶

[28]:

fig = I.plt.figure()
ax = fig.add_subplot(projection='3d')
fig.patch.set_facecolor(ax.get_facecolor())
ax.view_init(azim=40, elev=40)
ax.grid(False)

for cell_type, pts in somata.sample(2000).groupby("cell_type"):
    x, y, z = pts[["x", "y", "z"]].values.T
    ax.scatter(
        x, y, z,
        label=cell_type)

ax.legend(
    bbox_to_anchor=(1, 0.5),
    loc='center left',
    frameon=False,
    ncol=2,
    title="Cell type")
I.plt.show()

Plot per column¶

Since the barrel cortex is well-organized into recognizable columns, we can colorize the cells based on which column they belong to:

[29]:

fig = I.plt.figure()
ax = fig.add_subplot(projection='3d')
fig.patch.set_facecolor(ax.get_facecolor())
ax.view_init(azim=60, elev=40)
ax.grid(False)

for column, pts in somata.sample(2000).groupby("column"):
    x, y, z = pts[["x", "y", "z"]].values.T
    ax.scatter(
        x, y, z,
        label=column)

ax.legend(
    bbox_to_anchor=(1, 0.5),
    loc='center left',
    frameon=False,
    ncol=2,
    title="Column")
I.plt.show()

Inspecting the output format¶

In the directory generated by the singlecell_input_mapper, there are the following files:

[30]:

db['anatomical_constraints'].listdir()

[30]:

['Loader.json',
 'metadata.json',
 '86_L5_CDK20041214_nr3L5B_dend_PC_neuron_transform_registered_C2center.hoc',
 '86_L5_CDK20041214_nr3L5B_dend_PC_neuron_transform_registered_C2center_scaled_diameters.hoc',
 'example_embedding_86_C2_center.con',
 'example_embedding_86_C2_center.syn']

The most important files for making single-cell simulations are the .con file and the .syn file. These files are the relevant output of the map_singlecell_inputs pipeline.

The .con file maps presynaptic cells to a synapse. Not just a celltype, but all individual cells with a specific cell ID are mapped to an individual synapse with synase ID.

[31]:

con_file = db['anatomical_constraints'].get_file('.con')
con_file_path = db['anatomical_constraints'].join(con_file)
con_file_path

[31]:

'/home/bgmeulem/isf_tutorial_output/network_modeling/db/anatomical_constraints/example_embedding_86_C2_center.con'

[32]:

with open(con_file_path) as f:
    print(f.read()[:300])

# Anatomical connectivity realization file; only valid with synapse realization:
# 86_C2_center.syn
# Type - cell ID - synapse ID

L6cc_A3 0       0
L6cc_A3 1       1
L6cc_A3 2       2
L6cc_A3 3       3
L6cc_A3 4       4
L6cc_A3 4       5
L6cc_A3 5       6

The .syn file specifies the exact position of each synapse on the hoc morphology:

[33]:

syn_file_path = db['anatomical_constraints'].join(db['anatomical_constraints'].get_file('.syn'))
with open(syn_file_path) as f:
    print(f.read()[:300])

# Synapse distribution file
# corresponding to cell: 86_C2_center
# Type - section - section.x

VPM_E1  112     0.138046479525
VPM_E1  130     0.305058053119
VPM_E1  130     0.190509288017
VPM_E1  9       0.368760777084
VPM_E1  110     0.0
VPM_E1  11      0.120662910562

Here, section referes to the ID of the section in the cell object. x specifies, where along that section the synapse is placed. If x is 0, this is the beginning of the section, if x is one, this is the end of the section.

Compatibility with other dense connectome models.¶

ISF is not limited to this particular way of generating synapse locations. As long as a .syn and .con file can be distilled, any dense connectome model is compatible with ISF. All further workflows in ISF simply rely on these .syn and .con files, not on the input data that was necessary to generate these. in other words, the barrel cortex model as downloaded in the beginning of this tutorial is simply one example of a network model. Many other dense connectome models may be used.

Recap¶

This tutorial provided a full pipeline on how to generate synapse locations for a given neuron morphology (i.e. an “anatomical network embedding”), and how to parallellize these for multiple morphologies.

In order to assign activity to these synapse realizations, consult the the next tutorial.