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

[1]:

from os.path import exists as path_exists, join as path_join, expandvars
assert path_exists(expandvars(path_join("$ISF_HOME", "barrel_cortex"))), "Please download the barrel cortex model first by running pixi r download_bc_model"

To get started, adapt the desired output directory below:

[2]:

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

[3]:

%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/docs+0.g07c756230.dirty
[INFO] ISF: Current pid: 56117
[INFO] l5pt: Loaded mechanisms in NEURON namespace.

Warning: no DISPLAY environment variable.
--No graphics will be displayed.

[ATTENTION] ISF: The source folder has uncommited changes!



[INFO] ISF: Loaded modules with __version__ attribute are:
IPython: 8.12.2, Interface: heads/docs+0.g07c756230.dirty, PIL: 9.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.15, 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.4.0, 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: 25.0, pandas: 1.1.3, 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, 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, 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: 6.0.2, 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)

[4]:

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

[4]:

['presynaptic_somata_soma_locations.csv',
 '86_C2_center.hoc',
 '89_C2.hoc',
 '86_C2_center.swc',
 '86_C2_center_scaled_diameters.hoc',
 'example_embedding_86_C2_center']

[5]:

path_to_hoc = I.os.path.join(
    anatomy_dir,
    '86_C2_center.hoc')
path_to_scaled_hoc = I.os.path.join(
    anatomy_dir,
    '86_C2_center_scaled_diameters.hoc')

assert I.os.path.exists(path_to_hoc)
assert I.os.path.exists(path_to_scaled_hoc)

We copy these morphology files to our DataBase.

[6]:

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'])

[6]:

'/gpfs/soma_fs/home/meulemeester/isf_tutorial_output/network_modeling/db/anatomical_constraints/86_C2_center_scaled_diameters.hoc'

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:

The neuron morphology (the .hoc file that we already copied over)
Density of total Post-Synaptic Target sites (PST) across the dendritic tree of the postynaptic cell (used for normalization)
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.
Bouton densities across the entire brain area of interest.
Spatial distribution of cells, depending on their type.

For a morphology .hoc file, a network embedding is fully defined by a .syn and .con file. The former to specify the locations fo the synapses onto the morphology, and the latter to specify their presynaptic origin.

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 very long to compute (up to 4 hours), but you can continue with a precomputed result.

Computing it yourself¶

Here, we only use \(1\) sample. This is 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 \(1\) is however sufficient for the sake of this tutorial.

[7]:

celltype = 'L5tt'  # Layer 5 thick-tufted
path = db['anatomical_constraints'].join('86_C2_center.hoc')

I.logger.setLevel("INFO") # verbose output
I.map_singlecell_inputs(
    cellName=path,
    cellTypeName=celltype,
    nrOfSamples=1  # <--- Increase this number to get more samples. 50 is sufficient for convergence, but 2 is faster
)
I.logger.setLevel("WARNING")

[INFO] map_singlecell_inputs: Loading cell morphology
[INFO] ISF: Reading hoc file: /gpfs/soma_fs/home/meulemeester/isf_tutorial_output/network_modeling/db/anatomical_constraints/86_C2_center.hoc
[INFO] map_singlecell_inputs: Loading spreadsheets and bouton/PST densities...
[INFO] map_singlecell_inputs:     Loading numberOfCells spreadsheet /gpfs/soma_fs/scratch/meulemeester/project_src/in_silico_framework/barrel_cortex/nrCells.csv
[INFO] map_singlecell_inputs: Loading bouton densities from folder /gpfs/soma_fs/scratch/meulemeester/project_src/in_silico_framework/barrel_cortex/singleaxon_boutons_ascii (may take a while)
[INFO] map_singlecell_inputs: Creating network embedding for  /gpfs/soma_fs/home/meulemeester/isf_tutorial_output/network_modeling/db/anatomical_constraints/86_C2_center.hoc
[INFO] network_embedding: ---------------------------
[INFO] network_embedding: Created 536162 presynaptic cells in total
[INFO] network_embedding: ---------------------------
[INFO] network_embedding: Computing synapse densities
[INFO] synapse_mapper: ---------------------------
[INFO] network_embedding: Generating network embedding sample 0 of 1
[INFO] network_embedding: ---------------------------
[INFO] network_embedding: Creating anatomical connectivity map for output...
[INFO] network_embedding: ---------------------------
[INFO] network_embedding: ---------------------------
[INFO] network_embedding: Calculating results summary
[INFO] network_embedding:     Computing path length to soma for all synapses...
[INFO] network_embedding: ---------------------------
[INFO] network_embedding: Computing parameter distribution for 1 samples in population...
[INFO] network_embedding: ---------------------------
[INFO] network_embedding: Found representative sample with all parameters within +-2 SD
[INFO] network_embedding: Closest sample within +-2 SD (ID 0) has minimum distance 0.0 ...
[INFO] network_embedding: ---------------------------
[INFO] network_embedding: Writing output files for representative sample (id=0)
[INFO] network_embedding: ---------------------------
[INFO] network_embedding: Writing output files...

Directory Name is  /gpfs/soma_fs/home/meulemeester/isf_tutorial_output/network_modeling/db/anatomical_constraints/86_C2_center_synapses_20250821-1543_56117/
CSV file name is  /gpfs/soma_fs/home/meulemeester/isf_tutorial_output/network_modeling/db/anatomical_constraints/86_C2_center_synapses_20250821-1543_56117/86_C2_center_summary_20250821-1543_56117

---------------------------
[INFO] network_embedding: Done generating network embedding!
[INFO] network_embedding: ---------------------------
[INFO] map_singlecell_inputs: Runtime: 9.4 minutes

[23]:

# Adapt the path if you have generated a new anatomical model
path_to_anatomical_model = sorted(I.glob.glob(db['anatomical_constraints'].join('*_center_synapses*')))[-1]
print("Saved embedding results to: ", path_to_anatomical_model)

Saved embedding results to:  /gpfs/soma_fs/home/meulemeester/isf_tutorial_output/network_modeling/db/anatomical_constraints/86_C2_center_synapses_20250821-1543_56117

Copying a pre-computed result instead¶

[ ]:

# 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'].join("example_embedding_86_C2_center"))
path_to_anatomical_model = db['anatomical_constraints'].join("example_embedding_86_C2_center")

Parallelizing the generation of anatomical models¶

We can use a dask Client to create our .syn and .con files for one morphology per process. 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.

[29]:

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

client = I.get_client(timeout=10)
client

[29]:

Client

Scheduler: tcp://10.100.2.201:38786
Dashboard: http://10.100.2.201:38787/status

Cluster

Workers: 40
Cores: 40
Memory: 4.00 TB

[ ]:

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()

Let’s move the results of our synapse embedding to the database for convenience during this tutorial

[34]:

# Let's first unpack the re
I.shutil.copytree(path_to_anatomical_model, I.os.path.dirname(path_to_anatomical_model), dirs_exist_ok=True)
I.shutil.rmtree(path_to_anatomical_model)

Inspecting the network model¶

Visualization¶

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

[8]:

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

[8]:

	Type	cell ID	x	y	z	cell_type	column
0	L45Peak_D1	0	-626.213	20.344	-58.270	L45Peak	D1
1	L45Peak_D1	1	-320.834	13.905	-248.013	L45Peak	D1
2	L45Peak_D1	2	-471.762	100.233	-114.434	L45Peak	D1
3	L45Peak_D1	3	-423.770	94.010	-156.002	L45Peak	D1
4	L45Peak_D1	4	-294.765	145.499	-241.157	L45Peak	D1
...	...	...	...	...	...	...	...
19699	L34_D3	292	485.718	83.024	122.828	L34	D3
19700	L34_D3	293	483.525	-213.608	190.240	L34	D3
19701	L34_D3	294	586.649	-160.437	98.395	L34	D3
19702	L34_D3	295	275.697	189.550	319.151	L34	D3
19703	L34_D3	296	588.751	-7.494	253.711	L34	D3

19704 rows × 7 columns

[9]:

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

Plot per cell type¶

[10]:

%matplotlib inline

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()

../../_images/tutorials_2._network_models_2.1_Anatomical_embedding_34_0.png

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:

[11]:

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()

../../_images/tutorials_2._network_models_2.1_Anatomical_embedding_36_0.png

Inspecting the output format¶

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

[24]:

db['anatomical_constraints'].listdir()

[24]:

['Loader.json',
 '86_C2_center_synapses_20250821-1543_56117',
 'NumberOfConnectedCells.csv',
 '86_C2_center.hoc',
 'example_embedding_C2_center',
 'metadata.json',
 '86_C2_center_scaled_diameters.hoc']

The most important files for making single-cell simulations are the .syn and .con files. 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.

[35]:

con_file = db['anatomical_constraints'].get_file('.con')
con_file_path = db['anatomical_constraints'].join(con_file)
with open(con_file_path) as f:
    print(f.read()[:300])

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

L2_Alpha        0       0
L2_Alpha        1       1
L2_Alpha        2       2
L2_Alpha        3       3
L2_Alpha        4       4
L2_Alpha        5       5
L2_Alpha        6       6
L2_Alpha        7       7
L2_Alpha        8       8
L2_Alpha        9       9
L2_Alpha        1

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

[36]:

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.hoc
# Type - section - section.x

L2_Alpha        24      0.2942857620525994
L2_Alpha        25      0.3212573404349538
L2_Alpha        25      0.2742595025486565
L2_Alpha        23      0.012963438177888158
L2_Alpha        27      0.5271702850374373
L2_Alpha        37      0.31520261960770024
L2_Alpha        24

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 is present (no matter which pipeline created them), any dense connectome model is compatible with ISF. All further workflows in ISF simply rely on these .syn and .con files, but 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.