Data requirements for Udvary et al. (2022)¶

[78]:

import Interface as I
from pathlib import Path
from getting_started import example_data_dir
import os

You can download the barrel cortex model from Harvard Dataverse:

pixi run download_bc_model

[79]:

# Check if the model is downloaded
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"

[54]:

def print_head(fn, head=10):
    """Convenience function to print the first n lines of a file"""
    with open(fn) as input_file:
        head = [next(input_file) for _ in range(head)]
    print("".join(head))

The methods in Udvary et al. (2022) require the following empirical data:

• Amount of cells per cell type across the neuropil

• 3D density distribution of synapses across the neuropil

• Density of postsynaptic targets on the postsynaptic cell type, in units of length and area

• Locations of boutons in the neuropil, per cell type

The files that contain this data can be defined as default settings in the user configuration. Below, we inspect the content of each file to familiarize the user with the format of these data.

Amount of cells per cell type across the neuropil¶

The file format is simply a csv file:

[55]:

bc_data_dir = Path(example_data_dir).parent.parent / "barrel_cortex"
n_cells_fn = bc_data_dir / "nrCells.csv"

[56]:

n_cells_df = I.pd.read_csv(n_cells_fn, sep="\t", index_col=0)
n_cells_df

[56]:

We see how the amount of cells across the neuropil is defined for each celltype (see the "CELL TYPE" index), and for each anatomical area. For the barrel cortex, that corresponds to individual columns, index from \([A-E][1-4]\) and \([\alpha-\delta]\).

[57]:

n_cells = n_cells_df.sum().sum()
print(f"Counted {n_cells} cells in total")

Counted 536162 cells in total

The counting strategy for the data outlined above is described in Oberlaender et al. (2009)

Let’s visualize the distribution of these cells across the barrel cortex:

[58]:

%matplotlib inline
I.plt.style.use("fivethirtyeight")

def heatmap_cell_count(n_cells_df, cell_type):
    fig = I.plt.figure()
    n_cells_this_type = n_cells_df.loc[cell_type].drop(
        ["Alpha", "Beta", "Gamma", "Delta"]
    ).pipe(  # unravel index "A1-E4" into index A-E and columns 1-4
        lambda s: s.set_axis(
            I.pd.MultiIndex.from_arrays(arrays=[s.index.str[0], s.index.str[1].astype(int)])
        )
    ).unstack()

    I.sns.heatmap(
        data=n_cells_this_type,
    )
    I.plt.title(
        "Cell count for cell type: {}".format(cell_type)
    )
    I.plt.xlabel("Row")
    I.plt.ylabel("Column")
    return fig

[59]:

heatmap_cell_count(n_cells_df=n_cells_df, cell_type="L6cc");

3D density distribution of synapses across the neuropil¶

Let’s check the file format again:

[60]:

from singlecell_input_mapper.udvary2022.reader import read_scalar_field

pst_fn = str(bc_data_dir / "PST" / "EXNormalizationPSTs.am")
print_head(pst_fn, 25)

# AmiraMesh 3D ASCII 2.0


define Lattice 139 120 75

Parameters {
    Content "139x120x75 float, uniform coordinates",
    Spacing 50 50 50,
    QueryString "PSTWithCellTypesInSelection(neuron,XYZ,50.0000)",
    BoundingBox -3975 2925 -1725 4225 -2725 975,
    CoordType "uniform"
}

Lattice { float Data } @1

# Data section follows
@1
0.000000000000000e+00
0.000000000000000e+00
0.000000000000000e+00
0.000000000000000e+00
0.000000000000000e+00
0.000000000000000e+00
0.000000000000000e+00
0.000000000000000e+00

PST density fields are defined as AMIRA lattice files. These is simply a VTK-like definition of a 3D histogram. You can read this in ISF as follows:

[61]:

sf = read_scalar_field(fname=pst_fn)
mesh = sf.mesh
mesh.shape

[61]:

(139, 120, 75)

You could visualize this 3D density field in matplotlib. However, it it better suited to use dedicated 3D visualization software, as matplotlibs 3D capabilities are somewhat limited. For example, you can export this mesh to VTK to visualize it in e.g. ParaView:

[62]:

from visualize.vtk import convert_amira_lattice_to_vtk

convert_amira_lattice_to_vtk(
    pst_fn,
    outname="exc_pst"
)

After some processing, you can inspect the 3D density as 3D isosurfaces:

Density of postsynaptic targets, in units of dendrite length or membrane area¶

This is again a simple csv file:

[63]:

conn_df = I.pd.read_csv(
    bc_data_dir / "ConnectionsV8.csv",
    sep="\t"
)
conn_df.head(12)

[63]:

We can see that e.g. excitatory neurons range from \(1.04\) to \(1.68\) spines per \(\mu m\) dendritic length. For example:

[64]:

post_ct = "L2"
conn_df[conn_df['POSTSYNAPTIC_CELLTYPE'] == post_ct]

[64]:

This means that we can expect about \(1.68\) excitatory PSTs per \(\mu m\) of L2 apical dendrite, and \(0.03\) inhibitory PSTs per \(\mu m^2\) of L2 soma membrane. More info for how these values were calculated can be found in Udvary et al. (2022)

Locations of boutons in the neuropil, per cell type¶

[65]:

boutons_dir = bc_data_dir / "singleaxon_boutons_ascii"
I.os.listdir(boutons_dir)

[65]:

['C2',
 'E1',
 'E4',
 'B2',
 'D1',
 'D4',
 'A4',
 'Delta',
 'A2',
 'C1',
 'C4',
 'Beta',
 'E3',
 'B1',
 'B4',
 'Alpha',
 'D3',
 'A1',
 'C3',
 'Gamma',
 'E2',
 'B3',
 'D2',
 'A3']

We see that the locations of boutons are organized per anatomical area. For the barrel cortex, that means per column. Let’s inspect a single column.

[66]:

boutons_c2_dir = boutons_dir / "C2"
I.os.listdir(boutons_c2_dir)

[66]:

['SymLocal1',
 'L45Sym',
 'L34',
 'SymLocal2',
 'L45Peak',
 'L5st',
 'L6cc',
 'L5tt',
 'L6ccinv',
 'SymLocal6',
 'L4ss',
 'SymLocal3',
 'L1',
 'VPM',
 'SymLocal4',
 'L4py',
 'L6ct',
 'L23Trans',
 'L4sp',
 'L2',
 'L56Trans',
 'SymLocal5']

Here, they are again organized per cell type.

[67]:

cell_type = "L6cc"
ct_dir = boutons_c2_dir / cell_type
I.os.listdir(ct_dir)

[67]:

['L6ccaxon1boutonsC2.am',
 'L6ccaxon4boutonsC2.am',
 'L6ccaxon2boutonsC2.am',
 'L6ccaxon0boutonsC2.am',
 'L6ccaxon3boutonsC2.am',
 'L6ccaxon5boutonsC2.am']

Finally, we see the actual data files. Let’s inspect one of them:

[68]:

boutons_fn = [e for e in ct_dir.glob("*.am")][0]
print_head(boutons_fn, 25)

# AmiraMesh 3D ASCII 2.0


define Lattice 39 34 42

Parameters {
    Content "39x34x42 float, uniform coordinates",
    Spacing 50 50 50,
    QueryString "boutons(axon WHERE (CellTypeIs('L6ccaxon1') AND PrimaryColumnIs('C2')),XYZ,50.0000)",
    BoundingBox -925 975 -375 1275 -1475 575,
    CoordType "uniform"
}

Lattice { float Data } @1

# Data section follows
@1
0.000000000000000e+00
0.000000000000000e+00
0.000000000000000e+00
0.000000000000000e+00
0.000000000000000e+00
0.000000000000000e+00
0.000000000000000e+00
0.000000000000000e+00

This is again an AMIRA mesh file. We can visualize like before:

[69]:

convert_amira_lattice_to_vtk(
    boutons_fn,
    "boutons"
)

This shows a bouton density for a cell type, located at some anatomical area in the neuropil (in this case: a L6cc in the C2 column).

We can also project onto an axis:

[75]:

sf = read_scalar_field(str(boutons_fn))

proj_axis = 1
projection = I.np.sum(
    sf.mesh,
    axis=proj_axis
)

I.sns.heatmap(
    I.np.flip(projection.T, 0)
)
I.plt.xlabel("X")
I.plt.ylabel("Z")

[75]:

Text(26.78333333333333, 0.5, 'Z')

Additional data: representative soma locations¶

We can visualize a sample of these soma locations to get a sense of our network model.

[71]:

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

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

[71]:

Plot per cell type¶

[73]:

%matplotlib inline
I.plt.style.use('ggplot')

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:

[74]:

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