scidash / neuronunit

A package for data-driven validation of neuron and ion channel models using SciUnit
http://neuronunit.scidash.org
38 stars 24 forks source link

NeuroMLDBModel current amplitudes that are multiple of rheobase #222

Open rgerkin opened 5 years ago

rgerkin commented 5 years ago

@JustasB The code for loading waveforms by current amplitude has the following problem:

import quantities as pq
from neuronunit.neuromldb import NeuroMLDBModel
model = NeuroMLDBModel(model_id='NMLCL000086')
model.fetch_waveform_list();
model.get_waveform_by_current(100*pq.pA)

leads to:

    136     def get_waveform_current_amplitude(self, waveform):
--> 137         return float(waveform["Waveform_Label"].replace(" nA", "")) * pq.nA
    138 
    139 

ValueError: could not convert string to float: '0.75xRB'

Currents that are expressed as a multiple of rheobase will need their own handling. Is the rheobase stored somewhere in the class or elsewhere in the database? If so, I can add a conversion from "NUMxRB" to the true value.

russelljjarvis commented 5 years ago

@rgerkin

I think you are using the wrong type of waveform/protocol combination.

If you use the logic below you can query the NML-DB server to get appropriate waveforms that are amenable to

def get_waveform_current_amplitude(waveform):
    long_squares = [ w for w in wlist if w['Protocol_ID'] == 'LONG_SQUARE' ]

    #variable_names = [  w["Variable_Name"] for w in wlist if w["Protocol_ID"] == "LONG_SQUARE" ]
    applied_current_injections = [ w for w in wlist if w["Protocol_ID"] == "LONG_SQUARE" and w["Variable_Name"] == "Current" ]

    #import pdb; pdb.set_trace()
    currents = [ w for w in wlist if w["Protocol_ID"] == "LONG_SQUARE" and w["Variable_Name"] == "Voltage" ]

    in_current_filter = [ w for w in wlist if w["Protocol_ID"] == "SQUARE" and w["Variable_Name"] == "Voltage" ]
    rheobases = []
    for wl in long_squares:
        wid = wl['ID']
        url = str("https://neuroml-db.org/api/waveform?id=")+str(wid)
        waves = requests.get(url)
        temp = json.loads(waves.text)
        if temp['Spikes'] >= 1:
            rheobases.append(temp)
    if len(rheobases) == 0:
        return None

    in_current = []
    for check in in_current_filter:
        amp = get_waveform_current_amplitude(check)
        if amp < 0 * pq.nA:
            in_current.append(amp)
    rheobase_current = get_waveform_current_amplitude(rheobases[0])
    druckmann2013_standard_current = get_waveform_current_amplitude(currents[-2])
druckmann2013_strong_current = get_waveform_current_amplitude(currents[-1])

https://github.com/russelljjarvis/neuronunit/blob/barcelona/get_allen_features_from_nml_db.py#L153-L169

rgerkin commented 5 years ago

@russelljjarvis Justas had pointed me specifically to get_waveform_by_current so I was trying to produce a minimal example using that method.

russelljjarvis commented 5 years ago

The code I am proposing uses get_waveform_by_current under the hood too.

rgerkin commented 5 years ago

@russelljjarvis It does? I see that get_waveform_by_current calls get_waveform_current_amplitude but not the other way around. Is there somewhere else in your code where something that calls get_waveform_by_current is used?

russelljjarvis commented 5 years ago

Sorry no, you are right, I used to call those methods indirectly with the three lines below but then I realized that they would re-download waveforms every time a waveform was needed which would slow down the batch process I was designing. Therefore I get the waveforms myself once using code that is very similar, but then I hard assign those waveforms to the model.

The functionality is there by I bypassed it, because I needed it to work with a different pre-existing code for querying the NML-DB waveforms.

get_waveform_by_current is more like a private method that only works under particular circumstances. Its better if the object calls this method itself, by calling these three methods from the object:

self.standard = self.model.nmldb_model.get_druckmann2013_standard_current()
self.strong = self.model.nmldb_model.get_druckmann2013_strong_current()
self.ir_currents = self.model.nmldb_model.get_druckmann2013_input_resistance_currents()

These lines https://github.com/russelljjarvis/AllenDruckmannData/blob/master/dm_test_interoperable.py#L126-#L128

The code below is acting on a neuroml-db static backend model. Since I am supplying a pre-made model, and models type is not None.

        if type(model) is type(None):
            self.model = model_dict[model_id]
            self.model_id = model_id
            if self.model_id not in self.predicted:
                self.predicted[self.model_id] = [None for i in range(38)] # There are 38 tests
            self.standard = self.model.nmldb_model.get_druckmann2013_standard_current()
            self.strong = self.model.nmldb_model.get_druckmann2013_strong_current()
            self.ir_currents = self.model.nmldb_model.get_druckmann2013_input_resistance_currents()

        else:
            self.model = model
            self.standard = model.druckmann2013_standard_current
            self.strong = model.druckmann2013_strong_current
            self.ir_currents = model.druckmann2013_input_resistance_currents

https://github.com/russelljjarvis/AllenEFELDruckmanData/blob/master/neuromldb.py#L109

What I should have said instead is if these lines of code are executed instead:

            self.standard = self.model.nmldb_model.get_druckmann2013_standard_current()
            self.strong = self.model.nmldb_model.get_druckmann2013_strong_current()
            self.ir_currents = self.model.nmldb_model.get_druckmann2013_input_resistance_currents()

They properly set up the neuroml-db static model class in a way, that makes it able to call get_waveform_by_current itself without syntactic errors.

It's a bit complicated to understand, but when I call get_waveform_by_current. I am still adding currents to the models lookup table, because the last thing that happens is that method is this call:

 return self.fetch_waveform_as_AnalogSignal(w["ID"])
russelljjarvis commented 5 years ago

https://github.com/russelljjarvis/AllenEFELDruckmanData/blob/master/get_allen_features_from_nml_db.py#L441-L446

    def not_necessary_for_program_completion(DMTNMLO):
        current = DMTNMLO.model.nmldb_model.get_druckmann2013_standard_current()
        DMTNMLO.model.nmldb_model.get_waveform_by_current(current)
        temp0 = np.mean(DMTNMLO.model.nmldb_model.get_waveform_by_current(DMTNMLO.model.nmldb_model.get_druckmann2013_strong_current()))
        temp1 = np.mean(DMTNMLO.model.nmldb_model.get_waveform_by_current(DMTNMLO.model.nmldb_model.get_druckmann2013_standard_current()))
        assert temp0 != temp1
        return
    _ = not_necessary_for_program_completion(DMTNMLO)    

This code works as an example.

JustasB commented 5 years ago

It looks like I really need to add more documentation on the API. I'll work on that when I find spare moments.

For the question:

Is the rheobase stored somewhere in the class or elsewhere in the database?

In this line: https://github.com/scidash/neuronunit/blob/1c0e4b2221cc32938b39dbf4088cdc516a54d321/neuronunit/neuromldb.py#L38

the data variable will contain the contents of the model API call results e.g. https://neuroml-db.org/api/model?id=NMLCL001129

For cell models, that data variable will have a Rheobase_High key, which stores the value of the rheobase current.