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}
 }