NeuroML / pyNeuroML

A single package in Python unifying scripts and modules for reading, writing, simulating and analysing NeuroML2/LEMS models.
https://docs.neuroml.org/Userdocs/Software/pyNeuroML.html
GNU Lesser General Public License v3.0
36 stars 30 forks source link

[Bug] Output export to Neuron are specific on ions, not channels #196

Closed christian-oreilly closed 1 year ago

christian-oreilly commented 1 year ago

When creating two different types of channels using the same variety of ions ion, e.g.,

    add_channel_density(cell, nml_cell_doc,
                        cd_id="Kdrfast_soma",
                        cond_density="73.37 mS_per_cm2",
                        ion_channel="Kdrfast",
                        ion_chan_def_file=path + "Kdrfast.channel.nml",
                        erev="-77mV",
                        ion="k",
                        group="soma_group")

    add_channel_density(cell, nml_cell_doc,
                        cd_id="KvAolm_soma",
                        cond_density="4.95 mS_per_cm2",
                        ion_channel="KvAolm",
                        ion_chan_def_file=path + "KvAolm.channel.nml",
                        erev="-77mV",
                        ion="k",
                        group="soma_group")

and then trying to plot the current density for these channels, e.g.,

    simulation.add_column_to_output_file("output0", column_id="Kdrfast_soma", 
                                         quantity="pop/0/pyr/0/biophys/membraneProperties/Kdrfast_soma/iDensity/")
    simulation.add_column_to_output_file("output0", column_id="KvAolm_soma", 
                                         quantity="pop/0/pyr/0/biophys/membraneProperties/KvAolm_soma/iDensity/")

the entity recorded is specified as the channel type (as specified in the quantity argument), not the ion variety. The NeuroML file is created accordingly

        <OutputFile id="output0" fileName="amygdala_sim.dat">
            <OutputColumn id="pop_0_pyr_0_v" quantity="pop/0/pyr/0/v"/> 
            <OutputColumn id="Kdrfast_soma" quantity="pop/0/pyr/0/biophys/membraneProperties/Kdrfast_soma/iDensity/"/> 
            <OutputColumn id="Nav_soma" quantity="pop/0/pyr/0/biophys/membraneProperties/Nav_soma/iDensity/"/> 
            <OutputColumn id="HCNolm_soma" quantity="pop/0/pyr/0/biophys/membraneProperties/HCNolm_soma/iDensity/"/> 
            <OutputColumn id="KvAolm_soma" quantity="pop/0/pyr/0/biophys/membraneProperties/KvAolm_soma/iDensity/"/> 
        </OutputFile>

but when simulated with NEURON, it ends up recording the currents that are ion-specific (i.e, potassium), not channel-specific (i.e., Kdrfast_soma, KvAolm_soma):

        # ######################   File to save: amygdala_sim.dat (output0)
        # Column: pop/0/pyr/0/biophys/membraneProperties/Kdrfast_soma/iDensity/
        h(' objectvar v_Kdrfast_soma_output0 ')
        h(' { v_Kdrfast_soma_output0 = new Vector() } ')
        h(' { v_Kdrfast_soma_output0.record(&a_pop[0].soma_0.ik(0.5)) } ')
        h.v_Kdrfast_soma_output0.resize((h.tstop * h.steps_per_ms) + 1)
[...]
        # Column: pop/0/pyr/0/biophys/membraneProperties/KvAolm_soma/iDensity/
        h(' objectvar v_KvAolm_soma_output0 ')
        h(' { v_KvAolm_soma_output0 = new Vector() } ')
        h(' { v_KvAolm_soma_output0.record(&a_pop[0].soma_0.ik(0.5)) } ')
        h.v_KvAolm_soma_output0.resize((h.tstop * h.steps_per_ms) + 1)

and, sure enough, the two recorded signals for Kdrfast_soma and KvAolm_soma end up being identical! I assume this is a bug and not intentional behavior.

sanjayankur31 commented 1 year ago

@pgleeson : this one looks like one for you :)

sanjayankur31 commented 1 year ago

(Note: probably belongs to one of the java packages, probably org.neuroml.export---this does the writing to NEURON, but we can move it there later)

christian-oreilly commented 1 year ago

Yes, since there are different pieces of software interacting (pyNeuroML, libNeuroML, java backend, ...) I was not completely sure where to submit this. The issue seemed to me to be due to LEMSSimulation which is part of pyNeuroML but I did not look further to see if the issue was in an underlying lib.

pgleeson commented 1 year ago

Good point. Will look into this...

pgleeson commented 1 year ago

@christian-oreilly As you guessed, the issue is down to ik, ina being used to record the current, which works fine if there is just one K or Na based channel (e.g. classic HH), but fails if there are 2 or more (just returns the sum). This will require an update of the mod file generated bu jnml from NML channels (which needs to be updated for other changes) and will be worked on over the next 1-2 weeks. Hope this is ok.

christian-oreilly commented 1 year ago

Sure, @pgleeson, that would actually be more than OK but great!

pgleeson commented 1 year ago

@christian-oreilly I have a local version of this working, see below where different K currents are plotted independently. I need to do a few more tests to make sure it doesn't break anything else, and will make a new release of pynml in the next few days.

Screenshot from 2022-11-03 12-22-21

pgleeson commented 1 year ago

Sorry for the delay in getting back about this @christian-oreilly, but the latest pyNeuroML should support this now: https://pypi.org/project/pyNeuroML/#history

See this example for multiple K currents which can be recorded/plotted independently: https://github.com/NeuroML/NeuroML2/blob/development/LEMSexamples/LEMS_NML2_Ex15_CaDynamics.xml#L365

Note though, the units for the current have been corrected to use consistent units in NEURON: https://github.com/NeuroML/org.neuroml.export/pull/89/files#diff-64ca13373b09c711e7b7a15647df2c332806b9bb9f02d230e8bc51543bdecc78, and when they're saved to file they're still saved in SI units.

Note also, that the direction of the current being saved is reversed, i.e. currents into the cell like Na are now positive, to match the default behaviour in jLEMS (as the current handling is defined in the LEMS language).

Let me know if you have questions or issues.

christian-oreilly commented 1 year ago

@pgleeson That looks great! Thanks for having a look at that. It will be very useful to debug models with multiple currents.