NetPyNE simulation results for 300 biophys cell network

salvadord commented 6 years ago

Raster plot generated by netpyne for the 300_cells SONATA example
Simulation run output (run time etc.):

Running simulation for 1500.0 ms... Done; run time = 390.51 s; real-time ratio: 0.00.

Gathering data... Done; gather time = 6.61 s.

Analyzing... Cells: 400 Connections: 45245 (113.11 per cell) Synaptic contacts: 148320 (370.80 per cell) Spikes: 26026 (43.38 Hz) Simulated time: 1.5 s; 1 workers Run time: 390.51 s

Potential reasons of why the raster looks different: -- Incorrect mapping of sec_id to sections -- Incorrect mapping of population cell ids to gids -- Incorrect importing of the cell dynamic params (channel properties, distribution, etc)
NetPyNE code to import SONATA format here
NetPyNE example to import the 300 cell model, run the simulation and plot the raster here

wvangeit commented 6 years ago

Are you replacing the axon with an AIS stub ? (as far as I remember this is necessary in the 300 cell network). I don't see it in your code

salvadord commented 6 years ago

I'm just reading the SONATA saved files -- is this axon replacement something that needs to be done 'manually' in a script after loading? Do you have example code where you make this replacement?

wvangeit commented 6 years ago

yes, the descriptions of the algorithms you can find in https://github.com/AllenInstitute/sonata/blob/master/docs/SONATA_DEVELOPER_GUIDE.md

when you search for e.g. aibs_perisomatic

I think the network here uses aibs_perisomatic. Some implementation of this you can find in the bmtk: https://github.com/AllenInstitute/bmtk/blob/develop/bmtk/simulator/bionet/default_setters/cell_models.py#L82

salvadord commented 6 years ago

I added the axon replacement but I'm now the raster looks less similar to your version.

Question: the new axon sections are added to the end of the cell template all list -- before, the axons were close to the beginning of the list. This means the sec_id of most sections changes. Is this the expected behavior?

EDIT: Just saw in this other issue that the indexing of the all list after replacing was indeed causing reproducibility problems -- https://github.com/AllenInstitute/sonata/issues/15#issuecomment-403027891

salvadord commented 6 years ago

With axon replacement:
With nseg based on dL:
With both axon replacement + nseg based on dL:

salvadord commented 6 years ago

With axon replacement (using wvangeit method to preserve sec_id: https://github.com/AllenInstitute/sonata/issues/15#issuecomment-403045728) + nseg based on dL (seems like getting closer):

Cells: 400 Connections: 45245 (113.11 per cell) Synaptic contacts: 148320 (370.80 per cell) Spikes: 17900 (29.83 Hz) Simulated time: 1.5 s; 1 workers Run time: 898.95 s

wvangeit commented 6 years ago

That's good news. Could I suggest that you first try to run the 5-cell iclamp and 9-cell examples mentioned github.com/AllenInstitute/sonata/issues/15 It separates you can run into (electrical models, synapses models, synapse locations, etc).

salvadord commented 5 years ago

Raster of 300 cell example without local connectivity (only external inputs)
Voltage traces of 300 cell example without local connectivity (only external inputs)
Simulation run output log

Running simulation for 1500.0 ms... Done; run time = 801.80 s; real-time ratio: 0.00.

Gathering data... Done; gather time = 5.86 s.

Analyzing... Cells: 400 Connections: 24032 (60.08 per cell) Synaptic contacts: 60430 (151.07 per cell) Spikes: 6498 (10.83 Hz) Simulated time: 1.5 s; 1 workers Run time: 801.80 s

salvadord commented 5 years ago

@kaeldai - in bmtk how can I check the list of presyn inputs of a cell, and assuming these are vecstims with a list of spike times, how can I check those spike times interactively from python terminal? I want to compare directly with the netpyne version to make sure cells are getting the right inputs. thanks!

salvadord commented 5 years ago

Pasting here BBP results for reference:

salvadord commented 5 years ago

Latest netpyne 300_cells raster, after matching 9_cells and 9_cells_iclamp:

model_output_raster_axonv2_dl_300cells

Still some discrepancies, which could potentially come from differences in the morphology 3d points I observed in the 9 cells example. Will investigate more.

salvadord commented 5 years ago

BMTK 300 cells raster:

NetPyNE 300 cells raster:

Comparison BMTK/NetPyNE 300 cells raster: comparison_raster_withconns

Comparison BMTK/NetPyNE 300 cells traces: comparison_traces_withconns

salvadord commented 5 years ago

Interestingly if I remove intrinsic conns (internal->internal) and leave only external inputs (external->internal), I get a perfect match: comparison_raster comparison_traces

I checked and external connections project to different cell sections (not just soma) and they seem to target most cells. So not sure where the issue could be with the intrinsic connections -- any suggestions?

Also, @kaeldai is there a way to get the info (sec, loc, weight, delay, presyn gid) for all conns in a postsyn cell?

antonarkhipov commented 5 years ago

I think in cases like this Michael Hines recommends to look at the first spike where the two simulations disagree. That seems to happen very early on here for inhibitory neuron spiking. It may be helpful to look at the voltage traces of those cells (around GID 270-280).

On Wed, Mar 6, 2019 at 1:57 PM Salvador Dura-Bernal < notifications@github.com> wrote:

Interestingly if I remove intrinsic conns (internal->internal) and leave only external inputs (external->internal), I get a perfect match: [image: comparison_raster] https://user-images.githubusercontent.com/8440562/53916507-529adf00-4030-11e9-96da-889ea9fbb1be.png [image: comparison_traces] https://user-images.githubusercontent.com/8440562/53916508-529adf00-4030-11e9-8213-55deec6c4647.png

I checked and external connections project to different cell sections (not just soma) and they seem to target most cells. So not sure where the issue could be with the intrinsic connections -- any suggestions?

Also, @kaeldai https://github.com/kaeldai is there a way to get the info (sec, loc, weight, delay, presyn gid) for all conns in a postsyn cell?

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub https://github.com/AllenInstitute/sonata/issues/54#issuecomment-470293605, or mute the thread https://github.com/notifications/unsubscribe-auth/ANDGJ9FSZok0jiix16W92FyJ8dMO6Ezrks5vUDm2gaJpZM4XIjlH .

salvadord commented 5 years ago

Ok, will try that. Though what surprises me is that just with the extrinsic conns the results are identical, even for those inhib cells, so I think it's unlikely that the issue is due to the cell biophysics. Seems to be related to the local conns, perhaps those involving inhib cells.

wvangeit commented 5 years ago

Just a guess based on your info. Could the spike detection method be different ?

salvadord commented 5 years ago

Ok found the issue! The network includes 390 self-connections i.e. with same pre and post gid. For example, conn indices 1755-1759 are all connections from cell 9 to cell 9. In netpyne self-connections are omitted by default.

Once I allowed for self-connections in netpyne I got a much better match: comparison_raster comparison_traces

btw, should the network include self-connections?

antonarkhipov commented 5 years ago

That's really cool. Great to see the improvement.

I think self-connections do occur sometimes in real networks, so in principle there's no reason for SONATA to prohibit them. However, self-connections in our examples is a weird thing. That's probably just an artifact of a simple random connectivity algorithm. We may need to fix that before finalizing everything.

For time being, though, let's see if we can get to the bottom of the discrepancy. The latest result looks much better, but there's still mismatch between NetPyNE and BMTK. One thing I noticed is that it looks like cell with GID 300 (or is it 299?) fires spikes with NetPyNE, but not BMTK, and the first spike there is very early on. What's going on there? Or is it a plotting artifact?

On Thu, Mar 7, 2019 at 1:29 PM Salvador Dura-Bernal < notifications@github.com> wrote:

Ok found the issue! The network includes 390 self-connections i.e. with same pre and post gid. For example, conn indices 1755-1759 are all connections from cell 9 to cell 9. In netpyne self-connections are omitted by default.

Once I allowed for self-connections in netpyne I got a much better match: [image: comparison_raster] https://user-images.githubusercontent.com/8440562/53990396-31042b00-40f6-11e9-8ee5-8947f49c7064.png [image: comparison_traces] https://user-images.githubusercontent.com/8440562/53990402-3497b200-40f6-11e9-93d9-2b66acff289e.png

btw, should the network include self-connections?

— You are receiving this because you commented. Reply to this email directly, view it on GitHub https://github.com/AllenInstitute/sonata/issues/54#issuecomment-470702210, or mute the thread https://github.com/notifications/unsubscribe-auth/ANDGJwWOaDfBTfpTd5TPoPGgV6jCxvXKks5vUYTRgaJpZM4XIjlH .

salvadord commented 5 years ago

I noticed the gid 300 spikes, but it's just spikes from the external inputs which I shouldn't be plotting (should go up to gid 299 not 300) -- I already uploaded the fixed figure.

At this point, the biophysics (based on h.psection(sec)) and the connections (based on h.List('NetCon')) are identical in both models. So the only remaining options are differences in morphology or some global NEURON options (eg. cvode.cache_efficient etc). I will try to check these before the next meeting.

salvadord commented 5 years ago

The discrepancy in the morphology 3d points between netpyne and bmtk was because netpyne adds the x,y,z location of the cell to the 3d points. This means if you used e.g. the NEURON GUI (Interviews) to represent all the cells they would show in the correct locations; whereas for in the bmtk representation they would all show on top of each other.

However, this should not affect the cell output since the relative positions of cell sections remains identical in both cases. I removed this feature from netpyne (such that 3d points don't include the cell location) and compared the .psection() output, rasters and traces to bmtk. For some reason, this reduced the spike time (raster) and voltage trace discrepancies but there is still some remaining.

Note: I also corrected the orientation of the soma, which I had flipped for visualization purposes.

comparison_raster

comparison_traces_270-280

salvadord commented 5 years ago

Also, noticed that bmtk does not recenter 3d points so the soma corresponds to the cell location. So e.g. for Pvalb_470522102_m.swc (cells 240 to 269) their x,y,z location does not correspond to the soma location, which has an offset of 237.4944, 233.8336, 35.28.

However, removing the soma recentering in netpyne did not improve the activity discrepancy.

salvadord commented 5 years ago

The only remaining difference I can find are small differences in the 3d point values in the order of <0.01 um. These could be rounding errors arising from the way that netpyne vs bmtk sets the 3d points from the swc files. Despite being small differences, given that they are present in most sections, over time they could accumulate and potentially cause the activity discrepancies.

I compared the swc, netpyne and bmtk values for 3d points, and generally it seems the netpyne values are closer to those in the original swc file, e.g. for cell gid=0 (Scnn1a_473845048_m.swc) below is the comparison of the 8 first 3d points (excluding soma) obtained via .psection() (some differences are highlighted in bold):

SWC: id,type,x,y,z,r,pid 2 3 232.5043 230.7745 35.28 0.2669 1 3 3 231.4495 231.199 35.28 0.2796 2 4 3 230.3513 231.5181 35.28 0.2796 3 5 3 229.2084 231.5159 35.28 0.2542 4 6 3 228.0793 231.3282 35.28 0.2542 5 7 3 226.9616 231.0846 35.28 0.2669 6 8 3 225.8176 231.0617 35.28 0.3305 7 9 3 224.6816 230.9313 35.2803 0.4576 8

netpyne: 2 (232.50430297851562, 230.77450561523438, 35.279998779296875, 0.5338000059127808) 3 (231.44949340820312, 231.19900512695312, 35.279998779296875, 0.5591999888420105) 4 (230.35130310058594, 231.51809692382812, 35.279998779296875, 0.5591999888420105) 5 (229.20840454101562, 231.51589965820312, 35.279998779296875, 0.508400022983551) 6 (228.0792999267578, 231.3282012939453, 35.279998779296875, 0.508400022983551) 7 (226.9615936279297, 231.0845947265625, 35.279998779296875, 0.5338000059127808) 8 (225.81759643554688, 231.06170654296875, 35.279998779296875, 0.6610000133514404) 9 (224.6815948486328, 230.93130493164062, 35.28030014038086, 0.9151999950408936)

bmtk: (3, (232.50450134277344, 230.77499389648438, 35.279998779296875, 0.5338000059127808)) (4, (231.4495086669922, 231.19900512695312, 35.279998779296875, 0.5591999888420105)) (5, (230.35150146484375, 231.51800537109375, 35.279998779296875, 0.5591999888420105)) (6, (229.20849609375, 231.51600646972656, 35.279998779296875, 0.508400022983551)) (7, (228.07949829101562, 231.3280029296875, 35.279998779296875, 0.508400022983551)) (8, (226.96250915527344, 231.0850067138672, 35.279998779296875, 0.5338000059127808)) (9, (225.81849670410156, 231.06199645996094, 35.279998779296875, 0.6610000133514404)) (10, (224.68251037597656, 230.93099975585938, 35.28030014038086, 0.9151999950408936))

salvadord commented 5 years ago

Belos is a snippet of the code I use to extract the 3d points from the swc file:

                # morphology
                morphology_file = self.subs(self.network_config['components']['morphologies_dir']) +'/'+info['morphology'] + '.swc'
                cellMorph = EmptyCell()
                swcData = h.Import3d_SWC_read()
                swcData.input(morphology_file)
                swcSecs = h.Import3d_GUI(swcData, 0)
                swcSecs.instantiate(cellMorph)

                # replace axon with AIS stub
                if self.replaceAxon:
                    fix_axon_peri(cellMorph)

                # extract netpyne parameters
                secs, secLists, synMechs, globs = neuronPyHoc.getCellParams(cellMorph)

I didn't yet find the corresponding bmtk code to do this. @kaeldai if you can point me to this code, I can try to use the same in netpyne, to test whether this morph differences are actually causing the discrepancy in spk times.

kaeldai commented 5 years ago

Turns out import3d.hoc uses different precision depending on if you instantiate with a python class (NetPyNE) or a hoc template (bmtk). I have more info about it here. If you call instantiate() with a hoc template/reference it cut off at 3 to 4 sig figs, whereas a full 64-bit float has a little over 15.

I've tested with both full and limited precision - it doesn't seem to affect the single cell examples but it can lead to differences in the 300 cell network.

kaeldai commented 5 years ago

The code for loading the swc using Hoc Templates in BMTK is here. You'll need to download the following template.

In general the code should look like:

h.load_file(str('path/to/Biophys1.hoc'))
hobj = h.Biophys1(morphology_file.swc)

# Replace the axon
for sec in hobj.axon:
    h.delete_section(sec=sec)
h.execute('create axon[2]', hobj)

for sec in hobj.axon:
    sec.L = 30
    sec.diam = 1
    hobj.axonal.append(sec=sec)
    hobj.all.append(sec=sec)
hobj.axon[0].connect(hobj.soma[0], 0.5, 0)
hobj.axon[1].connect(hobj.axon[0], 1, 0)
h.define_shape()

salvadord commented 4 years ago

discrepancy might be due to netpyne not creating individual synapse objects for connections (netcons) that target the same section and location -- should retest this model with the flag oneSynPerNetcon=True to check if this was source of discrepancy

AllenInstitute / sonata

NetPyNE simulation results for 300 biophys cell network #54