From 9f9ad417d2ae4ebd23aeb2707a3c4a9ba2ed2a27 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Aur=C3=A9lien=20Jaquier?= <72930209+AurelienJaquier@users.noreply.github.com> Date: Thu, 11 Jul 2024 17:14:37 +0200 Subject: [PATCH] fix impedance (#402) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-authored-by: Jaquier Aurélien Tristan --- CHANGELOG.rst | 9 ++++++++ docs/source/eFeatures.rst | 14 +++++++++++++ efel/DependencyV5.txt | 1 + efel/cppcore/BasicFeatures.cpp | 2 +- efel/cppcore/FillFptrTable.cpp | 1 + efel/cppcore/Subthreshold.cpp | 38 ++++++++++++++++++++++++++++++++-- efel/cppcore/Subthreshold.h | 3 +++ efel/cppcore/cfeature.cpp | 1 + efel/pyfeatures/pyfeatures.py | 10 +++++++-- efel/units/units.json | 1 + tests/featurenames.json | 1 + tests/test_allfeatures.py | 2 +- tests/test_basic.py | 29 ++++++++++++++++++++++++++ tests/test_pyfeatures.py | 2 +- 14 files changed, 107 insertions(+), 7 deletions(-) diff --git a/CHANGELOG.rst b/CHANGELOG.rst index 96278fc1..63fe2b62 100644 --- a/CHANGELOG.rst +++ b/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 `_, and this project adheres to `Semantic Versioning `_. +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 ---------------- diff --git a/docs/source/eFeatures.rst b/docs/source/eFeatures.rst index 3f5ecaa7..395f7821 100644 --- a/docs/source/eFeatures.rst +++ b/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 ~~~~~~~~~~~~~~~~~~ diff --git a/efel/DependencyV5.txt b/efel/DependencyV5.txt index 18dbbd68..63dc610a 100644 --- a/efel/DependencyV5.txt +++ b/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 diff --git a/efel/cppcore/BasicFeatures.cpp b/efel/cppcore/BasicFeatures.cpp index f61dc171..b6d4cd31 100644 --- a/efel/cppcore/BasicFeatures.cpp +++ b/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; diff --git a/efel/cppcore/FillFptrTable.cpp b/efel/cppcore/FillFptrTable.cpp index d3ab0967..67c53eff 100644 --- a/efel/cppcore/FillFptrTable.cpp +++ b/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; diff --git a/efel/cppcore/Subthreshold.cpp b/efel/cppcore/Subthreshold.cpp index bbeff88f..e2e9fd0c 100644 --- a/efel/cppcore/Subthreshold.cpp +++ b/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 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& currents = doubleFeatures.at("I"); + const vector& 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 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 diff --git a/efel/cppcore/Subthreshold.h b/efel/cppcore/Subthreshold.h index 4e3a00bc..4ee8adc9 100644 --- a/efel/cppcore/Subthreshold.h +++ b/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); diff --git a/efel/cppcore/cfeature.cpp b/efel/cppcore/cfeature.cpp index 36fd7845..a2d0af0a 100644 --- a/efel/cppcore/cfeature.cpp +++ b/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"; diff --git a/efel/pyfeatures/pyfeatures.py b/efel/pyfeatures/pyfeatures.py index 96d0bdee..83738e85 100644 --- a/efel/pyfeatures/pyfeatures.py +++ b/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) diff --git a/efel/units/units.json b/efel/units/units.json index c584bf10..ba48c410 100644 --- a/efel/units/units.json +++ b/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", diff --git a/tests/featurenames.json b/tests/featurenames.json index f97a0514..d69c648c 100644 --- a/tests/featurenames.json +++ b/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", diff --git a/tests/test_allfeatures.py b/tests/test_allfeatures.py index 21c414bb..bb61b2de 100644 --- a/tests/test_allfeatures.py +++ b/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(): diff --git a/tests/test_basic.py b/tests/test_basic.py index d139d2d1..9a712604 100644 --- a/tests/test_basic.py +++ b/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""" diff --git a/tests/test_pyfeatures.py b/tests/test_pyfeatures.py index a7f0e8cf..39711e86 100644 --- a/tests/test_pyfeatures.py +++ b/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} }