New example: voltage clamp (#408)

---------

Co-authored-by: Jaquier Aurélien Tristan <aurelien.jaquier@epfl.ch>

docs/examples_to_rst.sh

@@ -9,13 +9,16 @@ rm docs/source/load_nwb.rst
 rm -rf docs/source/load_nwb_files
 rm docs/source/extrafeats_example.rst
 rm docs/source/multiprocessing_example.rst
+rm docs/source/voltage_clamp.rst
+rm -rf docs/source/voltage_clamp_files
 
 # convert
 jupyter nbconvert --to rst examples/sonata-network/sonata-network.ipynb
 jupyter nbconvert --to rst examples/nmc-portal/L5TTPC2.ipynb
 jupyter nbconvert --to rst examples/neo/load_nwb.ipynb
 jupyter nbconvert --to rst examples/extracellular/extrafeats_example.ipynb
 jupyter nbconvert --to rst examples/parallel/multiprocessing_example.ipynb
+jupyter nbconvert --to rst examples/voltage_clamp/voltage_clamp.ipynb
 
 # move
 mv examples/sonata-network/sonata-network.rst docs/source/
@@ -25,4 +28,6 @@ mv examples/nmc-portal/L5TTPC2_files docs/source/
 mv examples/neo/load_nwb.rst docs/source/
 mv examples/neo/load_nwb_files docs/source/
 mv examples/extracellular/extrafeats_example.rst docs/source/
-mv examples/parallel/multiprocessing_example.rst docs/source/
+mv examples/parallel/multiprocessing_example.rst docs/source/
+mv examples/voltage_clamp/voltage_clamp.rst docs/source/
+mv examples/voltage_clamp/voltage_clamp_files docs/source/

docs/source/eFeatures.rst

@@ -1710,7 +1710,7 @@ time_constant
 The extraction of the time constant requires a voltage trace of a cell in a hyper- polarized state.
 Starting at stim start find the beginning of the exponential decay where the first derivative of V(t) is smaller than -0.005 V/s in 5 subsequent points.
 The flat subsequent to the exponential decay is defined as the point where the first derivative of the voltage trace is bigger than -0.005
-and the mean of the follwowing 70 points as well.
+and the mean of the follwowing 70 ms as well.
 If the voltage trace between the beginning of the decay and the flat includes more than 9 points, fit an exponential decay.
 Yield the time constant of that decay.
 
@@ -1802,16 +1802,19 @@ decay_time_constant_after_stim
 - **Units**: ms
 - **Pseudocode**: ::
 
-    time_interval = t[numpy.where(t => decay_start_after_stim &
-                       t < decay_end_after_stim)] - t[numpy.where(t == stim_end)]
-    voltage_interval = abs(voltages[numpy.where(t => decay_start_after_stim &
-                                    t < decay_end_after_stim)]
-                           - voltages[numpy.where(t == decay_start_after_stim)])
+    interval_indices = numpy.where(
+        (time >= interval_start) & (time < interval_end))
+    stim_start_index = get_index(time, stim_start)
+    interval_time = time[interval_indices] - stim_end
+    interval_voltage = abs(
+        voltage[interval_indices] -
+        voltage[stim_start_index])
 
-    log_voltage_interval = numpy.log(voltage_interval)
-    slope, _ = numpy.polyfit(time_interval, log_voltage_interval, 1)
+    # fit
+    log_interval_voltage = numpy.log(interval_voltage)
+    slope, _ = numpy.polyfit(interval_time, log_interval_voltage, 1)
 
-    decay_time_constant_after_stim = -1. / slope
+    tau = -1. / slope
 
 multiple_decay_time_constant_after_stim
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -2085,6 +2088,126 @@ with impedance_max_freq being a setting with 50.0 as a default value.
     else:
         return None
 
+activation_time_constant
+~~~~~~~~~~~~~~~~~~~~~~~~
+
+`Python efeature`_ : Time constant for an ion channel activation trace.
+Fits for stim_start to trace maximum interval as A - B * exp(-t/tau).
+
+**Attention!** For *voltage clamp* data, user should pass the current response as voltage to efel.
+See voltage clamp example for more details.
+
+- **Required features**: time, voltage, stim_start, stim_end
+- **Units**: ms
+- **Pseudocode**: ::
+
+    def exp_fit(t, tau, A0, A1) -> np.ndarray | float:
+        return A0 + A1 * np.exp(-t / tau)
+
+    # isolate stimulus interval
+    stim_start_idx = np.flatnonzero(time >= stim_start)[0]
+    stim_end_idx = np.flatnonzero(time >= stim_end)[0]
+    time_interval = time[stim_start_idx:stim_end_idx]
+    voltage_interval = voltage[stim_start_idx:stim_end_idx]
+
+    # keep trace going from stim_start to voltage max
+    max_idx = np.argmax(voltage_interval)
+    time_interval = time_interval[:max_idx + 1]
+    voltage_interval = voltage_interval[:max_idx + 1]
+
+    # correct time so that it starts from 0
+    time_interval -= time_interval[0]
+
+    # fit
+    popt, _ = curve_fit(
+        exp_fit,
+        time_interval,
+        voltage_interval,
+        p0=(1., voltage_interval[-1], voltage_interval[0] - voltage_interval[-1]),
+        bounds=((0, -np.inf, -np.inf), (np.inf, np.inf, 0)),  # positive tau, negative A1
+        nan_policy="omit",
+    )
+    time_constant =  np.array([abs(popt[0])])
+
+deactivation_time_constant
+~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+`Python efeature`_ : Time constant for an ion channel deactivation trace.
+Fits for stim_start to stim_end as A + B * exp(-t/tau).
+
+**Attention!** For *voltage clamp* data, user should pass the current response as voltage to efel.
+See voltage clamp example for more details.
+
+- **Required features**: time, voltage, stim_start, stim_end
+- **Units**: ms
+- **Pseudocode**: ::
+
+    def exp_fit(t, tau, A0, A1) -> np.ndarray | float:
+        return A0 + A1 * np.exp(-t / tau)
+
+    # isolate stimulus interval
+    interval_indices = np.where((time >= stim_start) & (time < stim_end))
+    time_interval = time[interval_indices]
+    voltage_interval = voltage[interval_indices]
+
+    # correct time so that it starts from 0
+    time_interval -= time_interval[0]
+
+    # fit
+    popt, _ = curve_fit(
+        exp_fit,
+        time_interval,
+        voltage_interval,
+        p0=(1., voltage_interval[-1], max(0, voltage_interval[0] - voltage_interval[-1])),
+        bounds=((0, -np.inf, 0), np.inf),  # positive tau, positive A1
+        nan_policy="omit",
+    )
+    time_constant =  np.array([abs(popt[0])])
+
+inactivation_time_constant
+~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+`Python efeature`_ : Time constant for an ion channel inactivation trace.
+Fits for trace maximum to stim end interval as A + B * exp(-t/tau).
+Depends on inactivation_tc_end_skip setting, which removes a given number of data points at the end of the trace,
+right before stim_end. This is useful to remove artifacts that would bias the fit. Default is 10 data points.
+
+**Attention!** For *voltage clamp* data, user should pass the current response as voltage to efel.
+See voltage clamp example for more details.
+
+- **Required features**: time, voltage, stim_start, stim_end, inactivation_tc_end_skip (default = 10)
+- **Units**: ms
+- **Pseudocode**: ::
+
+    def exp_fit(t, tau, A0, A1) -> np.ndarray | float:
+        return A0 + A1 * np.exp(-t / tau)
+
+    # isolate stimulus interval
+    stim_start_idx = np.flatnonzero(time >= stim_start)[0]
+    stim_end_idx = np.flatnonzero(time >= stim_end)[0]
+    time_interval = time[stim_start_idx:stim_end_idx - end_skip]
+    voltage_interval = voltage[stim_start_idx:stim_end_idx - end_skip]
+
+    # keep trace going from voltage max to stim end
+    # remove end of trace to remove artifacts due to stimulus change
+    max_idx = np.argmax(voltage_interval)
+    time_interval = time_interval[max_idx:]
+    voltage_interval = voltage_interval[max_idx:]
+
+    # correct time so that it starts from 0
+    time_interval -= time_interval[0]
+
+    # fit
+    popt, _ = curve_fit(
+        exp_fit,
+        time_interval,
+        voltage_interval,
+        p0=(1., voltage_interval[-1], voltage_interval[0] - voltage_interval[-1]),
+        bounds=((0, -np.inf, 0), np.inf),  # positive tau, positive A1
+        nan_policy="omit",
+    )
+    time_constant = np.array([abs(popt[0])])
+
 Extracellular features
 ----------------------
 

docs/source/examples.rst

@@ -10,4 +10,5 @@ Examples
    load_nwb
    nmc-portal
    sonata-network
-   extrafeats_example
+   extrafeats_example
+   voltage_clamp

efel/pyfeatures/pyfeatures.py

@@ -71,7 +71,10 @@
     'inv_third_ISI',
     'inv_fourth_ISI',
     'inv_fifth_ISI',
-    'inv_last_ISI'
+    'inv_last_ISI',
+    'activation_time_constant',
+    'deactivation_time_constant',
+    'inactivation_time_constant',
 ]
 
 
@@ -352,3 +355,148 @@ def phaseslope_max() -> np.ndarray | None:
         return np.array([np.max(phaseslope)])
     except ValueError:
         return None
+
+
+def exp_fit(t, tau, A0, A1) -> np.ndarray | float:
+    """Exponential function used in exponential fitting.
+
+    Args:
+        t (ndarray or float): time series
+        tau (float): time constant
+        A0 (float): constant added to the exponential
+        A1 (float): constant multiplying the exponential
+    """
+    return A0 + A1 * np.exp(-t / tau)
+
+
+def activation_time_constant() -> np.ndarray | None:
+    """Time constant for an ion channel activation trace.
+    Fits for stim_start to trace maximum interval as A - B * exp(-t/tau)."""
+    from scipy.optimize import curve_fit
+
+    stim_start = _get_cpp_data("stim_start")
+    stim_end = _get_cpp_data("stim_end")
+    voltage = get_cpp_feature("voltage")
+    time = get_cpp_feature("time")
+
+    if voltage is None or time is None:
+        return None
+
+    # isolate stimulus interval
+    stim_start_idx = np.flatnonzero(time >= stim_start)[0]
+    stim_end_idx = np.flatnonzero(time >= stim_end)[0]
+    time_interval = time[stim_start_idx:stim_end_idx]
+    voltage_interval = voltage[stim_start_idx:stim_end_idx]
+
+    # keep trace going from stim_start to voltage max
+    max_idx = np.argmax(voltage_interval)
+    time_interval = time_interval[:max_idx + 1]
+    voltage_interval = voltage_interval[:max_idx + 1]
+
+    # correct time so that it starts from 0
+    time_interval -= time_interval[0]
+
+    # fit
+    try:
+        popt, _ = curve_fit(
+            exp_fit,
+            time_interval,
+            voltage_interval,
+            p0=(1., voltage_interval[-1], voltage_interval[0] - voltage_interval[-1]),
+            # positive tau, negative A1
+            bounds=((0, -np.inf, -np.inf), (np.inf, np.inf, 0)),
+            nan_policy="omit",
+        )
+    except (ValueError, RuntimeError):
+        return None
+
+    return np.array([abs(popt[0])])
+
+
+def deactivation_time_constant() -> np.ndarray | None:
+    """Time constant for an ion channel deactivation trace.
+    Fits for stim_start to stim_end as A + B * exp(-t/tau)."""
+    from scipy.optimize import curve_fit
+
+    stim_start = _get_cpp_data("stim_start")
+    stim_end = _get_cpp_data("stim_end")
+    voltage = get_cpp_feature("voltage")
+    time = get_cpp_feature("time")
+
+    if voltage is None or time is None:
+        return None
+
+    # isolate stimulus interval
+    interval_indices = np.where((time >= stim_start) & (time < stim_end))
+    time_interval = time[interval_indices]
+    voltage_interval = voltage[interval_indices]
+
+    # correct time so that it starts from 0
+    time_interval -= time_interval[0]
+
+    # fit
+    try:
+        popt, _ = curve_fit(
+            exp_fit,
+            time_interval,
+            voltage_interval,
+            p0=(
+                1., voltage_interval[-1], max(
+                    0, voltage_interval[0] - voltage_interval[-1]
+                )
+            ),
+            bounds=((0, -np.inf, 0), np.inf),  # positive tau, positive A1
+            nan_policy="omit",
+        )
+    except (ValueError, RuntimeError):
+        return None
+
+    return np.array([abs(popt[0])])
+
+
+def inactivation_time_constant() -> np.ndarray | None:
+    """Time constant for an ion channel inactivation trace.
+    Fits for trace maximum to stim end interval as A + B * exp(-t/tau)."""
+    from scipy.optimize import curve_fit
+
+    stim_start = _get_cpp_data("stim_start")
+    stim_end = _get_cpp_data("stim_end")
+    voltage = get_cpp_feature("voltage")
+    time = get_cpp_feature("time")
+    # used to remove end of trace to remove artifacts due to stimulus change
+    end_skip = _get_cpp_data("inactivation_tc_end_skip")
+
+    if voltage is None or time is None:
+        return None
+
+    # isolate stimulus interval
+    stim_start_idx = np.flatnonzero(time >= stim_start)[0]
+    stim_end_idx = np.flatnonzero(time >= stim_end)[0]
+    time_interval = time[stim_start_idx:stim_end_idx - end_skip]
+    voltage_interval = voltage[stim_start_idx:stim_end_idx - end_skip]
+
+    # keep trace going from voltage max to stim end
+    # remove end of trace to remove artifacts due to stimulus change
+    max_idx = np.argmax(voltage_interval)
+    time_interval = time_interval[max_idx:]
+    voltage_interval = voltage_interval[max_idx:]
+
+    # correct time so that it starts from 0
+    if time_interval.size < 1:
+        return None
+    time_interval -= time_interval[0]
+
+    # fit
+    try:
+        popt, _ = curve_fit(
+            exp_fit,
+            time_interval,
+            voltage_interval,
+            p0=(1., voltage_interval[-1], voltage_interval[0] - voltage_interval[-1]),
+            bounds=((0, -np.inf, 0), np.inf),  # positive tau, positive A1
+            nan_policy="omit",
+        )
+    except (ValueError, RuntimeError):
+        return None
+
+    return np.array([abs(popt[0])])

efel/settings.py

@@ -65,6 +65,8 @@ class Settings:
         sahp_start (float): SAHP start (default: 5.0).
         ignore_first_ISI (bool): Ignore first ISI (default: True).
         impedance_max_freq (float): Impedance maximum frequency (default: 50.0).
+        inactivation_tc_end_skip (int): number of data points to skip before
+            stim end for inactivation_time_constant feature
     """
 
     Threshold: float = -20.0
@@ -99,6 +101,7 @@ class Settings:
     ignore_first_ISI: bool = True
     impedance_max_freq: float = 50.0
     AP_phaseslope_range: int = 2
+    inactivation_tc_end_skip: int = 10
 
     def set_setting(self,
                     setting_name: str,

efel/units/units.json

@@ -156,5 +156,8 @@
     "neg_peak_diff": "s",
     "pos_peak_diff": "s",
     "neg_image": "constant",
-    "pos_image": "constant"
+    "pos_image": "constant",
+    "activation_time_constant": "ms",
+    "deactivation_time_constant": "ms",
+    "inactivation_time_constant": "ms"
 }