Voltage clamp trace and eFEL

This notebook explains how to use traces from voltage clamp in eFEL. First, you need to load your modules and your trace data. Here, we will use experimental data from Channelpedia: https://channelpedia.epfl.ch/expdata/details/9430

Ranjan R, Logette E, Marani M, Herzog M, Tache V, Scantamburlo E, Buchillier V and Markram H (2019) A Kinetic Map of the Homomeric Voltage-Gated Potassium Channel (Kv) Family. Front. Cell. Neurosci. 13:358. doi: 10.3389/fncel.2019.00358

Blue Brain Project Portal (https://portal.bluebrain.epfl.ch) and Channelpedia (https://channelpedia.epfl.ch)

import h5py
import numpy as np

from matplotlib import pyplot as plt
from scipy.optimize import curve_fit

import efel

Activation traces

We will start to load the data form the nwb file. In this notebook, we will use the traces from three experiments: Activation, Deactivation and Inactivation.

data_path = "../../tests/testdata/basic/rCell9430.nwb"

data = {}
exp_types = ["Activation", "Deactivation", "Inactivation"]
with h5py.File(data_path, "r") as content:
    for exp_type in exp_types:
        data[exp_type] = {}
        reps = content["acquisition"]["timeseries"][exp_type]["repetitions"]
        for rep_name, rep in reps.items():
            data[exp_type][rep_name] = {
                "dt": np.asarray(rep["x_interval"], dtype=np.float32),
                "current": np.asarray(rep["data"], dtype=np.float32),
            }

Now we will turn this data into the traces list expected by eFel. Here, we will start with the traces from the activation experiment. For this example, we will only use the 1st repetition. You can change the repetition by renaming the rep_name variable into ‘repetition2’ or ‘repetition3’. Notice that we have to give the time in ms.

Attention! eFEL was designed mainly to extract features from voltage timeseries. The voltage clamp data are currents. To use the eFEL features with voltage clamp data, we have to give the current recording to the “V” key in the trace dict. eFEL will treat the current trace as if it were a voltage trace.

traces = []
rep_name = "repetition1"
exp_type = "Activation"
for idx in range(len(data[exp_type][rep_name]["dt"])):
    trace = {}
    i = data[exp_type][rep_name]["current"][:,idx]
    t = np.arange(i.size) * data[exp_type][rep_name]["dt"][idx]
    # efel expects ms: s -> ms
    t = t * 1000.0
    trace["T"] = t
    trace["V"] = i  # trick: input current as if it was voltage
    trace["stim_start"] = [99]
    trace["stim_end"] = [600]
    traces.append(trace)

We can plot the traces to see how they look:

for t in traces:
    plt.plot(t["T"], t["V"], c="red")
plt.xlabel("Stimulus Voltage (mV)")
plt.ylabel("Response Current (nA)")
plt.title("Activation traces")

Text(0.5, 1.0, 'Activation traces')

We can now extract the desired features from eFel. Even if there is ‘voltage’ in the feature names, since we gave the current trace instead of the voltage trace, it will still compute these features for the current trace. So here, ‘maximum_voltage’ will actually compute the maximum of the current, ‘voltage_base’ will actually compute the current base and ‘steady_state_voltage_stimend’ will actually compute the current steady state at the end of the stimulus. While those three features are usually used on voltage traces, the last one, ‘activation_time_constant’ has been written specifically for this kind of voltage clamp activation trace.

feature_names = ["maximum_voltage", "voltage_base", "steady_state_voltage_stimend", "activation_time_constant"]
feats = efel.get_feature_values(traces, feature_names)

/Users/ajaquier/Documents/eFEL/efel/pyfeatures/pyfeatures.py:398: OptimizeWarning: Covariance of the parameters could not be estimated
  popt, _ = curve_fit(

From these features, we can get several information. We can, for example, extract the ‘activation voltage’, defined as the voltage at which the maximum current of the trace is at least 10% of the maximum current among all of the traces.

# Stimulus Voltage used in the experiment - used later for plotting
stim_v = list(range(-90, 90, 10))

max_i = np.array([feat_dict["maximum_voltage"][0] for feat_dict in feats])
max_i = max_i / np.max(max_i)

act_v_idx = np.argwhere(max_i >= 0.1)[0][0]
act_v = stim_v[act_v_idx]
act_i = max_i[act_v_idx]
ylim = [np.min(max_i) - 0.05, 1.05]
print(f"Activation voltage is at {act_v} mV")

Activation voltage is at -20 mV

Here, we plot the I-V curve for the 1st repetition, and we mark the activation voltage with a red line

plt.plot(stim_v, max_i, '.')
plt.xlabel("Stimulus Voltage (mV)")
plt.ylabel("Peak current / Max current")
plt.plot([act_v, act_v], [ylim[0], act_i], color="red")
plt.ylim(ylim)
plt.title("Activation voltage")

# sigmoid fit
def sigmoid(v, delta, tau):
    return 1 / (1 + np.exp(-(v - delta) / tau))

popt, _ = curve_fit(sigmoid, stim_v, max_i)
v = np.linspace(stim_v[0], stim_v[-1], 200)
i = sigmoid(v, *popt)
plt.plot(v, i, "--", c="black")

[<matplotlib.lines.Line2D at 0x1338178e0>]

We can also show the evolution of the steady state current at the end of the stimulus across input voltage:

ssis = np.array([feat_dict["steady_state_voltage_stimend"][0] for feat_dict in feats])
plt.plot(stim_v, ssis, '.')
plt.xlabel("Stimulus Voltage (mV)")
plt.ylabel("Steady State Current (nA)")
plt.title("Activation steady state")

# sigmoid fit
def sigmoid_scaled(v, delta, tau, A):
    return A / (1 + np.exp(-(v - delta) / tau))

popt, _ = curve_fit(sigmoid_scaled, stim_v, ssis)
ss = sigmoid_scaled(v, *popt)
plt.plot(v, ss, "--", c="black")

[<matplotlib.lines.Line2D at 0x1358cdb40>]

We can also plot the current base. Since this feature is computed on the trace before the stimulus is applied, it should be constant across the traces.

i_base = np.array([feat_dict["voltage_base"][0] for feat_dict in feats])
plt.plot(stim_v, i_base, '.')
plt.xlabel("Stimulus Voltage (mV)")
plt.ylabel("Current Base (nA)")
plt.title("Activation current base")

Text(0.5, 1.0, 'Activation current base')

We can see that the current base is the same for all stimulus voltages, which is expected since it is computed before the stimulus is applied.

Finally, we can get the time constant of the activation curve (from stim_start to the maximum of the trace):

act_tau_efel = np.array([feat_dict["activation_time_constant"][0] for feat_dict in feats])
plt.plot(stim_v, act_tau_efel, '.')
plt.xlabel("Stimulus Voltage (mV)")
plt.ylabel("Activation time constant (ms)")
plt.xlim((-15, 85))
plt.ylim((0, 10))
plt.title("Activation time constant")

Text(0.5, 1.0, 'Activation time constant')

As a visual verification, we can do an exponential fit to the traces, use our time constant as a parameter of the exponential and plot the exponential on top of the trace.

def exp_fit(t, tau_1, A1, A0):
        return A0 - A1 * np.exp(-t / tau_1)

for i, trace in enumerate(traces):
        stim_start_idx = np.argwhere(trace["T"] >= trace["stim_start"][0])[0][0]
        stim_end_idx = np.argwhere(trace["V"] >= np.max(trace["V"]))[0][0]
        # trace["stim_end"] - 1 to account for artefact
        if stim_end_idx < stim_start_idx or trace["T"][stim_end_idx] > trace["stim_end"][0] - 1:
                continue
        t_interval = trace["T"][stim_start_idx:stim_end_idx + 1]
        v_interval = trace["V"][stim_start_idx:stim_end_idx + 1]
        t0 = t_interval[0]
        t_interval_corrected = t_interval - t0

        popt, _ = curve_fit(exp_fit, t_interval_corrected, v_interval)

        v_fit = exp_fit(t_interval_corrected, act_tau_efel[i], popt[1], popt[2])
        plt.plot(trace["T"], trace["V"], c="red")
        plt.plot(t_interval, v_fit, "--", c="black", lw=3)
        t_tau = popt[0] + t0
        plt.plot(t_tau, exp_fit(popt[0], *popt), 'o', color='black', markersize=12)
plt.xlim(90, 170)

(90.0, 170.0)

Deactivation traces

Now we will also perform efeature extraction on traces from deactivation experiment. First, we load the data from the file. Once again, we will put the current data under the ‘V’ key representing the voltage, so that eFEL treats our trace like any voltage trace. Notice that here, we are setting the stimulus start at 401 and stimulus end at 599 because from 400 to 600 is the interval of interest for our time constant fit, and we start slightly after and end slightly before ‘real’ stimulus in order to prevent artifacts that show when there is a stimulus change that affect our computation

traces = []
rep_name = "repetition1"
exp_type = "Deactivation"
for idx in range(len(data[exp_type][rep_name]["dt"])):
    trace = {}
    i = data[exp_type][rep_name]["current"][:,idx]
    t = np.arange(i.size) * data[exp_type][rep_name]["dt"][idx]
    # efel expects ms: s -> ms
    t = t * 1000.0
    trace["T"] = t
    trace["V"] = i  # trick: input current as if it was voltage
    trace["stim_start"] = [401]
    trace["stim_end"] = [599]
    traces.append(trace)

As usual, we can start by plotting the traces:

for t in traces:
    plt.plot(t["T"], t["V"], c="red")
plt.xlabel("Stimulus Voltage (mV)")
plt.ylabel("Response Current (nA)")
plt.title("Deactivation traces")

Text(0.5, 1.0, 'Deactivation traces')

We can now extract the deactivation time constant using eFEL:

efel.reset()
stim_v = list(range(-80, 40, 10))
feature_names = ["deactivation_time_constant"]
feats = efel.get_feature_values(traces, feature_names)

deact_tau_efel = np.array([feat_dict["deactivation_time_constant"][0] for feat_dict in feats])
plt.plot(stim_v, deact_tau_efel, '.')
plt.xlabel("Stimulus Voltage (mV)")
plt.ylabel("Deactivation time constant (ms)")
plt.title("Deactivation time constant")
plt.xlim((-65, 15))
plt.ylim((-10, 200))

(-10.0, 200.0)

As a visual verification, we can do an exponential fit to the traces, use our time constant as a parameter of the exponential and plot the exponential on top of the trace.

for i, trace in enumerate(traces):
        interval_indices = np.where((trace["T"] >= trace["stim_start"][0]) & (trace["T"] < trace["stim_end"][0]))
        t_interval = trace["T"][interval_indices]
        v_interval = trace["V"][interval_indices]
        t0 = t_interval[0]
        t_interval_corrected = t_interval - t0

        popt, _ = curve_fit(exp_fit, t_interval_corrected, v_interval)

        v_fit = exp_fit(t_interval_corrected, deact_tau_efel[i], popt[1], popt[2])
        plt.plot(trace["T"], trace["V"], c="red")
        plt.plot(t_interval, v_fit, "--", c="black", lw=3)
        t_tau = popt[0] + t0
        plt.plot(t_tau, exp_fit(popt[0], *popt), 'o', color='black', markersize=12)
plt.xlim(390, 610)

/var/folders/9j/1kf4nxwd62v5vnvhcnf87h616469rk/T/ipykernel_12656/1772092254.py:2: RuntimeWarning: overflow encountered in exp
  return A0 - A1 * np.exp(-t / tau_1)

(390.0, 610.0)

Inactivation traces

Now we will also perform efeature extraction on traces from inactivation experiment. First, we load the data from the file. Once again, we will put the current data under the ‘V’ key representing the voltage, so that eFEL treats our trace like any voltage trace. For this one, we will have to extract the features in two steps, because they depend on which stimulus interval we are interested in. We will first look into the features we can extract from the stimulus interval between 1600 and 1700 ms.

traces = []
rep_name = "repetition1"
exp_type = "Inactivation"
for idx in range(len(data[exp_type][rep_name]["dt"])):
    trace = {}
    i = data[exp_type][rep_name]["current"][:,idx]
    t = np.arange(i.size) * data[exp_type][rep_name]["dt"][idx]
    # efel expects ms: s -> ms
    t = t * 1000.0
    trace["T"] = t
    trace["V"] = i  # trick: input current as if it was voltage
    trace["stim_start"] = [1600]
    trace["stim_end"] = [1700]
    traces.append(trace)

As usual, we will start by plotting the traces:

for t in traces:
    plt.plot(t["T"], t["V"], c="red")
plt.xlabel("Stimulus Voltage (mV)")
plt.ylabel("Response Current (nA)")
plt.title("Inactivation traces")

Text(0.5, 1.0, 'Inactivation traces')

We will now plot the I-V curve showing the maximum current in function of the input voltage. Remember that even though we are using the ‘maximum_voltage’ feature, since the trace we have given as input under ‘V’ is actually current, we will actually get the maximum current by using this feature:

efel.reset()
stim_v = list(range(-40, 80, 10))
feature_names = ["maximum_voltage"]
feats = efel.get_feature_values(traces, feature_names)

max_i = np.array([feat_dict["maximum_voltage"][0] for feat_dict in feats])
max_i = max_i / np.max(max_i)

plt.clf()
plt.plot(stim_v, max_i, '.')
plt.xlabel("Stimulus Voltage (mV)")
plt.ylabel("Peak current / Overall max current")
plt.title("Inactivation I-V curve")
plt.ylim((-0.05, 1.05))

(-0.05, 1.05)

We will now compute the inactivation time constant. But in order to do that, we need to change the stimulus start and stimulus end to 100 and 1600 respectively:

traces = []
rep_name = "repetition1"
exp_type = "Inactivation"
for idx in range(len(data[exp_type][rep_name]["dt"])):
    trace = {}
    i = data[exp_type][rep_name]["current"][:,idx]
    t = np.arange(i.size) * data[exp_type][rep_name]["dt"][idx]
    # efel expects ms: s -> ms
    t = t * 1000.0
    trace["T"] = t
    trace["V"] = i  # trick: input current as if it was voltage
    trace["stim_start"] = [100]
    trace["stim_end"] = [1600]
    traces.append(trace)

And now we can compute the inactivation time constant. It will select the trace from the trace maximum up to the end of the stimulus (minus a few data points to avoid artifacts).

feature_names = ["inactivation_time_constant"]
feats = efel.get_feature_values(traces, feature_names)
inact_tau_efel = np.array([feat_dict["inactivation_time_constant"][0] for feat_dict in feats])
plt.plot(stim_v, inact_tau_efel, '.')
plt.xlim((-25, 75))
plt.ylim((-5, 450))
plt.title("Inactivation time constant")
plt.xlabel("Stimulus voltage (mV)")
plt.ylabel("Time constant (ms)")

Text(0, 0.5, 'Time constant (ms)')

Note that if you want to change the number of data points to remove from the end of the trace, you can change it by changing the setting of inactivation_tc_end_skip. Setting it to 0, for example, will keep the artifacts in and can bias the computation, depending on the trace.

efel.set_setting("inactivation_tc_end_skip", 0)
feature_names = ["inactivation_time_constant", "steady_state_voltage_stimend"]
feats = efel.get_feature_values(traces, feature_names)
inact_tau_efel = np.array([feat_dict["inactivation_time_constant"][0] for feat_dict in feats])
plt.plot(stim_v, inact_tau_efel, '.')
plt.xlim((-25, 75))
plt.title("Inactivation time constant")
plt.xlabel("Stimulus voltage (mV)")
plt.ylabel("Time constant (ms)")

/Users/ajaquier/Documents/eFEL/efel/pyfeatures/pyfeatures.py:477: OptimizeWarning: Covariance of the parameters could not be estimated
  popt, _ = curve_fit(

Text(0, 0.5, 'Time constant (ms)')

As a visual verification, we can do an exponential fit to the traces, use our time constant as a parameter of the exponential and plot the exponential on top of the trace.

for i, trace in enumerate(traces):
        # - 1 to account for artifacts
        interval_indices = np.where((trace["T"] >= trace["stim_start"][0]) & (trace["T"] < trace["stim_end"][0] - 1))
        t_interval = trace["T"][interval_indices]
        v_interval = trace["V"][interval_indices]
        stim_start_idx = np.argwhere(v_interval >= np.max(v_interval))[0][0]
        t_interval = t_interval[stim_start_idx:]
        v_interval = v_interval[stim_start_idx:]
        t0 = t_interval[0] # for testing
        t_interval_corrected = t_interval - t_interval[0]

        popt, _ = curve_fit(exp_fit, t_interval_corrected, v_interval)

        v_fit = exp_fit(t_interval_corrected, inact_tau_efel[i], popt[1], popt[2])
        plt.plot(trace["T"], trace["V"], c="red")
        plt.plot(t_interval, v_fit, "--", c="black", lw=3)
        t_tau = popt[0] + t0
        plt.plot(t_tau, exp_fit(popt[0], *popt), 'o', color='black', markersize=12)
# plt.xlim(390, 610)

/var/folders/9j/1kf4nxwd62v5vnvhcnf87h616469rk/T/ipykernel_12656/1772092254.py:2: RuntimeWarning: overflow encountered in exp
  return A0 - A1 * np.exp(-t / tau_1)
/var/folders/9j/1kf4nxwd62v5vnvhcnf87h616469rk/T/ipykernel_12656/870918093.py:12: OptimizeWarning: Covariance of the parameters could not be estimated
  popt, _ = curve_fit(exp_fit, t_interval_corrected, v_interval)

Finally, the steady state of the Inactivation curve can be plotted as follow:

ssis = np.array([feat_dict["steady_state_voltage_stimend"][0] for feat_dict in feats])
plt.plot(stim_v, ssis, '.')
plt.xlabel("Stimulus Voltage (mV)")
plt.ylabel("Steady State Current (nA)")
plt.title("Inactivation steady state")

popt, _ = curve_fit(sigmoid_scaled, stim_v, ssis)
v = np.linspace(stim_v[0], stim_v[-1], 200)
ss = sigmoid_scaled(v, *popt)
plt.plot(v, ss, "--", c="black")

[<matplotlib.lines.Line2D at 0x135f9bcd0>]

Source

voltage_clamp.rst