Test circuit for running simulations

[arsenius7]

One of the tasks that we have at hand is running a simulation of a test SONATA circuit.

Unfortunately, due to some internal format simulations, we can not fully support all of SONATA features at the moment. Some of the issues were mentioned in the emails previously; I believe having a ticket would facilitate the discussion. Should you believe that emails are a better communication channel in this case, please let me know.

@kaeldai, can you please provide a test SONATA circuit simulation which: 1) Uses only biophysical neurons, not a mix of biophysical + point_process 2) Defines the electrical models in NeuroML rather than JSON

Also, can you please clarify how virtual neurons are participating in the simulation? Why do they have only two coordinates in positions and not three? Are these coordinates of any importance for the sake of the simulation?

I guess more questions will pop up later.

[arsenius7]

'Sonata_example_04_16_18' content is a bit self-contradictory for 'v1' node population. Based on HDF5 contents, GIDs 100-448 are point neurons (they belong to node group '1' marked as 'point'); but 'model_type' in CSV file marks them as 'biophysical'. Should we ignore 'v1/node_group_id' metadata? What it is there for then?

[arsenius7]

Also, I'm not sure how changing v1/1 node group to 'biophysical' affected the edges part. The incoming edges for GIDs 100-448 should have become 'biophysical' as well, right? (i.e., they should now have 'sec_id' and 'sec_x'). Are the all synapses effectively 'on soma' for these neurons?

[wvangeit]

@kaeldai Could you also inject a step current (or maybe even better several) in each of the cell models (1 per nml file), and provide the resulting voltage traces (without any synapses added to the model). This would be very useful to validate if our e-model converter works correctly.

[kaeldai]

It shouldn't be a problem but we're pretty busy this week so I probably won't be able to get it done until at least the end of the week. Assuming we use ~7 biophysical models we have in been using in the past - do you want 7 individual networks and simulations of 1 cell, or I can also create a single unconnected network with 7 cells?

[wvangeit]

The latter is ok. As long as there are no synapses activated we can use it to test the e-model. (We don't even necessarily need a network in sonata format. If you give us, for every nml file, the current that was injected and the voltage response we can compare the results.

[arsenius7]

@kaeldai , we are concerned there is a possible mismatch in section numbering (see #16). Could you please provide x, y, z coordinates for all of the synapses in 1-cell circuit ("absolute ones", i.e. after the cell is placed in the circuit and rotated accordingly)? That could serve us as an indirect check that we are on the same page wrt to cell rotation and section ID order. Thanks.

[arsenius7]

As far as I understand, output/cell_vars.h5:/v/data is soma voltage for the only biophysical cell in 1-cell example. When I compare it with input spike times (coming from inputs/virtual_spike_trains.h5), I get the following picture:

In my (very limited) understanding of cell electrical behavior, the voltage trace should be more "noisy", reacting to every input spike to some degree (until the accumulated voltage reaches some threshold). Is there some factor which I am missing which causes the trace to be that smooth?

[wvangeit]

@kaeldai so more specifically about the trace above: Could it be that we don't see EPSPs in your trace for all the presynaptic events because e.g.:

• your synapses have a very high failure rate ?
• a significant amount of synapses have a 0 weight ?
• do we have the wrong list of presynaptic cells ?
• ...

[arsenius7]

Could you please provide x, y, z coordinates for all of the synapses in 1-cell circuit ("absolute ones", i.e. after the cell is placed in the circuit and rotated accordingly)?

Not necessarily xyz-coordinates, in fact. Path distance + section type for each synapse would be also useful.

[CellAssembly]

Thanks all. For the membrane trace, we are looking into it. It is possible that most connections are so distal that they have negligible effect on the soma except for one or two that are right at the soma which is causing the neuron to "instantly" fire in response to those inputs. May have to send another example.

[kaeldai]

@arsenius7 @CellAssembly It turns out my new code was not reading the Sonata input spike-trains file correctly. I fixed it and regenerated the network and now we get something more realistic.

I also updated network files so that in the edges.h5 file has positions of every synapse plus arc-distance from soma center. The new files can be found on the shared google drive as 1_cell_synaptic.tar.gz.

@wvangeit Also on the google drive there I've upload 5_cells_iclamp.tar.gz. It contains sonata network and simulation of 5 different biophysical cells each receiving (the same) current clamp input. I used a series of 3 increasing step currents with a 500 ms delay.

[arsenius7]

I also updated network files so that in the edges.h5 file has positions of every synapse plus arc-distance from soma center.

Thanks, @kaeldai , that was super helpful.

Alas, I could not match the synapse positions by x, y, z coordinates only; and I could not find arc-distance in the edges.h5 file you provided.

Could you please include these as well (+ideally associated section type for each synapse position, i.e. whether it's basal dendrite, apical dendrite etc).

I have tried to figure it out from build_network.py script, but btmk branch I see on GitHub does not have SWCReader class used there. If you can provide btmk branch used for running that script, that might also help.

Thanks.

[kaeldai]

@arsenius7

Sorry, I got a little careless trying to send out the example before going home on Friday. I updated the example on the drive so the edges.h5 now includes a distance from the soma and swc type (3=basal dendrite, 4=apical dendrite). I also fixed the bmtk code in case you want to play around with the build script.

Note: we use the built-in h.distance() method in NEURON to determine distance. NEURON tries to interpolate the center of a straight segment which can be a little different than calculating the precise synaptic coordinates from the swc morphology. Functionally it doesn't matter but if you plot it you may see some floating synapses.

Also because segments can be up to 40 micrometers in length some of the distances are outside the [50, 150] um range. As long as part of the segment falls within the range there's a chance that we'll place a synapse on it even if the center falls outside.

[arsenius7]

It turns out my new code was not reading the Sonata input spike-trains file correctly. I fixed it and regenerated the network and now we get something more realistic.

@kaeldai , did that affect 450-cells circuit you provided earlier in April? If that's the case, can you please provide an updated version of that one? Thanks.

[kaeldai]

The 450 cell circuit I sent in April wouldn't have had that issue - the spike trains were stored in nwb files.

[arsenius7]

@kaeldai , where can I find GID mapping for 450-cell circuit? Also, it looks like gids and timestamps in output/spikes.h5 are not of same length. Do I get it right that the last value in gids should be ignored? Thanks.

[arsenius7]

Also, is the network built in a way that it should exhibit some specific spiking behavior (say, particular subgroup constantly spiking while other cells being "silent")? Or is it more like a "random" example?

[kaeldai]

@arsenius7

For these examples there is just a one-to-one correspondence between the gids and the non-virtual nodes. I did however generate a gids h5 file anyways.

That extra row in gids is an old bug. I've reran the simulation and confirmed an extra row was being added to the gids dataset.

As far as the example network - it was built by taking a 1% random sampling of an existing 45,000 cell cortical layer model. The original network has specific spiking behavior, but I doubt you'd see anything meaningful in the sub-sample. I only use it for testing purposes and just assume the spiking is random. We'd need to talk to @CellAssembly or possibly someone from Costa's group if you want a more biologically relevant example

[arsenius7]

It's more about understanding the expectations, actually.

For instance, there seems to be a certain pattern in the raster plot of output/spikes.h5:

As far as I understand, it's mostly LIF_inh cells which are spiking (node_type_id=100000102). Is it something we aim to reproduce?

[hernando]

Instead of a google drive, can we use the github repository to store circuit and simulation examples?

[arsenius7]

@kaeldai , can you please comment how should we interpret syn_weight for point_process neurons converted to biophysical?

These appear to be on a different scale:

• for biophysical ones (edge group 0 in all edge populations): 0.1 - 1
• for former point_type (edge group 1 in all edge populations): 10 - 100

Thanks.

[kaeldai]

@arsenius7 What file are you looking at? In NEURON the weight depends on the post-synaptic neuron model. For point_process models (IntFire1, IntFire2, etc), the syn_weight should be in the range [-1.0, +1.0]. It's similar to the weight in a standard artifical neural network.

For biophysical models the weight is roughly similar to the synaptic conductance and has a range from (0.0, +\inf). To determine if a synapse is inhibitory or excitation depends on the erev value of the synapse.

[kaeldai]

I'm fine with uploading the examples to github. However I would recommend we do so in a separate repo, especially if we want to share compartmental simulation results. Github can become problematic when repos starting getting into the GBs range.

[arsenius7]

syn_weight should be in the range [-1.0, +1.0]

Sorry, forgot to undo x1000 multiplication (in the current BBP format we store conductance is nanosiemens). So for biophysical the range was (1e-4; 1e-3); and for point_process (1e-2; 1e-1).

Anyway, once these cells are marked as biophysical in the nodes files, as is now the case in 'Sonata_example_04_16_18', we should treat the corresponding edges as biophysical as well, and interpret syn_weight accordingly, no?

Maybe we should remove former point_process nodes (node group 1 in v1 population) to avoid any misinterpretation like that and generate a test circuit with only biophysical + virtual neurons?

Former point_process nodes (node types 100000101, 100000102) would also affect the visualisation, since they are given some morphology w/o configuring its rotation angles.

[hernando]

@kaeldai I was thinking of something in the order of KB if possible, 100s of cells in the circuit including all populations and just a dozen cells and a few timestamps per report.

[kaeldai]

I PR'd an example recurrent network of 300 cells into the examples/ directory. You should use this instead of the 450 cell example, as this one shouldn't have the issue with invalid biophys syn_weight.

Also as requested there is a smaller 9 cell example that contains both inhibitory and excitatory inputs.

[arsenius7]

Thanks, @kaeldai.

Could you please clarify what does model_processing=NONE stand for?

According to the spec:

For model_processing="fullaxon", the biophysical neuron will construct and simulate the full axon. This is the default behaviour if model_processing is undefined for a given node.

Is NONE treated as fullaxon for the circuits / simulations provided?

[kaeldai]

Currently for NML based files our simulator implements the equivalent of "aibs_perisomatic", it cuts the axon and replaces it two 30 um cylinders.

[kaeldai]

I will go ahead and update model_processing before the PR gets merged.

[arsenius7]

OK, thanks.

[arsenius7]

Status update.

We have run simulations at our side by converting 9-cells and 300-cells circuits to BBP internal format.

The resulting voltage traces for 9-cells simulation:

Spike raster plot for 300-cells simulation:

There are still some discrepancies; we will continue to investigate the reasons for those. Please let us know your opinion wrt the plots provided. Thanks.

[CellAssembly]

Thanks @arsenius7. This is exciting to see. One quick thing to check would be the time-step in both simulations. Otherwise, we may need to check mod-files and/or synaptic mechanisms.

[kaeldai]

If it helps I should be able to run the network with synapses turned off and just inject a series of current clamp? It could help us determine if the difference is cause by the synapses or something else.

[wvangeit]

@CellAssembly We're now using Exp2Syn for synapses and exactly the same MOD files as you provided for the ion channels. The time step should be ok (dt = 0.1, correct ?) @kaeldai We also have simulated the 5-iclamp circuit you provided. The e-models are 'almost' acting the same, except for some accommodation in gid 2. Didn't manage to trace the reason for that down yet. I was suspecting the calcium, but there isn't much difference there either. But yes, it could be useful if you also could inject step currents in a disconnected version of the 9-clamp circuit, so that we compare these particular cells too.

[wvangeit]

Btw, could somebody give an answer to this issue: https://github.com/AllenInstitute/sonata/issues/22. We've already seen that changes in the nseg can generate differences of a magnitude as seen in the previous comments. At the moment we're using this as equation to calculate the nseg: nseg = 1 + 2*int(L/dL) Btw, I also saw that in the json files for the e-models in the circuits above, the code to replace the axon sets nseg too. Do you know what the final nseg is for the AIS ? Is it the one set in the AIS-replace code, or the one calculated by using dL (in the simulation config) ?

[kaeldai]

@wvangeit We use a slightly different equation: nseg = 1 + 2int(L/(2dL)). I'm not sure why use use 2*dL but I can try to ask around.

Also another issue that can have a significant effect is the way the axon is cut. Here is the method we use to cut the axon. If you want to send us your method I can integrate it into our simulator and check the results.

[wvangeit]

@kaeldai ok, thx for that equation. Mmm, but that's interesting, it actually means that with a dL = 20 it's the same equation we had been using all along :-) nseg = 1 + 2*int(L/40)

About the axon cutting. That is already being taken care of. Unless there is a bug, I should be using exactly the same.

I guess that if we want to get the results any closer I'm going to have to start instantiating the model using bionet on our systems.

[wvangeit]

After a debugging session I found the problem that caused the biggest discrepancy in the traces above. In the files you provided in https://github.com/AllenInstitute/sonata/issues/15#issuecomment-386763831 there is a difference in the values for e_pas between the nml and json files for cell 473863510 (in nml: -85.6125717163 mV, in json: -85.07815551757812 mV). The latter values seems to be the one you used to create the output files, but since we based our files on the nml format, we used the former.

The tiny differences in the other traces were caused by us using 'nseg = 1 + 2int(L/20)' instead of 'nseg = 1 + 2int(L/40)' (since we were not sure about the definition of dL before)

So I'm now getting a perfect match in all the traces fo the 5 cell circuit:

[CellAssembly]

This is fantastic @wvangeit ! Really great detective work and wonderful to see the perfect match. I think we just need to do a sanity check with a raster plot for the full network and, if that works, then this issue can be closed (if everyone agrees).

[wvangeit]

After another Sherlock Holmes session, we discovered the following. When running the 9-cell network (from sonata repo master) using bmtk, we get the following info for the inhibitory synapses of the first gid (from a print statement I added to bmtk):

secname: Biophys1[0].dend[27], seg_x: 0.5, syn.e: -70.0
secname: Biophys1[0].dend[25], seg_x: 0.5, syn.e: -70.0
secname: Biophys1[0].dend[18], seg_x: 0.5, syn.e: -70.0
secname: Biophys1[0].dend[37], seg_x: 0.5, syn.e: -70.0
secname: Biophys1[0].dend[29], seg_x: 0.5, syn.e: -70.0
secname: Biophys1[0].dend[26], seg_x: 0.5, syn.e: -70.0
secname: Biophys1[0].dend[18], seg_x: 0.5, syn.e: -70.0
secname: Biophys1[0].dend[25], seg_x: 0.5, syn.e: -70.0
secname: Biophys1[0].dend[29], seg_x: 0.5, syn.e: -70.0
secname: Biophys1[0].dend[18], seg_x: 0.5, syn.e: -70.0

As you can see there is not a single apical section on that cell with an inhibitory synapse. However, when we instantiate this cell in our simulators, we get 1 out of the 10 synapses placed on the apical dendrites (the 4th synapse on apic[3] to be precise, the 9 others on basal). It seems the issue might be on your side, because you also provided us with the type of the dendrites on which the synapses are placed in the inhvirt_cortex_edges.h5 file. These types match exactly with the ones in our simulator. We also did some point checks for the excitatory synapses on that cell, and they match with our simulator too.

So it seems for some reason the bmtk is not correctly loading the section ids. One issue that I think might already partly explain part of the problem (but doesnt solve it fully), is the fix_axon code you use. If I understand the code correctly, it removes all the axonal sections, and adds 2 new sections add the 'end' of the neuron's 'all' sectionlist. However, since you use this 'all' list to index the sectionids, removing part of the list in the middle might generate problems. (unless you fix this in some place code I'm not aware of) The code on our side makes sure that after the replace_axon code is executed, the sectionid still point to the correct sections as seen by tools that didnt remove the axonal sections (like our touch detector).

[wvangeit]

I discovered this line in the bmtk: https://github.com/AllenInstitute/bmtk/blob/develop/bmtk/simulator/bionet/biocell.py#L134 That seems to be main issue.

for sec in self.hobj.all:
            for _ in sec:
                secs.append(sec)  # section to which segments belongs

This double loop adds an extra section to the list for every 'compartment'. So it means sections appear several times in the list, and since this list is used to index the sections, we get a wrong match. (I think the fix_axon issue still needs to be fixed too to get a full match)

[wvangeit]

To confirm the theory, I patched bmtk to remove the double loop, and I fixed the fix axon function using:

    for i, sec in enumerate(hobj.axon):                                                                                                                                   
        if i < 2:                                                                                                                                                         
            sec.L=30                                                                                                                                                      
            sec.diam=1                                                                                                                                                    
        else:                                                                                                                                                             
            sec.L=1e-6                                                                                                                                                    
            sec.diam=1                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                
    h.define_shape()         

It sets the L of the other axonal sections to 1e-6, which obviously is a hack (to be clear, I didnt remove a single axonal section), but at least for the first gid it solves all the issues: (the tiny difference in spike height you might see in 1 of the spikes is caused by a recording issue in our simulator, it doesnt affect the simulation itself)

However, there are still issues with the other gids as can be seen from the plot below. I wonder if it could have anything to do with at least 1 of these morphologies not having an axon at all.

[wvangeit]

I can confirm that at least in case of gid=4 at least a part of the problem seems to have been caused by my temporary hack of fix_axon above. I'll let you first fix the sectionid handling in the bmtk. There is a chance after this the match for all cells is perfect.

[arsenius7]

Maybe the traces we obtain on our side for 9-cells circuit could be of some use for the debugging: soma.h5.gz

Sorry, HDF5 layout does not follow SONATA spec; but should be pretty straightforward (each column in /data corresponds to soma voltage trace for the corresponding GID). Please let me know if you'd like to see those traces in some other format.

@wvangeit, can you please confirm that these are the same traces you used for the plots above, to avoid any miscommunication on our side? Thanks.

[wvangeit]

@arsenius7 if this file is from the sim at yourhomedir/workspace/tmp/NSETM-504/n9/20180705/simulations/002/BlueConfig, then yes, I used these for the plot.

[arsenius7]

Yes, it's that one.

[wvangeit]

I also patched the hoc files on our end to reflect my hacked version of the fix_axon I added to bmtk. I now get a perfect match for the 9-cell network. Will now try to rerun the 300-cell network.

fyi, the current hoc code I'm using for our replace_axon function (EDIT small fix below, recalculate nseg in case nSec > 0)

proc replace_axon(){ local nSec, i1, i2                                                                                                                                   
  // preserve the number of original axonal sections                                                                                                                      
  nSec = sec_count(axonal)                                                                                                                                                
  // Try to grab info from original axon                                                                                                                                  
  if(nSec == 0) { //No axon section present                                                                                                                               
    i1 = i2 = 0                                                                                                                                                           
  } else if(nSec == 1) {                                                                                                                                                  
    access axon[0]                                                                                                                                                        
    i1 = i2 = v(0.0001)                                                                                                                                                   
  } else {                                                                                                                                                                
    access axon[0]                                                                                                                                                        
    i1 = i2 = v(0.0001)                                                                                                                                                   
    access soma distance() //to calculate distance from soma                                                                                                              
    count = 0                                                                                                                                                             
    forsec axonal{                                                                                                                                                        
      if (count == 1){                                                                                                                                                    
        i2 = v(0.0001)                                                                                                                                                    
        break                                                                                                                                                             
      }                                                                                                                                                                   
    }                                                                                                                                                                     
  }                
  // get rid of the old axon                                                                                                                                              
  if (nSec == 0) {                                                                                                                                                        
      forsec axonal{                                                                                                                                                      
        delete_section()                                                                                                                                                  
      }                                                                                                                                                                   
      create axon[2]                                                                                                                                                      
      access axon[0] {                                                                                                                                                    
        L = 30                                                                                                                                                            
        diam = 1                                                                                                                                                          
        nseg = 1 + 2*int(L/40)                                                                                                                                            
        all.append()                                                                                                                                                      
        axonal.append()                                                                                                                                                   
        v(0.0001) = i1                                                                                                                                                    
      }                                                                                                                                                                   
      access axon[1] {                                                                                                                                                    
        L = 30                                                                                                                                                            
        diam = 1                                                                                                                                                          
        nseg = 1 + 2*int(L/40)                                                                                                                                            
        all.append()                                                                                                                                                      
        axonal.append()                                                                                                                                                   
        v(0.0001) = i2                                                                                                                                                    
      }                                                                                                                                                                   
      nSecAxonal = 2                                                                                                                                                      
      soma[0] connect axon[0](0), .5                                                                                                                                      
      axon[0] connect axon[1](0), 1                                                                                                                                       
   } else {                                                                                                                                                               
      isec = 0                                                                                                                                                            
      forsec axonal {                                                                                                                                                     
        if (isec < 2) {                                                                                                                                                   
            L = 30                                                                                                                                                        
            diam = 1                                                                                                                                                      
        } else {                                                                                                                                                          
            L = 1e-6                                                                                                                                                      
            diam = 1                                                                                                                                                      
        }
        nseg = 1 + 2*int(L/40)
        v(0.0001) = i1
        isec = isec + 1                                                                                                                                                   
      }                                                                                                                                                                   

   }
}                            

[wvangeit]

I reran the 300 cell network, we're getting very close now:

It 'seems' the difference now originates from a single spike triggered in the BBP version and not in the AIBS version (at around 30 ms in plot below). When I look at the voltage traces of that particular neuron, they start diverging close to threshold, without any incoming spikes present. I wonder if it might be due to roundoff errors.

[wvangeit]

@kaeldai did you already have time to look up the code to make the bmtk dump the first time point to the voltage output files ? (so that I can check if it corresponds exactly with v_init)

[salvadord]

See here results for 300_cells in NetPyNE: https://github.com/AllenInstitute/sonata/issues/54