Check voltage (#106)

* Fix to sAHP

* Warn the user when no spikes are detected while the voltage is above the
detection threshold

* fix AttributeError: can't set attribute

* fix IndexError

* update the efel setting in case of autothreshold

* add test_auto_threshold_detection

Co-authored-by: Anil Tuncel <anil.tuncel@epfl.ch>

bluepyefe/ecode/DeHyperPol.py

@@ -103,7 +103,7 @@ def interpret(self, t, current, config_data, reader_data):
 
         # Converting back to ms
         for name_timing in ["ton", "tmid", "toff"]:
-            self.timing_index_to_ms(name_timing, t)
+            self.index_to_ms(name_timing, t)
         self.tend = len(t) * self.dt
 
     def generate(self):

bluepyefe/ecode/HyperDePol.py

@@ -104,7 +104,7 @@ def interpret(self, t, current, config_data, reader_data):
 
         # Converting back to ms
         for name_timing in ["ton", "tmid", "toff"]:
-            self.timing_index_to_ms(name_timing, t)
+            self.index_to_ms(name_timing, t)
         self.tend = len(t) * self.dt
 
     def generate(self):

bluepyefe/ecode/negCheops.py

@@ -94,7 +94,7 @@ def interpret(self, t, current, config_data, reader_data):
 
         # Converting back to ms
         for name_timing in ["ton", "t1", "t2", "t3", "t4", "toff"]:
-            self.timing_index_to_ms(name_timing, t)
+            self.index_to_ms(name_timing, t)
         self.tend = len(t) * self.dt
 
     def generate(self):

bluepyefe/ecode/posCheops.py

@@ -99,7 +99,7 @@ def interpret(self, t, current, config_data, reader_data):
 
         # Converting back to ms
         for name_timing in ["ton", "t1", "t2", "t3", "t4", "toff"]:
-            self.timing_index_to_ms(name_timing, t)
+            self.index_to_ms(name_timing, t)
         self.tend = len(t) * self.dt
 
     def generate(self):

bluepyefe/ecode/ramp.py

@@ -89,7 +89,7 @@ def interpret(self, t, current, config_data, reader_data):
 
         # Converting back to ms
         for name_timing in ["ton", "toff"]:
-            self.timing_index_to_ms(name_timing, t)
+            self.index_to_ms(name_timing, t)
         self.tend = len(t) * self.dt
 
     def generate(self):

bluepyefe/ecode/sAHP.py

@@ -111,7 +111,7 @@ def interpret(self, t, current, config_data, reader_data):
 
         # Converting back to ms
         for name_timing in ["ton", "tmid", "tmid2", "toff"]:
-            self.timing_index_to_ms(name_timing, t)
+            self.index_to_ms(name_timing, t)
         self.tend = len(t) * self.dt
 
     def generate(self):

bluepyefe/ecode/sineSpec.py

@@ -105,7 +105,7 @@ def interpret(self, t, current, config_data, reader_data):
 
         # Converting back to ms
         for name_timing in ["ton", "toff"]:
-            self.timing_index_to_ms(name_timing, t)
+            self.index_to_ms(name_timing, t)
         self.tend = len(t) * self.dt
 
     def generate(self):

bluepyefe/protocol.py

@@ -130,6 +130,9 @@ def append(self, recording):
                     recording.files
                 )
 
+        if recording.auto_threshold is not None:
+            self.feature_targets[i]._auto_thresholds.append(recording.auto_threshold)
+
         self.recordings.append(recording)
 
     def as_dict(self):

bluepyefe/recording.py

@@ -1,3 +1,5 @@
+from __future__ import annotations
+
 """Recording class"""
 
 """
@@ -60,7 +62,6 @@ def __init__(self, config_data, reader_data, protocol_name):
 
         self.location = None
         self.efeatures = {}
-        self.spikecount = None
 
         self.t = None
         self.current = None
@@ -74,6 +75,8 @@ def __init__(self, config_data, reader_data, protocol_name):
             )
 
         self.export_attr = None
+        self.auto_threshold = None
+        self.peak_time = None
 
     @property
     def time(self):
@@ -85,6 +88,13 @@ def time(self, value):
         """Setter for an alias of the time attribute"""
         self.t = value
 
+    @property
+    def spikecount(self) -> int | None:
+        if self.peak_time is None:
+            return None
+        else:
+            return len(self.peak_time)
+
     def set_timing_ecode(self, name_timings, config_data):
         """Used by some of the children classes to check that the timing of
         the ecode is provided by the user and assign it to the correct
@@ -111,12 +121,15 @@ def set_amplitudes_ecode(self, amp_name, config_data, reader_data, value):
         else:
             setattr(self, amp_name, value)
 
-    def timing_index_to_ms(self, name_timing, time_series):
+    def index_to_ms(self, name_timing, time_series):
         """Used by some of the children classes to translate a timing attribute
          from index to ms."""
 
         setattr(self, name_timing, time_series[int(round(getattr(self, name_timing)))])
 
+    def ms_to_index(self, timing):
+        return int(round(timing / self.dt))
+
     def get_params(self):
         """Returns the eCode parameters"""
         return {attr: getattr(self, attr) for attr in self.export_attr}
@@ -212,6 +225,12 @@ def call_efel(self, efeatures, efel_settings=None):
             efel_settings = {}
 
         settings = {"stimulus_current": self.amp}
+
+        if "Threshold" not in settings and self.auto_threshold is not None:
+            logger.warning(f"Threshold was not provided and was automatically"
+                           f" set to {self.auto_threshold}")
+            settings["Threshold"] = self.auto_threshold
+
         for setting in efel_settings:
             if setting not in ['stim_start', 'stim_end']:
                 settings[setting] = efel_settings[setting]
@@ -249,7 +268,6 @@ def compute_efeatures(
             efel_settings (dict): eFEL settings in the form
                 {setting_name: setting_value}.
         """
-
         if efeature_names is None:
             efeature_names = efeatures
 
@@ -266,16 +284,34 @@ def compute_efeatures(
 
             self.efeatures[efeature_name] = numpy.nanmean(value)
 
-    def compute_spikecount(self, efel_settings=None):
+    def compute_spikecount(self, efel_settings=None, offset_voltage=20.):
         """Compute the number of spikes in the trace"""
         if not efel_settings:
             efel_settings = {}
 
         tmp_settings = {'strict_stiminterval': True}
-        tmp_settings.update(efel_settings)
-        efel_vals = self.call_efel(['peak_time'], tmp_settings)
+        if efel_settings is not None:
+            tmp_settings.update(efel_settings)
+
+        if tmp_settings.get("Threshold", None) is None:
+            idx_ton = self.ms_to_index(self.ton)
+            idx_toff = self.ms_to_index(self.toff)
+            step_voltage = numpy.median(self.voltage[idx_ton:idx_toff])
+            thresh = step_voltage + offset_voltage
+            tmp_settings["Threshold"] = thresh
+            efel_vals = self.call_efel(['peak_time'], tmp_settings)
+            self.peak_time = efel_vals[0]['peak_time']
+            self.auto_threshold = thresh
+
+        else:
+            efel_vals = self.call_efel(['peak_time'], tmp_settings)
+            self.peak_time = efel_vals[0]['peak_time']
 
-        self.spikecount = len(efel_vals[0]['peak_time'])
+            if self.spikecount == 0 and numpy.max(self.voltage) > tmp_settings["Threshold"]:
+                logger.warning(
+                    f"No spikes were detected in recording {self.files} but the "
+                    "voltage goes higher than the spike detection threshold."
+                )
 
     def get_plot_amplitude_title(self):
         return " ({:.01f}%)".format(self.amp_rel)
@@ -307,6 +343,13 @@ def plot(
         axis_current.plot(gen_t, gen_i, c="C1", ls="--")
         axis_voltage.plot(self.t, self.voltage, c="C0")
 
+        if self.peak_time is not None:
+            max_v = numpy.max(self.voltage)
+            for pt in self.peak_time:
+                axis_voltage.plot(
+                    [pt, pt], [max_v - 30., max_v + 10.], c="C3", ls="--", lw=0.7
+                )
+
         if display_xlabel:
             axis_voltage.set_xlabel("Time (ms)")
         if display_ylabel:

bluepyefe/target.py

@@ -72,6 +72,7 @@ def __init__(
 
         self._values = []
         self._files = []
+        self._auto_thresholds = []
 
     @property
     def mean(self):
@@ -119,9 +120,18 @@ def clear(self):
         self._values = []
         self._files = []
 
+    def add_effective_threshold(self):
+        """If auto threshold detection was used during feature extraction,
+        update the efel settings with the Threshold that was actually used"""
+
+        if self._auto_thresholds:
+            self.efel_settings["Threshold"] = numpy.median(self._auto_thresholds)
+
     def as_dict(self):
         """Returns the target in the form of a dictionary"""
 
+        self.add_effective_threshold()
+
         return {
             "efeature_name": self.efeature_name,
             "feature": self.efel_feature_name,
@@ -137,6 +147,8 @@ def as_dict(self):
     def as_legacy_dict(self, save_files_used=False):
         """Returns the target in the form of a dictionary in a legacy format"""
 
+        self.add_effective_threshold()
+
         std = self.std
         if std == 0.0:
             logger.warning(

bluepyefe/tools.py

@@ -25,7 +25,7 @@
 
 DEFAULT_EFEL_SETTINGS = {
     'strict_stiminterval': True,
-    'Threshold': -30.,
+    'Threshold': -20.,
     'interp_step': 0.025
 }
 

tests/exp_data/hippocampus-portal/data-provenance.txt

@@ -0,0 +1,7 @@
+Hippocampus Hub, (2021) a website and data portal developed and operated jointly by the Institute of Biophysics, National Research Council (Consiglio Nazionale delle Ricerche [CNR]), Italy
+ and the Blue Brain Project (BBP), École polytechnique fédérale de Lausanne (EPFL), Switzerland.
+https://www.hippocampushub.eu
+
+Applicable terms and conditions:
+Consiglio Nazionale delle Ricerche, Istituto di Biofisica — https://www.hippocampushub.eu/build/about#terms_and_conditions
+EPFL/ Blue Brain Project (BBP) —  https://www.hippocampushub.eu/model/terms-of-use/

tests/test_recording.py

@@ -2,6 +2,9 @@
 
 import unittest
 
+from numpy.testing import assert_array_almost_equal
+
+import bluepyefe.cell
 import bluepyefe.recording
 from bluepyefe.reader import igor_reader
 from bluepyefe.ecode.step import Step
@@ -52,5 +55,23 @@ def test_in_target(self):
         self.assertFalse(self.recording.in_target(90, 2))
 
 
+def test_auto_threshold_detection():
+    """Test the auto_threshold detection in Recording.compute_spikecount."""
+    cell = bluepyefe.cell.Cell(name="MouseNeuron")
+    file_metadata = {
+                    "filepath": "./tests/exp_data/hippocampus-portal/99111002.nwb",
+                    "i_unit": "A",
+                    "v_unit": "V",
+                    "t_unit": "s",
+                    "ljp": 0.0,
+                    "protocol_name": "Step",
+                }
+
+    cell.read_recordings(protocol_data=[file_metadata], protocol_name="Step")
+    cell.recordings[1].compute_spikecount()
+    assert cell.recordings[1].spikecount == 2
+    assert_array_almost_equal(cell.recordings[1].peak_time, [85.4, 346.1])
+
+
 if __name__ == "__main__":
     unittest.main()