NetPyNE simulation for one_cell soma_cm2_hh model (sim_tests)

[salvadord]

I'm trying this example: https://github.com/AllenInstitute/sonata/tree/sim_tests_dontmerge/examples/sim_tests/biophysical/one_cell/soma_cm2_hh

But I get an error reading the neuroml biophys params, specifically this line: https://github.com/AllenInstitute/sonata/blob/sim_tests_dontmerge/examples/sim_tests/shared_components/templates/nml/soma_cm2_hh.nml#L9 (because 'condDensity' is '80' without units). This is the error:

/usr/local/lib/python3.7/site-packages/pyneuroml/pynml.py in convert_to_units(nml2_quantity, unit, verbose)
    300         if un.symbol == unit:
    301 
--> 302             new_value = si_value / (un.scale * pow(10, un.power)) - un.offset
    303             if not un.dimension == dim:
    304                 raise Exception("Cannot convert %s to %s. Dimensions of units (%s/%s) do not match!"% \

TypeError: unsupported operand type(s) for /: 'NoneType' and 'float'

The code I'm using is here (which is adapted from the neuroml import code used in netpyne): https://github.com/Neurosim-lab/netpyne/blob/sonata/netpyne/conversion/sonataImport.py#L651

I have pyneuroml version '0.3.13'

@pgleeson - any suggestion? maybe I need to update some of the neuroml conversion code?

[pgleeson]

condDensity certainly would need units, but as discussed with @wvangeit in another issue (I forget where, Werner?) that line is not required, as the -80 (mV) for el_hh is set in the previous line...

[salvadord]

thx, I'll try removing the line.

so @kaeldai, this doesn't seem to be an issue for bionet? can you point me to the bionet code that reads the biophys params from the neuroml file? thx

[salvadord]

I removed the problematic line as suggested by padraig: https://github.com/AllenInstitute/sonata/blob/sim_tests_dontmerge/examples/sim_tests/shared_components/templates/nml/soma_cm2_hh.nml#L9

But now getting a different error:
AttributeError: 'nrn.Mechanism' object has no attribute 'gbar'

@pgleeson - this seems to come from incorrectly reading the hh 'gnabar' and 'gkbar' params, even though I'm using the same code that is used by netpyne to import neuroml. Any suggestions?

In [5]: sim.net.params.cellParams['biophysical_soma_cm_100']['secs']
Out[5]: 
{'soma_0': {'geom': {'L': 1.0,
   'nseg': 1,
   'diam': 1.0,
   'Ra': 35.4,
   'cm': 1.0,
   'pt3d': [(-0.5, 0.0, 0.0, 1.0),
    (0.0, 0.0, 0.0, 1.0),
    (0.5, 0.0, 0.0, 1.0)]},
  'topol': {},
  'mechs': {'hh': {'gbar': 0.005, 'e': -80.0}},
  'ions': {'na': {'e': 50.0}, 'k': {'e': -85.0}}}}

[pgleeson]

The root of this seems to be specifying only part of the model in NeuroML. Ideally for this hh should be a NeuroML compliant ion channel representing one active conductance with one erev and one cond density. Here it's clearly referring to the NEURON hh mod file/mechanism with 3 independent conds/currents for na, k and pas... When you mix the two all bets are off...

[wvangeit]

Mmm, interesting. So if I understand this correctly, NeuroML will only be able to interact with MOD files that are complying to some standard. Is this standard defined somewhere, so that we could link from it from the SONATA format description ?
(And in this case we should also investigate if the incompatibility with hh might be a problem for future models in SONATA).

[pgleeson]

To be clear, NeuroML isn't designed to refer/interact with mod at all. Yes, there can be mod files generated from NML channel descriptions, but that's after everything including the cond densities/cells/channels are specified in valid NML.

Here the SONATA standard/BMTK is using an ad hoc combination of hoc & NML which is outside the scope of standard NeuroML, and it for that toolchain to decide what's allowed.

To summarise though, a mod file with just a single type of conductance (representing one voltage gated current) and a public parameter for conductance density and reversal potential (which is most of the mod files on ModelDB) should be fine, though the names for the parameters for cond dens/rev pot should probably be changed (maybe gbar and e?) so they can be found by the code linking the NML to the mod file.

[salvadord]

So @pgleeson what do you propose? Should I just add an ad-hoc case to read the 'hh' mechanism parameters in the current .nml file, or should the .nml file be modified/extended ?

@kaeldai could you please point me to the bionet code that is able to properly read+interpret this nml file?

[pgleeson]

@salvadord I don't see any reason why the hh mechanism should be used/encouraged in the standard examples, it would just lead to confusion as there are multiple independent conductances in there. It would be best to use separated Na/K mechanisms and possibly the inbuilt pas. Perhaps the mod files from the standard Allen/Hay model could be reused @wvangeit?

[salvadord]

A couple things:

• I noticed for the soma_cm2_hh the 'hh' na and k conductances were set to 0. So this would be equivalent to having just the 'pas' mech. I therefore upated the 'soma_cm2_hh.nml' file to have the following single channelDensity entry:
  <channelDensity id="g_pas_soma" segmentGroup="soma" ion="non_specific" ionChannel="pas" erev="-80 mV" condDensity="0.005 S_per_cm2"/>

• This allowed me to import the sonata model and run it in netpyne. However, the result I get is different from the one generated by bmtk in 'expected_output/membrane_potential.h5'. To check I also implemented the same model directly in Python/NEURON. The graph shows the comparison of bmtk vs netpyne vs NEURON/Python:
  

• The NEURON/Python code is shown below:

# NEURON version
from neuron import h, gui
from matplotlib import pyplot as plt

h.v_init = -80

# create soma section
soma = h.Section(name='soma')

# set geometry
soma.L = soma.diam = 1 # Makes a soma of 500 microns squared.
soma.cm = 1.0

# Insert active Hodgkin-Huxley current in the soma
soma.insert('hh')
soma.gnabar_hh = 0.0  # Sodium conductance in S/cm2
soma.gkbar_hh = 0.0  # Potassium conductance in S/cm2
soma.ena = 50
soma.ek = .85
soma.gl_hh = 0.005    # Leak conductance in S/cm2
soma.el_hh = -80     # Reversal potential in mV

# add current clamp stimulus
stim = h.IClamp(soma(0.5))
stim.amp = 0.002  # input current in nA
stim.delay = 20  # turn on after this time in ms
stim.dur = 30  # duration of 1 ms

# record soma voltage and time
t_vec = h.Vector()
v_vec_soma = h.Vector()
v_vec_soma.record(soma(0.5)._ref_v) # change recoding pos
t_vec.record(h._ref_t)

# run simulation
h.tstop = 60 # ms
h.run()  

plot voltage vs time
plt.figure(figsize=(8,4)) # Default figsize is (8,6)
plt.plot(t_vec, v_vec_soma, 'b', label='soma')
plt.xlabel('time (ms)')
plt.ylabel('mV')
plt.legend()
plt.show()