Confirm that NeuronUnits jNeuromLB backend really can be regarded as a ground truth.

[russelljjarvis]

@rgerkin I am not sure if the significance of this result was understood. From my perspective it meant that jNeuroML models could not be regarded as a ground truth until this issue was resolved. Perhaps the issue was resolved, since I could not run the final notebooks it was not possible to evaluate the status of this issue. @ChihweiLHBird, since your working installation of NU is likely to be a lot closer to Rick's I was wondering if you would be able to externally verify Rick's fixes at the very bottom of this thread?

[russelljjarvis]

backend/model/test interaction has built in assumptions now about model.set_stop_time and `test.protocol['sim_length']``

Consequences of this bug:

• A waveform that should last a whole 1second lasts for 500 ms.
• The NU tests cannot tell that waves are contracted in time.
• Scores for Capacitance Test Time Constant Test and Spike Width test are artificially improved due to the time contraction distortion.

The notebook contents need translating to be compatible with the newer versions of NU. Code below reproduces minimal version of chapter3 containing the bug

![NeuronUnit Logo](https://raw.githubusercontent.com/scidash/assets/master/logos/neuronunit-logo-text.png)
# Chapter 3
Back to [Chapter 2](chapter2.ipynb)

%matplotlib inline
import os,sys
import numpy as np
import matplotlib.pyplot as plt
import quantities as pq
import sciunit
import neuronunit
from neuronunit import aibs
from neuronunit.models.reduced import ReducedModel
import quantities as pq
from neuronunit import tests as nu_tests, neuroelectro
neuron = {'nlex_id': 'nifext_50'} # Layer V pyramidal cell
tests = []

dataset_id = 354190013  # Internal ID that AIBS uses for a particular Scnn1a-Tg2-Cre 
                        # Primary visual area, layer 5 neuron.
    
try:
    assert suite is not None
    assert len(suite) == 8
    

except:
    import quantities as pq
    from neuronunit import tests as nu_tests, neuroelectro
    neuron = {'nlex_id': 'nifext_50'} # Layer V pyramidal cell
    tests = []

    dataset_id = 354190013  # Internal ID that AIBS uses for a particular Scnn1a-Tg2-Cre 
                            # Primary visual area, layer 5 neuron.

    # Obtain the empirical rheobase current from the Allen Brain Insitute Cell Types database.  
    observation = aibs.get_observation(dataset_id,'rheobase')
    rheobase_test = nu_tests.RheobaseTest(observation=observation)
    tests += [rheobase_test]

    test_class_params = [(nu_tests.InputResistanceTest,None),
                         (nu_tests.TimeConstantTest,None),
                         (nu_tests.CapacitanceTest,None),
                         (nu_tests.RestingPotentialTest,None),
                         (nu_tests.InjectedCurrentAPWidthTest,None),
                         (nu_tests.InjectedCurrentAPAmplitudeTest,None),
                         (nu_tests.InjectedCurrentAPThresholdTest,None)
                        ]

    # Obtain all other parameters from neuroelectro.org.
    for cls,params in test_class_params:
        observation = cls.neuroelectro_summary_observation(neuron)
        tests += [cls(observation,params=params)]

    # A hook to update all tests after the RheobaseTest to use the rheobase current (to produce exactly one AP)
    def update_amplitude(test,tests,score):
        rheobase = score.prediction['value']
        for test in tests[1:]:
            if 'Injected' in test.name:
                # Set current injection to just suprathreshold
                test.params['injected_square_current']['amplitude'] = rheobase*1.01 

    hooks = {tests[0]:{'f':update_amplitude}}
    suite = sciunit.TestSuite(tests,name="vm_suite",hooks=hooks)
    suite[0].observation['mean'] = suite[0].observation['value']  
    suite[0].observation['std'] = suite[0].observation['value']    

    with open('suite.p','wb') as f:
        pickle.dump(suite,f)
# This example is from https://github.com/OpenSourceBrain/IzhikevichModel.
for rel_path in ['.','models','../models']:
    # Check to see if this is the neuronunit/docs directory.
    DOCS_PATH = os.path.abspath(os.path.join(os.getcwd(),rel_path)) 
    if DOCS_PATH.endswith('neuronunit/models'):
        break
assert DOCS_PATH.endswith('neuronunit/models'), ("Could not find the path to neuronunit/docs. "
                                               "Change the DOCS_PATH variable to the path to "
                                               "your neuronunit/docs directory.")
LEMS_MODEL_PATH = os.path.join(DOCS_PATH,'NeuroML2/LEMS_2007One.xml')
model = ReducedModel(LEMS_MODEL_PATH,name='vanilla',backend='jNeuroML')

suite[0].params
for t in suite:
    t.params2 = {}
    t.params2['injected_sqaure_current'] = t.params
    t.params = t.params2

suite[0].params
model = None
model = ReducedModel(LEMS_MODEL_PATH,name='vanilla',backend='jNeuroML')
try:
    model.inject_square_current(suite[0].params['injected_sqaure_current'])
except:
    model.inject_square_current(suite[0].params)

vm = model.get_membrane_potential()
len(vm)
vm.sampling_rate
vm.times[-1]/len(vm)
print(len(vm))
import matplotlib.pyplot as plt
plt.plot(vm.times,vm.magnitude)
plt.show()

suite[0].params = suite[0].params['injected_sqaure_current']
print(suite[0].params['delay'])
print(suite[0].params['duration'])

This last part reveals the bug:

sim_length = suite[0].params['delay'] + suite[0].params['duration']
print(vm.times[-1] == sim_length)
print(vm.times[-1],sim_length)
print(vm.times[-1] == sim_length/2.0)
print(vm.times[-1],sim_length/2.0)

For contrast you can compare with a newer version of Izhikevich model which has accurate Neo sampling period, and scores less well on those NU tests.

from neuronunit.optimisation.data_transport_container import DataTC
from neuronunit.optimisation.optimization_management import inject_and_plot_model, inject_and_not_plot_model

dtc = DataTC()
dtc.backend = "RAW"
fast_model = dtc.dtc_to_model()
dtc.attrs = fast_model.default_attrs
results = inject_and_not_plot_model(dtc)

print(len(results))
print(len(vm))
inject_and_plot_model(dtc)

[russelljjarvis]

The real scores for Spike width, Capacitance, and Time Constant would be much worse than these.

[rgerkin]

@russelljjarvis Are you sure it is always exactly one half of the correct length for any duration of injected current? When I run the model with jNeuroMLBackend, I get a truncated trace (which for some durations of current could be half of the desired value), but this results from the total simulation duration being too short. I can then update that total simulation duration to fix the problem. Can you check if that also works for you; if so, then the fix is to make sure that the simulation duration is always sufficiently long. This is supposed to happen automatically but apparently isn't working. Example below:

[russelljjarvis]

Wrong dt was just my best guess about why the model recording would be truncated. It was just speculation.

Okay I will check if setting model.set_stop_time(1500*pq.ms) helps later. Would setting the model stop time also guarantee that the tests where experiencing the right waveform?

[russelljjarvis]

Too much stuff has changed in sciunit to test this out. I would have to solve other compatibility bugs. Can you plot the vms that were used by the tests? If the vms used by tests are truncated then I think that is a more serious problem.

[rgerkin]

Not much has changed in sciunit, functionally (it has all been adding type hints, more documentation, and more unit tests), but if you experience any problems you can use an older commit or switch to the metacell branch which has been stable for a while.

In any case I sent you a notebook on Slack showing the behavior of the test and the plotted vm in this situation.

[russelljjarvis]

Okay so I am not sure does this mean the issue is fixed? It looks like I was not able to look at the notebooks you are talking about. Because our branches are too far apart.

This error in the NU test protocol truncating waveforms is propagating surprisingly far. It is increasingly looking like it affects optimization work.

Okay looking at notebook now. I had to reinstall sciunit very recently, so I am noticing compatibility differences I may have been previously immune to for unforeseeable reasons.