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) .. code:: ipython3 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. .. code:: ipython3 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. .. code:: ipython3 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: .. code:: ipython3 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") .. parsed-literal:: Text(0.5, 1.0, 'Activation traces') .. image:: voltage_clamp_files/voltage_clamp_8_1.png 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. .. code:: ipython3 feature_names = ["maximum_voltage", "voltage_base", "steady_state_voltage_stimend", "activation_time_constant"] feats = efel.get_feature_values(traces, feature_names) .. parsed-literal:: /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. .. code:: ipython3 # 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") .. parsed-literal:: 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 .. code:: ipython3 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") .. parsed-literal:: [] .. image:: voltage_clamp_files/voltage_clamp_14_1.png We can also show the evolution of the steady state current at the end of the stimulus across input voltage: .. code:: ipython3 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") .. parsed-literal:: [] .. image:: voltage_clamp_files/voltage_clamp_16_1.png 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. .. code:: ipython3 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") .. parsed-literal:: Text(0.5, 1.0, 'Activation current base') .. image:: voltage_clamp_files/voltage_clamp_18_1.png 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): .. code:: ipython3 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") .. parsed-literal:: Text(0.5, 1.0, 'Activation time constant') .. image:: voltage_clamp_files/voltage_clamp_21_1.png 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. .. code:: ipython3 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) .. parsed-literal:: (90.0, 170.0) .. image:: voltage_clamp_files/voltage_clamp_23_1.png 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 .. code:: ipython3 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: .. code:: ipython3 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") .. parsed-literal:: Text(0.5, 1.0, 'Deactivation traces') .. image:: voltage_clamp_files/voltage_clamp_27_1.png We can now extract the deactivation time constant using eFEL: .. code:: ipython3 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)) .. parsed-literal:: (-10.0, 200.0) .. image:: voltage_clamp_files/voltage_clamp_29_1.png 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. .. code:: ipython3 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) .. parsed-literal:: /var/folders/9j/1kf4nxwd62v5vnvhcnf87h616469rk/T/ipykernel_12656/1772092254.py:2: RuntimeWarning: overflow encountered in exp return A0 - A1 * np.exp(-t / tau_1) .. parsed-literal:: (390.0, 610.0) .. image:: voltage_clamp_files/voltage_clamp_31_2.png 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. .. code:: ipython3 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: .. code:: ipython3 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") .. parsed-literal:: Text(0.5, 1.0, 'Inactivation traces') .. image:: voltage_clamp_files/voltage_clamp_35_1.png 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: .. code:: ipython3 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)) .. parsed-literal:: (-0.05, 1.05) .. image:: voltage_clamp_files/voltage_clamp_37_1.png 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: .. code:: ipython3 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). .. code:: ipython3 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)") .. parsed-literal:: Text(0, 0.5, 'Time constant (ms)') .. image:: voltage_clamp_files/voltage_clamp_41_1.png 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. .. code:: ipython3 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)") .. parsed-literal:: /Users/ajaquier/Documents/eFEL/efel/pyfeatures/pyfeatures.py:477: OptimizeWarning: Covariance of the parameters could not be estimated popt, _ = curve_fit( .. parsed-literal:: Text(0, 0.5, 'Time constant (ms)') .. image:: voltage_clamp_files/voltage_clamp_43_2.png 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. .. code:: ipython3 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) .. parsed-literal:: /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) .. image:: voltage_clamp_files/voltage_clamp_45_1.png Finally, the steady state of the Inactivation curve can be plotted as follow: .. code:: ipython3 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") .. parsed-literal:: [] .. image:: voltage_clamp_files/voltage_clamp_47_1.png