Add BPAP validation

bluecellulab/analysis/analysis.py

@@ -3,14 +3,26 @@
     import efel
 except ImportError:
     efel = None
+from itertools import islice
+import logging
+import matplotlib.pyplot as plt
+import neuron
 import numpy as np
+import pathlib
 
+from bluecellulab import Cell
 from bluecellulab.analysis.inject_sequence import run_stimulus
 from bluecellulab.analysis.plotting import plot_iv_curve, plot_fi_curve
+from bluecellulab.analysis.utils import exp_decay
+from bluecellulab.simulation import Simulation
 from bluecellulab.stimulus import StimulusFactory
+from bluecellulab.stimulus.circuit_stimulus_definitions import Hyperpolarizing
 from bluecellulab.tools import calculate_rheobase
 
 
+logger = logging.getLogger(__name__)
+
+
 def compute_plot_iv_curve(cell,
                           injecting_section="soma[0]",
                           injecting_segment=0.5,
@@ -48,7 +60,7 @@
         post_delay (float, optional): The delay after the stimulation ends before the simulation stops
             (in ms). Default is 100.0 ms.
         threshold_voltage (float, optional): The voltage threshold (in mV) for detecting a steady-state
-            response. Default is -30 mV.
+            response. Default is -20 mV.
         nb_bins (int, optional): The number of discrete current levels between 0 and the maximum current.
             Default is 11.
         rheobase (float, optional): The rheobase current (in nA) for the cell. If not provided, it will
@@ -158,7 +170,7 @@
         max_current (float, optional): The maximum amplitude of the injected current (in nA).
             Default is 0.8 nA.
         threshold_voltage (float, optional): The voltage threshold (in mV) for detecting a steady-state
-            response. Default is -30 mV.
+            response. Default is -20 mV.
         nb_bins (int, optional): The number of discrete current levels between 0 and `max_current`.
             Default is 11.
         rheobase (float, optional): The rheobase current (in nA) for the cell. If not provided, it will
@@ -211,3 +223,170 @@
                   output_fname=output_fname)
 
     return np.array(list_amp), np.array(spike_count)
+
+
+class BPAP:
+    # taken from the examples
+
+    def __init__(self, cell: Cell) -> None:
+        self.cell = cell
+        self.dt = 0.025
+        self.stim_start = 1000
+        self.stim_duration = 1
+
+    @property
+    def start_index(self) -> int:
+        """Get the index of the start of the stimulus."""
+        return int(self.stim_start / self.dt)
+
+    @property
+    def end_index(self) -> int:
+        """Get the index of the end of the stimulus."""
+        return int((self.stim_start + self.stim_duration) / self.dt)
+
+    def get_recordings(self):
+        """Get the soma, basal and apical recordings."""
+        all_recordings = self.cell.get_allsections_voltagerecordings()
+        soma_rec = None
+        dend_rec = {}
+        apic_rec = {}
+        for key, value in all_recordings.items():
+            if "soma" in key:
+                soma_rec = value
+            elif "dend" in key:
+                dend_rec[key] = value
+            elif "apic" in key:
+                apic_rec[key] = value
+
+        return soma_rec, dend_rec, apic_rec
+
+    def run(self, duration: float, amplitude: float) -> None:
+        """Apply depolarization and hyperpolarization at the same time."""
+        sim = Simulation()
+        sim.add_cell(self.cell)
+        self.cell.add_allsections_voltagerecordings()
+        self.cell.add_step(start_time=self.stim_start, stop_time=self.stim_start + self.stim_duration, level=amplitude)
+        hyperpolarizing = Hyperpolarizing("single-cell", delay=0, duration=duration)
+        self.cell.add_replay_hypamp(hyperpolarizing)
+        sim.run(duration, dt=self.dt, cvode=False)
+
+    def amplitudes(self, recs) -> list[float]:
+        """Return amplitude across given sections."""
+        efel_feature_name = "maximum_voltage_from_voltagebase"
+        traces = [
+            {
+                'T': self.cell.get_time(),
+                'V': rec,
+                'stim_start': [self.stim_start],
+                'stim_end': [self.stim_start + self.stim_duration]
+            }
+            for rec in recs.values()
+        ]
+        features_results = efel.get_feature_values(traces, [efel_feature_name])
+        amps = [feat_res[efel_feature_name][0] for feat_res in features_results]
+
+        return amps
+
+    def distances_to_soma(self, recs) -> list[float]:
+        """Return the distance to the soma for each section."""
+        res = []
+        soma = self.cell.soma
+        for key in recs.keys():
+            section_name = key.rsplit(".")[-1].split("[")[0]  # e.g. "dend"
+            section_idx = int(key.rsplit(".")[-1].split("[")[1].split("]")[0])  # e.g. 0
+            attribute_value = getattr(self.cell.cell.getCell(), section_name)
+            section = next(islice(attribute_value, section_idx, None))
+            # section e.g. cADpyr_L2TPC_bluecellulab_x[0].dend[0]
+            res.append(neuron.h.distance(soma(0.5), section(0.5)))
+        return res
+
+    def get_amplitudes_and_distances(self):
+        soma_rec, dend_rec, apic_rec = self.get_recordings()
+        soma_amp = self.amplitudes({"soma": soma_rec})[0]
+        dend_amps = None
+        dend_dist = None
+        apic_amps = None
+        apic_dist = None
+        if dend_rec:
+            dend_amps = self.amplitudes(dend_rec)
+            dend_dist = self.distances_to_soma(dend_rec)
+        if apic_rec:
+            apic_amps = self.amplitudes(apic_rec)
+            apic_dist = self.distances_to_soma(apic_rec)
+
+        return soma_amp, dend_amps, dend_dist, apic_amps, apic_dist
+
+    def fit(self, soma_amp, dend_amps, dend_dist, apic_amps, apic_dist):
+        """Fit the amplitudes vs distances to an exponential decay function."""
+        from scipy.optimize import curve_fit
+
+        popt_dend = None
+        if dend_amps and dend_dist:
+            dist = [0] + dend_dist  # add soma distance
+            amps = [soma_amp] + dend_amps  # add soma amplitude
+            popt_dend, _ = curve_fit(exp_decay, dist, amps)
+
+        popt_apic = None
+        if apic_amps and apic_dist:
+            dist = [0] + apic_dist  # add soma distance
+            amps = [soma_amp] + apic_amps  # add soma amplitude
+            popt_apic, _ = curve_fit(exp_decay, dist, amps)
+
+        return popt_dend, popt_apic
+
+    def validate(self, soma_amp, dend_amps, dend_dist, apic_amps, apic_dist):
+        """Check that the exponential fit is decaying."""
+        validated = True
+        notes = ""
+        popt_dend, popt_apic = self.fit(soma_amp, dend_amps, dend_dist, apic_amps, apic_dist)
+        if popt_dend is None:
+            logger.debug("No dendritic recordings found.")
+            notes += "No dendritic recordings found.\n"
+        elif popt_dend[1] <= 0:
+            logger.debug("Dendritic fit is not decaying.")
+            validated = False
+            notes += "Dendritic fit is not decaying.\n"
+        else:
+            notes += "Dendritic validation passed: dendritic amplitude is decaying with distance relative to soma.\n"
+        if popt_apic is None:
+            logger.debug("No apical recordings found.")
+            notes += "No apical recordings found.\n"
+        elif popt_apic[1] <= 0:
+            logger.debug("Apical fit is not decaying.")
+            validated = False
+            notes += "Apical fit is not decaying.\n"
+        else:
+            notes += "Apical validation passed: apical amplitude is decaying with distance relative to soma.\n"
+
+        return validated, notes
+
+    def plot(self, soma_amp, dend_amps, dend_dist, apic_amps, apic_dist, show_figure=True, save_figure=False, output_dir="./", output_fname="bpap.pdf"):
+        """Plot the results of the BPAP analysis."""
+        popt_dend, popt_apic = self.fit(soma_amp, dend_amps, dend_dist, apic_amps, apic_dist)
+
+        outpath = pathlib.Path(output_dir) / output_fname
+        fig, ax1 = plt.subplots(figsize=(10, 6))
+        ax1.scatter([0], [soma_amp], color='black', label='Soma')
+        if dend_amps and dend_dist:
+            ax1.scatter(dend_dist, dend_amps, color='green', label='Dendrites')
+            if popt_dend is not None:
+                x = np.linspace(0, max(dend_dist), 100)
+                y = exp_decay(x, *popt_dend)
+                ax1.plot(x, y, color='darkgreen', linestyle='--', label='Dendritic Fit')
+        if apic_amps and apic_dist:
+            ax1.scatter(apic_dist, apic_amps, color='blue', label='Apical Dendrites')
+            if popt_apic is not None:
+                x = np.linspace(0, max(apic_dist), 100)
+                y = exp_decay(x, *popt_apic)
+                ax1.plot(x, y, color='darkblue', linestyle='--', label='Apical Fit')
+        ax1.set_xlabel('Distance to Soma (um)')
+        ax1.set_ylabel('Amplitude (mV)')
+        ax1.legend()
+        fig.suptitle('Back-propagating Action Potential Analysis')
+        fig.tight_layout()
+        if save_figure:
+            fig.savefig(outpath)
+        if show_figure:
+            plt.show()
+
+        return outpath

bluecellulab/analysis/utils.py

@@ -0,0 +1,7 @@
+"""Utility functions for analysis."""
+
+import numpy as np
+
+
+def exp_decay(x, a, b, c):
+    return a * np.exp(-b * x) + c