Simultaneously compute edge location and count onto biophysical postsynaptic neurons

[nicolomeneghetti]

Dear all,

I was trying to compute the distribution of synapses' location along the z-axis onto a family of post-synaptic neurons.

As far as I've understood, to do so, you are supposed to work with the 'add_properties' of the NetworkBuilder by doing something like this:

`conns = net.add_edges(
source={'type': 'basket'},
target={'type': 'pyramidal', 'model_type': 'biophysical'},
connection_rule=16,
syn_weight=0.00025,
delay=2.0,
dynamics_params='GABA_InhToExc.json',
model_template='Exp2Syn',
)

conns.add_properties(
['afferent_section_id', 'afferent_section_pos', 'afferent_swc_id', 'afferent_swc_pos'],
rule=rand_syn_locations,
rule_params={
'sections': ['basal', 'apical', 'axon'],
'distance_range': [30.0, 100.0],
'morphology_dir': 'components/morphologies'
},
dtypes=[np.int, np.float, np.int, np.float]
)`

Then you can work with the _edges.h5 file and the .swc morphology file to get the synapses’ location along the desired axis. So far (I guess) so good.

1. The first question is whether there is a simpler workaround to get the location of the synapses. In the bmtk.builder package page there is something that seems to resemble what I was looking for. Unfortunately the information there seems to be missing:

From: https://alleninstitute.github.io/bmtk/bmtk/bmtk.builder.html#:~:text=.float)-,We,-can%20also%20add
`def get_target_loc_and_dist(source, target):
...
return trg_loc, trg_dist

cm.add_properties(['syn_location', 'syn_distance'],
rule=get_target_loc_and_dist,
dtypes=[string, np.float])`

As you can see the get_target_loc_and_dist function is left empty.

2. The second question is more worrisome for me. When I run the add_edges of the bmtk builder without the add_properties:

net.add_edges( source={'type': 'basket'}, target={'type': 'pyramidal', 'model_type': 'biophysical'}, connection_rule=16, target_sections=['basal', 'apical', 'axon'], distance_range=[30.0, 100.0], .... )

I get into the resulting _edges.h5 file a key field called 'nsyns'. As far as I have understood the 'nsyns' contains the number of connections from source to target for each edge (am I wrong?). Unfortunately, however, if I run the code from above (the one with the 'add_properties') the 'nsyns' disappear. So my question is: can I get both the information? Specifically, am I able to get both the location of each edge and its numerosity (that if I stand correct is contained within the 'nsyns' field)?

Thanks in advance for any help.

Best regards,

Nicolò

[xpliu16]

Hi Nicolo,

Regarding 1. Unfortunately, get_target_loc_and_dist is just pseudocode for the user defined custom connectivity rule that is passed to add_properties to generate the synapse locations. Here, we are trying to illustrate that you can assign multiple properties simultaneously.

Regarding 2. `This is the correct behavior, but it is confusing. The implied nsyns after add_properties is 1, which is the default if no nsyns is specified. Since add_properties gives a unique property to each edge, each edge is its own entry in the table. So for instance, if previously we had 1 row representing synapses from neuron A to neuron B with nsyns 5, after add_properties, there would instead be 5 lines representing each synapse, with their individual locations.

In terms of being able to help with your overall goal, do you want to see the locations as relative to the neuron or in absolute space coordinates?

Thanks,
Ping

[nicolomeneghetti]

Dear Ping,

thanks first of all for your reply. Now the 'nsyns' actually make sense to me :)

My overall goal is actually to understand how bmtk place the synapses within a user-specified distance_range. Does it follow a uniform distribution (which is my guess)? Or a truncated normal? Or another one?

To try answering this question my goal is therefore to compute the distribution of the location in absolute space coordinates of all the synapses belonging to a particular edge type (say, for example, the distribution in absolute space coordinates of excitatory to excitatory synapses).

Thanks in advance,

Best regards,

Nicolò Meneghetti

[xpliu16]

Hi Nicolo,

With the rand_syn_locations function, the synapses are distributed uniformly per unit length, and then filtered by the distance criteria on total arc length distance from the soma. Of course, one may want a different targeting rule, such as Gaussian with tapering rather than a sharp cutoff, or targeting within a spatial box, but we do not currently have ready-made functions for other targeting rules.

We have made a modification to allow rand_syn_locations() to write the x,y,z coordinates to the HDF5 file which will be available in the next release. However, you can just copy the modified rand_syn_locations() function below into your build script rather than importing it from bmtk. You need to modify the call as follows:

conns.add_properties(
    ['afferent_section_id', 'afferent_section_pos', 'afferent_swc_id', 'afferent_swc_pos', 'coord_x', 'coord_y', 'coord_z'],
    rule=rand_syn_locations,
    rule_params={
        'sections': ['soma', 'basal', 'apical'], 
        'distance_range': [0.0, 1.0e20],
        'morphology_dir': 'components/morphologies',
        'return_coords': True
    },
    dtypes=[int, float, int, float, float, float, float]

Make sure to update the dtypes list as indicated for the new coordinate properties. Below is the modified rand_syn_locations():

def rand_syn_locations(src, trg, sections=('soma', 'apical', 'basal'), distance_range=(0.0, 1.0e20),
                       morphology_dir='./components/morphologies', return_coords=False):
    trg_swc = get_swc(trg, morphology_dir=morphology_dir, use_cache=True)
    sec_ids, seg_xs = trg_swc.choose_sections(sections, distance_range, n_sections=1)
    sec_id, seg_x = sec_ids[0], seg_xs[0]
    swc_id, swc_dist = trg_swc.get_swc_id(sec_id, seg_x)
    if return_coords:
        coords = trg_swc.get_coords(sec_id, seg_x)
        return [sec_id, seg_x, swc_id, swc_dist, coords[0], coords[1], coords[2]]
    else:
        return [sec_id, seg_x, swc_id, swc_dist]

You may need to import:
from bmtk.builder.bionet import SWCReader
from bmtk.builder.bionet import get_swc

Let us know if that solution works - if not, we can try to update the conda with a new release soon.

Note that the distribution will not be uniform in x,y, or z because they are uniformly distributed per length of neuronal processes. For instance, if you have a single long vertical process and then a big tuft at the end, the tuft will have a much higher density of synapses per y distance.

Also note that following Allen database conventions, y is the vertical axis spanning the depth of cortex, not z. Cell models initial orientations are not necessarily exactly vertical (they may need to be rotated).

Hope that helps!
-Ping

[nicolomeneghetti]

Dear Ping,

thanks for your answer. It actually works perfectly. Just a warning for other users: using the rand_syn_locations function above slows down significantly the network building process.

Just an additional question, if I may, before closing this issue. As you stated in your previous comment "Cell models initial orientations are not necessarily exactly vertical (they may need to be rotated)". Is there an automatic way to get the cell morphology automatically aligned along the depth- (i.e., y-) axis? Or in general a script to get the degree of orientation needed to align the morphology along the y-axis? As of now I am doing it by hand on VND, but I think an automatic way would be beneficial for many users in my opinion.

Thanks again for your availability,

Best regards,

Nicolò

[xpliu16]

Glad to hear it is working for you. I think making a function to align arbitrary biological cell types relative to the pia based on morphology would be rather difficult (although there are folks here studying the morphology of cell types rather quantitatively). However, we think we could do something for pyramidal cells that would be an initial guess, such as aligning the first section of the apical dendrite with respect to a user specified axis (e.g., 'y'). And we could facilitate rotation within VND and provide easy read out of the values. Is this the kind of solution that would be helpful? Did you have a specific implementation in mind?

Thanks,
Ping

[nicolomeneghetti]

Dear Ping,

I actually do not any specific implementation in mind, but what you are suggesting seems a good starting point.

Thanks again for your availability,

Nicolò

[nicolomeneghetti]

Dear Ping,

sorry for bothering you again.

As mentioned before, your fix (the modified "rand_syn_locations" function) works perfectly with a small network I built from one the bmtk GitHub example. It failed, however, when applied to the complete V1 model of Billeh et al. Just to reiterate, my goal is to understand the location of synapses in this model along the depth-axis (which is the y-axis in bmtk).

To try achieving this goal I modified the build_network.py available here: https://www.dropbox.com/sh/w5u31m3hq6u2x5m/AAAfeKZJ9ArMjPxBSK4FZXLca/simulations/v1_biophysical?dl=0&preview=build_network.py&subfolder_nav_tracking=1 as so:

bio_edges = net.add_edges(
source={'pop_name': src_type},
target={'node_type_id': node_type_id},
iterator='all_to_one',
connection_rule=connect_cells,
connection_params={'params': src_trg_params},
dynamics_params=row['params_file'],
model_template='exp2syn',
syn_weight=row['weight_max'],
delay=row['delay'],
weight_function=weight_fnc,
weight_sigma=weight_sigma,
)

            # I am sure there's a smarter way to do so, yet this works: 
            bio_edges_distance_range = row['distance_range']
            bio_edges_distance_range = bio_edges_distance_range.replace('[','')
            bio_edges_distance_range = bio_edges_distance_range.replace(']','')
            bio_edges_distance_range = [float(item) for item in bio_edges_distance_range.split(',')]
            
            bio_edges_target_sections = row['target_sections']
            bio_edges_target_sections = bio_edges_target_sections.replace('[','')
            bio_edges_target_sections = bio_edges_target_sections.replace(']','')
            bio_edges_target_sections = bio_edges_target_sections.replace('"','')
            bio_edges_target_sections = bio_edges_target_sections.replace(' ','')
            bio_edges_target_sections = [item for item in bio_edges_target_sections.split(',')]
            
            bio_edges.add_properties(
                ['afferent_section_id', 'afferent_section_pos', 'afferent_swc_id', 'afferent_swc_pos', 'coord_x', 'coord_y', 'coord_z'],
                # ['coord_x', 'coord_y', 'coord_z'],
                rule=rand_syn_locations,
                rule_params={
                    'sections': bio_edges_target_sections, 
                    'distance_range': bio_edges_distance_range,
                    'morphology_dir': '/V1 models/components/morphologies',
                    'return_coords': True
                },
                dtypes=[int, float, int, float, float, float, float])

Unfortunately this code returns the following error:

runfile('/home/nest/Desktop/bmtk github stuff/bmtk-develop/examples/V1 models/build_network.py', wdir='/home/nest/Desktop/bmtk github stuff/bmtk-develop/examples/V1 models')
2022-12-14 13:45:05,553 [WARNING] Adjusted calculated probability based on distance dependence is coming out to be greater than 1 for i23Sst and LIFi23Pvalb. Setting to 1.0
2022-12-14 13:45:05,555 [WARNING] Adjusted calculated probability based on distance dependence is coming out to be greater than 1 for LIFi23Sst and LIFi23Pvalb. Setting to 1.0
2022-12-14 13:45:05,563 [WARNING] Adjusted calculated probability based on distance dependence is coming out to be greater than 1 for i23Sst and i23Pvalb. Setting to 1.0
2022-12-14 13:45:05,563 [WARNING] Adjusted calculated probability based on distance dependence is coming out to be greater than 1 for i23Sst and i23Pvalb. Setting to 1.0
2022-12-14 13:45:05,564 [WARNING] Adjusted calculated probability based on distance dependence is coming out to be greater than 1 for i23Sst and i23Pvalb. Setting to 1.0
2022-12-14 13:45:05,571 [WARNING] Adjusted calculated probability based on distance dependence is coming out to be greater than 1 for LIFi23Sst and i23Pvalb. Setting to 1.0
2022-12-14 13:45:05,572 [WARNING] Adjusted calculated probability based on distance dependence is coming out to be greater than 1 for LIFi23Sst and i23Pvalb. Setting to 1.0
2022-12-14 13:45:05,572 [WARNING] Adjusted calculated probability based on distance dependence is coming out to be greater than 1 for LIFi23Sst and i23Pvalb. Setting to 1.0
2022-12-14 13:45:05,741 [WARNING] Adjusted calculated probability based on distance dependence is coming out to be greater than 1 for i4Sst and LIFi4Pvalb. Setting to 1.0
2022-12-14 13:45:05,749 [WARNING] Adjusted calculated probability based on distance dependence is coming out to be greater than 1 for i5Sst and i5Pvalb. Setting to 1.0
2022-12-14 13:45:05,749 [WARNING] Adjusted calculated probability based on distance dependence is coming out to be greater than 1 for i5Sst and i5Pvalb. Setting to 1.0
2022-12-14 13:45:05,750 [WARNING] Adjusted calculated probability based on distance dependence is coming out to be greater than 1 for i5Sst and i5Pvalb. Setting to 1.0
2022-12-14 13:45:05,751 [WARNING] Adjusted calculated probability based on distance dependence is coming out to be greater than 1 for i5Sst and i5Pvalb. Setting to 1.0
2022-12-14 13:45:05,751 [WARNING] Adjusted calculated probability based on distance dependence is coming out to be greater than 1 for i5Sst and i5Pvalb. Setting to 1.0
2022-12-14 13:45:05,752 [WARNING] Adjusted calculated probability based on distance dependence is coming out to be greater than 1 for i5Sst and i5Pvalb. Setting to 1.0
2022-12-14 13:45:05,752 [WARNING] Adjusted calculated probability based on distance dependence is coming out to be greater than 1 for i5Sst and i5Pvalb. Setting to 1.0
2022-12-14 13:45:05,753 [WARNING] Adjusted calculated probability based on distance dependence is coming out to be greater than 1 for i5Sst and i5Pvalb. Setting to 1.0
2022-12-14 13:45:05,753 [WARNING] Adjusted calculated probability based on distance dependence is coming out to be greater than 1 for i5Sst and i5Pvalb. Setting to 1.0
2022-12-14 13:45:05,754 [WARNING] Adjusted calculated probability based on distance dependence is coming out to be greater than 1 for i5Sst and LIFi5Pvalb. Setting to 1.0
2022-12-14 13:45:05,764 [WARNING] Adjusted calculated probability based on distance dependence is coming out to be greater than 1 for LIFi4Sst and LIFi4Pvalb. Setting to 1.0
2022-12-14 13:45:05,770 [WARNING] Adjusted calculated probability based on distance dependence is coming out to be greater than 1 for LIFi5Sst and i5Pvalb. Setting to 1.0
2022-12-14 13:45:05,770 [WARNING] Adjusted calculated probability based on distance dependence is coming out to be greater than 1 for LIFi5Sst and i5Pvalb. Setting to 1.0
2022-12-14 13:45:05,771 [WARNING] Adjusted calculated probability based on distance dependence is coming out to be greater than 1 for LIFi5Sst and i5Pvalb. Setting to 1.0
2022-12-14 13:45:05,772 [WARNING] Adjusted calculated probability based on distance dependence is coming out to be greater than 1 for LIFi5Sst and i5Pvalb. Setting to 1.0
2022-12-14 13:45:05,772 [WARNING] Adjusted calculated probability based on distance dependence is coming out to be greater than 1 for LIFi5Sst and i5Pvalb. Setting to 1.0
2022-12-14 13:45:05,773 [WARNING] Adjusted calculated probability based on distance dependence is coming out to be greater than 1 for LIFi5Sst and i5Pvalb. Setting to 1.0
2022-12-14 13:45:05,773 [WARNING] Adjusted calculated probability based on distance dependence is coming out to be greater than 1 for LIFi5Sst and i5Pvalb. Setting to 1.0
2022-12-14 13:45:05,774 [WARNING] Adjusted calculated probability based on distance dependence is coming out to be greater than 1 for LIFi5Sst and i5Pvalb. Setting to 1.0
2022-12-14 13:45:05,774 [WARNING] Adjusted calculated probability based on distance dependence is coming out to be greater than 1 for LIFi5Sst and LIFi5Pvalb. Setting to 1.0
2022-12-14 13:45:05,783 [WARNING] Adjusted calculated probability based on distance dependence is coming out to be greater than 1 for i4Sst and i4Pvalb. Setting to 1.0
2022-12-14 13:45:05,784 [WARNING] Adjusted calculated probability based on distance dependence is coming out to be greater than 1 for i4Sst and i4Pvalb. Setting to 1.0
2022-12-14 13:45:05,788 [WARNING] Adjusted calculated probability based on distance dependence is coming out to be greater than 1 for LIFi4Sst and i4Pvalb. Setting to 1.0
2022-12-14 13:45:05,788 [WARNING] Adjusted calculated probability based on distance dependence is coming out to be greater than 1 for LIFi4Sst and i4Pvalb. Setting to 1.0

One point section Import3d_Section[56064] ending at line 1141 has been removed
One point section Import3d_Section[56293] ending at line 1141 has been removed
NEURON: Arg out of range in user function
near line 0
{create axon[2]}
^
List[2349].object(228)
Traceback (most recent call last):

File "/home/nest/Desktop/bmtk github stuff/bmtk-develop/examples/V1 models/build_network.py", line 374, in
v1.build()

File "/home/nest/anaconda3/envs/allen_inst_bkmt/lib/python3.9/site-packages/bmtk/builder/network_builder.py", line 284, in build
self.adaptor.build(force=force)

File "/home/nest/anaconda3/envs/allen_inst_bkmt/lib/python3.9/site-packages/bmtk/builder/network_adaptors/network.py", line 538, in build
self.__build_edges()

File "/home/nest/anaconda3/envs/allen_inst_bkmt/lib/python3.9/site-packages/bmtk/builder/network_adaptors/network.py", line 519, in __build_edges
self._add_edges(conn_map, i)

File "/home/nest/anaconda3/envs/allen_inst_bkmt/lib/python3.9/site-packages/bmtk/builder/network_adaptors/dm_network.py", line 202, in _add_edges
pvals = rule(source_node, target_node)

File "/home/nest/Desktop/bmtk github stuff/bmtk-develop/examples/V1 models/build_network.py", line 26, in rand_syn_locations
swc_id, swc_dist = trg_swc.get_swc_id(sec_id, seg_x)

File "/home/nest/anaconda3/envs/allen_inst_bkmt/lib/python3.9/site-packages/bmtk/builder/bionet/swc_reader.py", line 404, in get_swc_id
filtered_swc = self.swc_map[(self.swc_map['type'] == sec_type) & (self.swc_map['nameindex'] == sec_nameindex)]

File "/home/nest/anaconda3/envs/allen_inst_bkmt/lib/python3.9/site-packages/bmtk/builder/bionet/swc_reader.py", line 227, in swc_map
name_index = int(self.hobj.nl.sections.object(sec_id).nameindex)

RuntimeError: hoc error

Any help on this regard would be greatly appreciated.

Thanks in advance,

Best regards,

Nicolò Meneghetti

[xpliu16]

Hi Nicolo,
I am able to replicate your error and have a sense of where in the code it might be coming from. But I'd like to meet with Kael after the holidays to get a better sense of the bigger context of why and to provide you with a solution.
Thanks,
Ping

[xpliu16]

Hi Nicolo,

The issue seems to relate to a NEURON function that removes anomalous sections with zero length after conversion from the SWC file and subsequent out-of-syncness of the internal NEURON element that we are using to map from SWC points to NEURON sections. We will investigate this further.

In the V1 model, there are a few morphology files with this anomaly, leading to the program throwing an error.

The good news is that this only affects the swc_map which is used to get the swc_id and swc_dist for each synapse in the h5 table. You do not need this to get the 3D coords of the synapses. So the easiest way to bypass this is to modify your rand_syn_locations as follows:

def rand_syn_locations(src, trg, sections=('soma', 'apical', 'basal'), distance_range=(0.0, 1.0e20),
                       morphology_dir='./components/morphologies', return_coords=False, dL=None):
    trg_swc = get_swc(trg, morphology_dir=morphology_dir, use_cache=True, dL=dL)

    sec_ids, seg_xs = trg_swc.choose_sections(sections, distance_range, n_sections=1)
    sec_id, seg_x = sec_ids[0], seg_xs[0]
    #swc_id, swc_dist = trg_swc.get_swc_id(sec_id, seg_x)
    swc_id = -1
    swc_dist = -1
    if return_coords:
        coords = trg_swc.get_coords(sec_id, seg_x)
        return [sec_id, seg_x, swc_id, swc_dist, coords[0], coords[1], coords[2]]
    else:
        return [sec_id, seg_x, swc_id, swc_dist]

Depending on your version of BMTK, you may still get a call to swc_map, in which case, please try with the latest Github version.

Note that based on my understanding of your modification, this is not the normal flow of the build code, but rather a modification to read out the synapse locations. So be careful when interpreting the edges h5 tables. For instance, if you applied this to V1_V1_edges and only the morphological (non-LIF) models, the afferent_section_id and coordinate fields will only apply to the morphological models while nsyns applies only to the LIF models.

Let us know if this does what you need.

[xpliu16]

Hi Nicolo,

I wanted to clarify some things: these issues should not affect the simulation using the standard V1 model building procedure. In that case, the synapse locations are not written to the SONATA file and are placed at runtime instead, which is part of why it's been hard to read them out! We started trying to create explicit placements for the use case of visualization of synapses in VND, and this turns out to be useful for situations like yours as well, where you need to read out the synapse placements. But that process is being ironed out, and thanks for reporting these issues.

Would it be alright if we closed this issue for now?

-Ping

[nicolomeneghetti]

Dear Ping,

sorry for the delay in the response. I have a very last question before closing the issue for good.

I tried, following your suggestion, to bypass the computation of the swc_map. For the sake of computational need, I build the network by reducing the overall numerosity of the network (which counts only the 4% of the original Allen model). In the returned "v1_v1_edges.h5" there's something I am not sure about.

Specifically in this h5 file I can find:

1. v1_h5_file['edges']['v1_to_v1']['source_node_id'] : which has a dimension of 238529
2. v1_h5_file['edges']['v1_to_v1']['0']['nsyns'] : which has a dimension of 95135
3. v1_h5_file['edges']['v1_to_v1']['1']['coord_y']: which has a dimension of 143394

I am not sure (1) the reason why source_node_id has a higher dimension of coord_y; (2) the reason of this overall dimensionality differences; (3) why it is included in the files the nsyns field in the first place. I did not indeed expect to find the nsyns, as (with the code of rand_syn_locations above) I was expecting that ALL the edges (despite their dimension in the nsyns) to be listed singularly in the coord_y field; (4) are all the edges listed individually in the coord_y field or should I weight their numerosity according to the nsyns field?

Thanks in advance for your costant support,

Best,

Nicolò

[xpliu16]

Yes, this surprised me too initially. As I said above, "For instance, if you applied this to V1_V1_edges and only the morphological (non-LIF) models, the afferent_section_id and coordinate fields will only apply to the morphological models while nsyns applies only to the LIF models." Because the model has morphological neurons in the core of the column and simplified LIF neurons with no morphology in the surround part of the column, if you apply rand_syn_locations in the build code only to the morphological neurons (which makes sense, since synapse locations would be meaningless for neurons without morphology), nsyns applies only to the LIF neurons, and coords only apply to the morphological neurons. Then it makes sense that 95135 + 143394 = 238529 total synapse entries (the total number of synapses is still greater than that due to multiplicity of nsyns). This is not a usual way for the model to be built, but I agree perhaps each entry should have a row and the not applicable entries should have placeholders, and will discuss with Kael. SONATA storage is intended for efficiency to reduce the size of the files, but we may not have considered all use cases! For now, you can just ignore nsyns and use the coords for your purposes of computing synapse locations.

Hope that helps!

[nicolomeneghetti]

Dear Ping,

it now makes completely sense to me. Thank you very much for your help. I think we can retain this issue as closed.

Thanks again for you help,

Best,

Nicolò