Fix AP_begin_indices-dependent features #minor

CHANGELOG.md

@@ -5,6 +5,13 @@ 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.5.0] - 2024-01
+
+### C++ changes
+- AP_end_indices, AP_rise_time, AP_fall_time, AP_rise_rate, AP_fall_rate do not take into account peaks before stim_start anymore
+- New test and test data for spontaneous firing case.
+  The data is provided by github user SzaBoglarka using cell https://modeldb.science/114047
+
 
 ## [5.4.0] - 2024-01
 

efel/DependencyV5.txt

@@ -34,7 +34,7 @@ LibV5:AP_end_indices        #LibV5:peak_indices #LibV1:interpolate
 LibV2:AP_fall_indices       #LibV5:peak_indices     #LibV5:AP_begin_indices     #LibV5:AP_end_indices #LibV1:interpolate 
 LibV2:AP_duration	        #LibV5:AP_begin_indices #LibV5:AP_end_indices #LibV1:interpolate 
 LibV2:AP_duration_half_width	#LibV2:AP_rise_indices	#LibV2:AP_fall_indices #LibV1:interpolate 
-LibV2:AP_rise_time	#LibV5:AP_begin_indices	#LibV5:peak_indices #LibV1:AP_amplitude #LibV1:interpolate 
+LibV2:AP_rise_time	#LibV5:AP_begin_indices	#LibV5:peak_indices #LibV1:interpolate 
 LibV2:AP_fall_time	#LibV5:peak_indices	#LibV5:AP_end_indices #LibV1:interpolate 
 LibV2:AP_rise_rate	#LibV5:AP_begin_indices	#LibV5:peak_indices #LibV1:interpolate 
 LibV2:AP_fall_rate	#LibV5:peak_indices	#LibV5:AP_end_indices #LibV1:interpolate 

efel/cppcore/LibV2.cpp

@@ -168,28 +168,38 @@ int LibV2::AP_duration_half_width(mapStr2intVec& IntFeatureData,
 static int __AP_rise_time(const vector<double>& t, const vector<double>& v,
                           const vector<int>& apbeginindices,
                           const vector<int>& peakindices,
-                          const vector<double>& apamplitude, double beginperc,
-                          double endperc, vector<double>& aprisetime) {
+                          double beginperc, double endperc,
+                          vector<double>& aprisetime) {
   aprisetime.resize(std::min(apbeginindices.size(), peakindices.size()));
+  // Make sure that we do not use peaks starting before the 1st AP_begin_index
+  // Because AP_begin_indices only takes into account peaks after stimstart
+  vector<int> newpeakindices;
+  if (apbeginindices.size() > 0) {
+    newpeakindices = peaks_after_stim_start(apbeginindices[0], peakindices);
+  }
   double begin_v;
   double end_v;
   double begin_indice;
   double end_indice;
+  double apamplitude;
   for (size_t i = 0; i < aprisetime.size(); i++) {
-    begin_v = v[apbeginindices[i]] + beginperc * apamplitude[i];
-    end_v = v[apbeginindices[i]] + endperc * apamplitude[i];
+    // do not use AP_amplitude feature because it does not take into account
+    // peaks after stim_end
+    apamplitude = v[newpeakindices[i]] - v[apbeginindices[i]];
+    begin_v = v[apbeginindices[i]] + beginperc * apamplitude;
+    end_v = v[apbeginindices[i]] + endperc * apamplitude;
 
     // Get begin indice
     size_t j = apbeginindices[i];
     // change slightly begin_v for almost equal case
     // truncature error can change begin_v even when beginperc == 0.0
-    while (j < peakindices[i] && v[j] < begin_v - 0.0000000000001) {
+    while (j < newpeakindices[i] && v[j] < begin_v - 0.0000000000001){
       j++;
     }
     begin_indice = j;
 
     // Get end indice
-    j = peakindices[i];
+    j = newpeakindices[i];
     // change slightly end_v for almost equal case
     // truncature error can change end_v even when beginperc == 0.0
     while (j > apbeginindices[i] && v[j] > end_v + 0.0000000000001) {
@@ -207,14 +217,13 @@ int LibV2::AP_rise_time(mapStr2intVec& IntFeatureData,
   // Fetching all required features in one go.
   const auto& doubleFeatures = getFeatures(
       DoubleFeatureData,
-      {"T", "V", "AP_amplitude", "rise_start_perc", "rise_end_perc"});
+      {"T", "V", "rise_start_perc", "rise_end_perc"});
   const auto& intFeatures =
       getFeatures(IntFeatureData, {"AP_begin_indices", "peak_indices"});
   vector<double> aprisetime;
   int retval = __AP_rise_time(
       doubleFeatures.at("T"), doubleFeatures.at("V"),
       intFeatures.at("AP_begin_indices"), intFeatures.at("peak_indices"),
-      doubleFeatures.at("AP_amplitude"),
       doubleFeatures.at("rise_start_perc").empty()
           ? 0.0
           : doubleFeatures.at("rise_start_perc").front(),
@@ -230,28 +239,34 @@ int LibV2::AP_rise_time(mapStr2intVec& IntFeatureData,
 
 // *** AP_fall_time according to E10 and E18 ***
 static int __AP_fall_time(const vector<double>& t,
+                          const double stimstart,
                           const vector<int>& peakindices,
                           const vector<int>& apendindices,
                           vector<double>& apfalltime) {
   apfalltime.resize(std::min(peakindices.size(), apendindices.size()));
+  // Make sure that we do not use peaks starting before stim start
+  // Because AP_end_indices only takes into account peaks after stim start
+  vector<int> newpeakindices = peaks_after_stim_start(stimstart, peakindices, t);
+
   for (size_t i = 0; i < apfalltime.size(); i++) {
-    apfalltime[i] = t[apendindices[i]] - t[peakindices[i]];
+    apfalltime[i] = t[apendindices[i]] - t[newpeakindices[i]];
   }
   return apfalltime.size();
 }
 int LibV2::AP_fall_time(mapStr2intVec& IntFeatureData,
                         mapStr2doubleVec& DoubleFeatureData,
                         mapStr2Str& StringData) {
-  const auto& doubleFeatures = getFeatures(DoubleFeatureData, {"T"});
+  const auto& doubleFeatures = getFeatures(DoubleFeatureData, {"T", "stim_start"});
   const auto& intFeatures =
       getFeatures(IntFeatureData, {"peak_indices", "AP_end_indices"});
 
   const vector<double>& t = doubleFeatures.at("T");
+  const double stim_start = doubleFeatures.at("stim_start")[0];
   const vector<int>& peakindices = intFeatures.at("peak_indices");
   const vector<int>& apendindices = intFeatures.at("AP_end_indices");
 
   vector<double> apfalltime;
-  int retval = __AP_fall_time(t, peakindices, apendindices, apfalltime);
+  int retval = __AP_fall_time(t, stim_start, peakindices, apendindices, apfalltime);
 
   if (retval > 0) {
     setVec(DoubleFeatureData, StringData, "AP_fall_time", apfalltime);
@@ -265,9 +280,13 @@ static int __AP_rise_rate(const vector<double>& t, const vector<double>& v,
                           const vector<int>& peakindices,
                           vector<double>& apriserate) {
   apriserate.resize(std::min(peakindices.size(), apbeginindices.size()));
+  vector<int> newpeakindices;
+  if (apbeginindices.size() > 0) {
+    newpeakindices = peaks_after_stim_start(apbeginindices[0], peakindices);
+  }
   for (size_t i = 0; i < apriserate.size(); i++) {
-    apriserate[i] = (v[peakindices[i]] - v[apbeginindices[i]]) /
-                    (t[peakindices[i]] - t[apbeginindices[i]]);
+    apriserate[i] = (v[newpeakindices[i]] - v[apbeginindices[i]]) /
+                    (t[newpeakindices[i]] - t[apbeginindices[i]]);
   }
   return apriserate.size();
 }
@@ -294,30 +313,34 @@ int LibV2::AP_rise_rate(mapStr2intVec& IntFeatureData,
 
 // *** AP_fall_rate according to E12 and E20 ***
 static int __AP_fall_rate(const vector<double>& t, const vector<double>& v,
+                          const double stimstart,
                           const vector<int>& peakindices,
                           const vector<int>& apendindices,
                           vector<double>& apfallrate) {
   apfallrate.resize(std::min(apendindices.size(), peakindices.size()));
+  vector<int> newpeakindices = peaks_after_stim_start(stimstart, peakindices, t);
+
   for (size_t i = 0; i < apfallrate.size(); i++) {
-    apfallrate[i] = (v[apendindices[i]] - v[peakindices[i]]) /
-                    (t[apendindices[i]] - t[peakindices[i]]);
+    apfallrate[i] = (v[apendindices[i]] - v[newpeakindices[i]]) /
+                    (t[apendindices[i]] - t[newpeakindices[i]]);
   }
   return apfallrate.size();
 }
 int LibV2::AP_fall_rate(mapStr2intVec& IntFeatureData,
                         mapStr2doubleVec& DoubleFeatureData,
                         mapStr2Str& StringData) {
-  const auto& doubleFeatures = getFeatures(DoubleFeatureData, {"T", "V"});
+  const auto& doubleFeatures = getFeatures(DoubleFeatureData, {"T", "V", "stim_start"});
   const auto& intFeatures =
       getFeatures(IntFeatureData, {"peak_indices", "AP_end_indices"});
 
   const vector<double>& t = doubleFeatures.at("T");
   const vector<double>& v = doubleFeatures.at("V");
+  const double stim_start = doubleFeatures.at("stim_start")[0];
   const vector<int>& peakindices = intFeatures.at("peak_indices");
   const vector<int>& apendindices = intFeatures.at("AP_end_indices");
 
   vector<double> apfallrate;
-  int retval = __AP_fall_rate(t, v, peakindices, apendindices, apfallrate);
+  int retval = __AP_fall_rate(t, v, stim_start, peakindices, apendindices, apfallrate);
 
   if (retval > 0) {
     setVec(DoubleFeatureData, StringData, "AP_fall_rate", apfallrate);

efel/cppcore/LibV5.cpp

@@ -536,31 +536,40 @@ int LibV5::AP_begin_indices(mapStr2intVec& IntFeatureData,
 }
 
 static int __AP_end_indices(const vector<double>& t, const vector<double>& v,
-                            const vector<int>& pi, vector<int>& apei,
-                            double derivativethreshold) {
+                            const double stimstart, const vector<int>& pi,
+                            vector<int>& apei, double derivativethreshold) {
+
   vector<double> dvdt(v.size());
   vector<double> dv;
   vector<double> dt;
+  vector<int> peak_indices;
   int max_slope;
   getCentralDifferenceDerivative(1., v, dv);
   getCentralDifferenceDerivative(1., t, dt);
   transform(dv.begin(), dv.end(), dt.begin(), dvdt.begin(),
             std::divides<double>());
 
-  apei.resize(pi.size());
-  vector<int> picopy(pi.begin(), pi.end());
-  picopy.push_back(v.size() - 1);
+  auto stimbeginindex = distance(t.begin(),
+      find_if(t.begin(), t.end(),
+          [stimstart](double time){ return time >= stimstart; }));
+
+   for (size_t i = 0; i < pi.size(); i++) {
+    if (pi[i] > stimbeginindex) {
+      peak_indices.push_back(pi[i]);
+    }
+  }
+  peak_indices.push_back(v.size() - 1);
 
-  for (size_t i = 0; i < apei.size(); i++) {
+  for (size_t i = 0; i < peak_indices.size() - 1; i++) {
     max_slope =
-        distance(dvdt.begin(), std::min_element(dvdt.begin() + picopy[i] + 1,
-                                                dvdt.begin() + picopy[i + 1]));
+        distance(dvdt.begin(), std::min_element(dvdt.begin() + peak_indices[i] + 1,
+                                                dvdt.begin() + peak_indices[i + 1]));
     // assure that the width of the slope is bigger than 4
-    apei[i] = distance(dvdt.begin(), find_if(dvdt.begin() + max_slope,
-                                             dvdt.begin() + picopy[i + 1],
+    apei.push_back(distance(dvdt.begin(), find_if(dvdt.begin() + max_slope,
+                                             dvdt.begin() + peak_indices[i + 1],
                                              [derivativethreshold](double x) {
                                                return x >= derivativethreshold;
-                                             }));
+                                             })));
   }
   return apei.size();
 }
@@ -570,14 +579,15 @@ int LibV5::AP_end_indices(mapStr2intVec& IntFeatureData,
                           mapStr2Str& StringData) {
   const auto& T = getFeature(DoubleFeatureData, "T");
   const auto& V = getFeature(DoubleFeatureData, "V");
+  const auto& stim_start = getFeature(DoubleFeatureData, "stim_start");
   const auto& peak_indices = getFeature(IntFeatureData, "peak_indices");
 
   vector<double> dTh;
   int retVal = getParam(DoubleFeatureData, "DownDerivativeThreshold", dTh);
   double downDerivativeThreshold = (retVal <= 0) ? -12.0 : dTh[0];
 
   vector<int> AP_end_indices;
-  retVal = __AP_end_indices(T, V, peak_indices, AP_end_indices,
+  retVal = __AP_end_indices(T, V, stim_start[0], peak_indices, AP_end_indices,
                             downDerivativeThreshold);
   if (retVal >= 0) {
     setVec(IntFeatureData, StringData, "AP_end_indices", AP_end_indices);

efel/cppcore/Utils.cpp

@@ -200,6 +200,31 @@ std::pair<size_t, size_t> get_time_index(const std::vector<double>& t,
   return std::pair<size_t, size_t>(startIndex, endIndex);
 }
 
+// expects an empty newpeakindices
+vector<int> peaks_after_stim_start(const double stimstart,
+                                   const vector<int>& peakindices,
+                                   const vector<double>& t) {
+    vector<int> newpeakindices;
+    for (size_t i = 0; i < peakindices.size(); i++) {
+        if (t[peakindices[i]] > stimstart) {
+            newpeakindices.push_back(peakindices[i]);
+        }
+    }
+    return newpeakindices;
+}
+
+// expects an empty newpeakindices
+vector<int> peaks_after_stim_start(const double stimstartindex,
+                                   const vector<int>& peakindices){
+  vector<int> newpeakindices;
+  for (size_t i = 0; i < peakindices.size(); i++) {
+    if (peakindices[i] > stimstartindex) {
+      newpeakindices.push_back(peakindices[i]);
+    }
+  }
+  return newpeakindices;
+}
+
 template <class T>
 double vec_mean(const vector<T>& v) {
   /*

efel/cppcore/Utils.h

@@ -59,7 +59,14 @@ double vec_mean(const vector<T> &v);
 
 std::pair<size_t, size_t> get_time_index(const std::vector<double> &t, double startTime,
                                      double endTime, double precisionThreshold);
-
+
+vector<int> peaks_after_stim_start(const double stimstart,
+                           const vector<int>& peakindices,
+                           const vector<double>& t);
+
+vector<int> peaks_after_stim_start(const double stimstartindex,
+                           const vector<int>& peakindices);
+
 template <class ForwardIterator>
 ForwardIterator first_min_element(ForwardIterator first, ForwardIterator last) {
   ForwardIterator lowest = first;