fix impedance (#402)

Co-authored-by: Jaquier Aurélien Tristan <[email protected]>

CHANGELOG.rst

@@ -5,6 +5,15 @@ All notable changes to this project will be documented in this file.
 The format is based on `Keep a Changelog <https://keepachangelog.com/en/1.0.0/>`_,
 and this project adheres to `Semantic Versioning <https://semver.org/spec/v2.0.0.html>`_.
 
+5.7.1 - 2024-07
+---------------
+
+- When stimulus starts at 0 for impedance (which can happen often), use steady_state_voltage_stimend to get holding voltage, and equivalent for holding current.
+  Since at the end of the stimulus, we expect the current and voltage to vary rapidly around the holding value, this is a good enough proxy to get it.
+- implemented test_steady_state_current_stimend, along with documentation and test
+- fixed current interpolation
+- fixed stim_end value in impedance test data file
+
 5.7.0 - 2024-07
 ----------------
 

docs/source/eFeatures.rst

@@ -1632,6 +1632,20 @@ steady_state_voltage_stimend
     end_time = stim_end
     steady_state_voltage_stimend = numpy.mean(voltage[numpy.where((t < end_time) & (t >= begin_time))])
 
+steady_state_current_stimend
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+`Subthreshold`_ : The average current during the last 10% of the stimulus duration.
+
+- **Required features**: t, I, stim_start, stim_end
+- **Units**: nA
+- **Pseudocode**: ::
+
+    stim_duration = stim_end - stim_start
+    begin_time = stim_end - 0.1 * stim_duration
+    end_time = stim_end
+    steady_state_current_stimend = numpy.mean(current[numpy.where((t < end_time) & (t >= begin_time))])
+
 steady_state_hyper
 ~~~~~~~~~~~~~~~~~~
 

efel/DependencyV5.txt

@@ -90,6 +90,7 @@ BasicFeatures:voltage #BasicFeatures:interpolate
 BasicFeatures:current #BasicFeatures:interpolate 
 BasicFeatures:time #BasicFeatures:interpolate 
 Subthreshold:steady_state_voltage_stimend #BasicFeatures:interpolate 
+Subthreshold:steady_state_current_stimend #BasicFeatures:interpolate 
 Subthreshold:voltage_base #BasicFeatures:interpolate
 Subthreshold:current_base #BasicFeatures:interpolate
 Subthreshold:decay_time_constant_after_stim #BasicFeatures:interpolate 

efel/cppcore/BasicFeatures.cpp

@@ -44,7 +44,7 @@ int BasicFeatures::interpolate(mapStr2intVec& IntFeatureData,
     I = getFeature(DoubleFeatureData, "I");
     LinearInterpolation(InterpStep, T, I, TIntrpolI, IIntrpol);
     setVec(DoubleFeatureData, StringData, "I", IIntrpol);
-    setVec(DoubleFeatureData, StringData, "T", TIntrpol);
+    setVec(DoubleFeatureData, StringData, "T", TIntrpolI);
   } catch (...) {
   }  // pass, it is optional
   return 1;

efel/cppcore/FillFptrTable.cpp

@@ -195,6 +195,7 @@ int FillFptrTable() {
 
   //******************  FptrTableST *****************************
   FptrTableST["steady_state_voltage_stimend"] = &Subthreshold::steady_state_voltage_stimend;
+  FptrTableST["steady_state_current_stimend"] = &Subthreshold::steady_state_current_stimend;
   FptrTableST["steady_state_hyper"] = &Subthreshold::steady_state_hyper;
   FptrTableST["steady_state_voltage"] = &Subthreshold::steady_state_voltage;
   FptrTableST["voltage_base"] = &Subthreshold::voltage_base;

efel/cppcore/Subthreshold.cpp

@@ -62,12 +62,12 @@ int Subthreshold::steady_state_voltage_stimend(mapStr2intVec& IntFeatureData,
   size_t start_index = distance(times.begin(), start_it);
   size_t stop_index = distance(times.begin(), stop_it);
 
-  double mean = accumulate(voltages.begin() + start_index,
-                           voltages.begin() + stop_index, 0.0);
   size_t mean_size = stop_index - start_index;
 
   vector<double> ssv;
   if (mean_size > 0) {
+    double mean = accumulate(voltages.begin() + start_index,
+                             voltages.begin() + stop_index, 0.0);
     mean /= mean_size;
     ssv.push_back(mean);
     setVec(DoubleFeatureData, StringData, "steady_state_voltage_stimend", ssv);
@@ -76,6 +76,40 @@ int Subthreshold::steady_state_voltage_stimend(mapStr2intVec& IntFeatureData,
   return -1;
 }
 
+// *** The average voltage during the last 90% of the stimulus duration. ***
+int Subthreshold::steady_state_current_stimend(mapStr2intVec& IntFeatureData,
+                                        mapStr2doubleVec& DoubleFeatureData,
+                                        mapStr2Str& StringData) {
+  const auto& doubleFeatures =
+      getFeatures(DoubleFeatureData, {"I", "T", "stim_end", "stim_start"});
+
+  const vector<double>& currents = doubleFeatures.at("I");
+  const vector<double>& times = doubleFeatures.at("T");
+  const double stimStart = doubleFeatures.at("stim_start")[0];
+  const double stimEnd = doubleFeatures.at("stim_end")[0];
+
+  double start_time = stimEnd - 0.1 * (stimEnd - stimStart);
+  auto start_it = find_if(times.begin(), times.end(),
+                          [start_time](double x) { return x >= start_time; });
+  auto stop_it = find_if(times.begin(), times.end(),
+                         [stimEnd](double x) { return x >= stimEnd; });
+  size_t start_index = distance(times.begin(), start_it);
+  size_t stop_index = distance(times.begin(), stop_it);
+
+  size_t mean_size = stop_index - start_index;
+
+  vector<double> ssc;
+  if (mean_size > 0) {
+    double mean = accumulate(currents.begin() + start_index,
+                             currents.begin() + stop_index, 0.0);
+    mean /= mean_size;
+    ssc.push_back(mean);
+    setVec(DoubleFeatureData, StringData, "steady_state_current_stimend", ssc);
+    return 1;
+  }
+  return -1;
+}
+
 // steady state of the voltage response during hyperpolarizing stimulus,
 // elementary feature for E29
 // *** steady_state_hyper

efel/cppcore/Subthreshold.h

@@ -28,6 +28,9 @@ namespace Subthreshold {
 int steady_state_voltage_stimend(mapStr2intVec& IntFeatureData,
                                  mapStr2doubleVec& DoubleFeatureData,
                                  mapStr2Str& StringData);
+int steady_state_current_stimend(mapStr2intVec& IntFeatureData,
+                                 mapStr2doubleVec& DoubleFeatureData,
+                                 mapStr2Str& StringData);
 int steady_state_hyper(mapStr2intVec& IntFeatureData,
                        mapStr2doubleVec& DoubleFeatureData,
                        mapStr2Str& StringData);

efel/cppcore/cfeature.cpp

@@ -176,6 +176,7 @@ void cFeature::fillfeaturetypes() {
   featuretypes["current"] = "double";
   featuretypes["time"] = "double";
   featuretypes["steady_state_voltage_stimend"] = "double";
+  featuretypes["steady_state_current_stimend"] = "double";
   featuretypes["voltage_base"] = "double";
   featuretypes["current_base"] = "double";
   featuretypes["decay_time_constant_after_stim"] = "double";

efel/pyfeatures/pyfeatures.py

@@ -145,11 +145,17 @@ def impedance():
     Z_max_freq = _get_cpp_data("impedance_max_freq")
     voltage_trace = get_cpp_feature("voltage")
     holding_voltage = get_cpp_feature("voltage_base")
-    normalized_voltage = voltage_trace - holding_voltage
+    # when stimulus starts at t=0, use steady_state_voltage_stimend as proxy
+    if holding_voltage is None:
+        holding_voltage = get_cpp_feature("steady_state_voltage_stimend")
+    normalized_voltage = voltage_trace - holding_voltage[0]
     current_trace = current()
     if current_trace is not None:
         holding_current = get_cpp_feature("current_base")
-        normalized_current = current_trace - holding_current
+        # when stimulus starts at t=0, use steady_state_current_stimend as proxy
+        if holding_current is None:
+            holding_current = get_cpp_feature("steady_state_current_stimend")
+        normalized_current = current_trace - holding_current[0]
         n_spikes = spike_count()
         if n_spikes < 1:  # if there is no spikes in ZAP
             fft_volt = np.fft.fft(normalized_voltage)

efel/units/units.json

@@ -109,6 +109,7 @@
     "steady_state_hyper": "mV",
     "steady_state_voltage": "mV",
     "steady_state_voltage_stimend": "mV",
+    "steady_state_current_stimend": "nA",
     "strict_burst_mean_freq": "Hz",
     "strict_burst_number": "constant",
     "strict_interburst_voltage": "mV",

tests/featurenames.json

@@ -120,6 +120,7 @@
     "steady_state_hyper",
     "steady_state_voltage",
     "steady_state_voltage_stimend",
+    "steady_state_current_stimend",
     "time",
     "time_constant",
     "time_to_first_spike",

tests/test_allfeatures.py

@@ -104,7 +104,7 @@ def load_data(filename):
         for x in all_featurenames
         if x
         not in db_featurenames
-        + ["current", "current_base", "impedance"]
+        + ["current", "current_base", "impedance", "steady_state_current_stimend"]
     ]
 
     with warnings.catch_warnings():

tests/test_basic.py

@@ -62,6 +62,7 @@
     / 'spiking_from_beginning_to_end.txt'
 )
 spontaneous_url = testdata_dir / 'basic' / 'spontaneous.txt'
+impedance_url = testdata_dir / 'basic' / 'impedance.txt'
 
 testdata_url = testdata_dir / 'allfeatures' / 'testdata.txt'
 
@@ -1804,6 +1805,34 @@ def test_steady_state_voltage_stimend():
     )
 
 
+def test_steady_state_current_stimend():
+    """basic: Test steady_state_current_stimend"""
+
+    import efel
+    efel.reset()
+
+    trace = {
+        'stim_start': [100.0],
+        'stim_end': [5100.0]
+    }
+    data = numpy.loadtxt(impedance_url)
+    trace['T'] = data[:, 0]
+    trace['V'] = data[:, 1]
+    trace['I'] = data[:, 2]
+
+    features = ['steady_state_current_stimend']
+
+    feature_values = \
+        get_feature_values(
+            [trace],
+            features)[0]
+
+    numpy.testing.assert_allclose(
+        -1.44246348e-10,
+        feature_values['steady_state_current_stimend'][0]
+    )
+
+
 def test_maximum_voltage_from_voltagebase():
     """basic: Test maximum_voltage_from_voltagebase"""
 

tests/test_pyfeatures.py

@@ -87,7 +87,7 @@
         'v_col': 2,
         'i_col': 3,
         'stim_start': 100.0,
-        'stim_end': 51000.0}
+        'stim_end': 5100.0}
 }