Add FI curve (#5)

* add FI curve

* support injecting and recording location in run_stimulus

* add injecting and recording location support for iv and fi curves

* add a helper function to validate a section/segment position

* enhance spikedetector with support for custom locations

* add unit tests for validate_section_and_segment

* add unit tests for run_stimulus

* improve robustness of 'create_netcon_spikedetector'

* add unit tests

* add support for calculating rheobase at specific sections/segments

* improve error handling in  method

* add unit tests

Author: ilkilic

bluecellulab/analysis/analysis.py

@@ -8,37 +8,72 @@
 from bluecellulab.stimulus import StimulusFactory
 from bluecellulab.tools import calculate_rheobase
 from bluecellulab.analysis.inject_sequence import run_stimulus
-from bluecellulab.analysis.plotting import plot_iv_curve
+from bluecellulab.analysis.plotting import plot_iv_curve, plot_fi_curve
 
 
-def compute_plot_iv_curve(cell, stim_start=100.0, duration=500.0, post_delay=100.0, threshold_voltage=-30, nb_bins=11):
-    """Compute and plot the IV curve from a given cell by injecting a
-    predefined range of currents.
+def compute_plot_iv_curve(cell,
+                          injecting_section="soma[0]",
+                          injecting_segment=0.5,
+                          recording_section="soma[0]",
+                          recording_segment=0.5,
+                          stim_start=100.0,
+                          duration=500.0,
+                          post_delay=100.0,
+                          threshold_voltage=-30,
+                          nb_bins=11):
+    """Compute and plot the Current-Voltage (I-V) curve for a given cell by
+    injecting a range of currents.
+
+    This function evaluates the relationship between the injected current amplitude and the resulting
+    steady-state membrane potential of a neuronal cell model. Currents are injected at a specified section
+    and segment, and the steady-state voltage at the recording location is used to construct the I-V curve.
 
     Args:
-        cell (bluecellulab.cell.Cell): The initialized cell model.
-        stim_start (float): Start time for current injection (in ms). Default is 100.0 ms.
-        duration (float): Duration of current injection (in ms). Default is 500.0 ms.
-        post_delay (float): Delay after the stimulation ends (in ms). Default is 100.0 ms.
-        nb_bins (int): Number of current injection levels. Default is 11.
+        cell (bluecellulab.cell.Cell): The initialized BlueCelluLab cell model.
+        injecting_section (str, optional): The name of the section where the stimulus is injected.
+            Default is "soma[0]".
+        injecting_segment (float, optional): The position along the injecting section (0.0 to 1.0)
+            where the stimulus is applied. Default is 0.5.
+        recording_section (str, optional): The name of the section where the voltage is recorded.
+            Default is "soma[0]".
+        recording_segment (float, optional): The position along the recording section (0.0 to 1.0)
+            where the voltage is recorded. Default is 0.5.
+        stim_start (float, optional): The start time of the current injection (in ms). Default is 100.0 ms.
+        duration (float, optional): The duration of the current injection (in ms). Default is 500.0 ms.
+        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.
+        nb_bins (int, optional): The number of discrete current levels between 0 and the maximum current.
+            Default is 11.
 
     Returns:
         tuple: A tuple containing:
-            - list_amp (np.ndarray): The predefined injected step current levels (nA).
-            - steady_states (np.ndarray): The corresponding steady-state voltages (mV).
+            - list_amp (np.ndarray): The injected current amplitudes (nA).
+            - steady_states (np.ndarray): The corresponding steady-state voltages (mV) recorded at the
+              specified location.
+
+    Raises:
+        ValueError: If the cell object is invalid, the specified sections/segments are not found, or if
+            the simulation results are inconsistent.
     """
-    rheobase = calculate_rheobase(cell)
+    rheobase = calculate_rheobase(cell=cell, section=injecting_section, segx=injecting_segment)
 
     list_amp = np.linspace(rheobase - 2, rheobase - 0.1, nb_bins)  # [nA]
 
     steps = []
     times = []
     voltages = []
     # inject step current and record voltage response
+    stim_factory = StimulusFactory(dt=0.1)
     for amp in list_amp:
-        stim_factory = StimulusFactory(dt=0.1)
         step_stimulus = stim_factory.step(pre_delay=stim_start, duration=duration, post_delay=post_delay, amplitude=amp)
-        recording = run_stimulus(cell.template_params, step_stimulus, section="soma[0]", segment=0.5)
+        recording = run_stimulus(cell.template_params,
+                                 step_stimulus,
+                                 section=injecting_section,
+                                 segment=injecting_segment,
+                                 recording_section=recording_section,
+                                 recording_segment=recording_segment)
         steps.append(step_stimulus)
         times.append(recording.time)
         voltages.append(recording.voltage)
@@ -57,7 +92,86 @@ def compute_plot_iv_curve(cell, stim_start=100.0, duration=500.0, post_delay=100
         steady_state = features_results[0]['steady_state_voltage_stimend']
         steady_states.append(steady_state)
 
-    # plot I-V curve
-    plot_iv_curve(list_amp, steady_states)
+    plot_iv_curve(list_amp,
+                  steady_states,
+                  injecting_section=injecting_section,
+                  injecting_segment=injecting_segment,
+                  recording_section=recording_section,
+                  recording_segment=recording_segment)
 
     return np.array(list_amp), np.array(steady_states)
+
+
+def compute_plot_fi_curve(cell,
+                          injecting_section="soma[0]",
+                          injecting_segment=0.5,
+                          recording_section="soma[0]",
+                          recording_segment=0.5,
+                          stim_start=100.0,
+                          duration=500.0,
+                          post_delay=100.0,
+                          max_current=0.8,
+                          nb_bins=11):
+    """Compute and plot the Frequency-Current (F-I) curve for a given cell by
+    injecting a range of currents.
+
+    This function evaluates the relationship between injected current amplitude and the firing rate
+    of a neuronal cell model. Currents are injected at a specified section and segment, and the number
+    of spikes recorded in the specified recording location is used to construct the F-I curve.
+
+    Args:
+        cell (bluecellulab.cell.Cell): The initialized BlueCelluLab cell model.
+        injecting_section (str, optional): The name of the section where the stimulus is injected.
+            Default is "soma[0]".
+        injecting_segment (float, optional): The position along the injecting section (0.0 to 1.0)
+            where the stimulus is applied. Default is 0.5.
+        recording_section (str, optional): The name of the section where spikes are recorded.
+            Default is "soma[0]".
+        recording_segment (float, optional): The position along the recording section (0.0 to 1.0)
+            where spikes are recorded. Default is 0.5.
+        stim_start (float, optional): The start time of the current injection (in ms). Default is 100.0 ms.
+        duration (float, optional): The duration of the current injection (in ms). Default is 500.0 ms.
+        post_delay (float, optional): The delay after the stimulation ends before the simulation stops
+            (in ms). Default is 100.0 ms.
+        max_current (float, optional): The maximum amplitude of the injected current (in nA).
+            Default is 0.8 nA.
+        nb_bins (int, optional): The number of discrete current levels between 0 and `max_current`.
+            Default is 11.
+
+    Returns:
+        tuple: A tuple containing:
+            - list_amp (np.ndarray): The injected current amplitudes (nA).
+            - spike_count (np.ndarray): The corresponding spike counts for each current amplitude.
+
+    Raises:
+        ValueError: If the cell object is invalid or the specified sections/segments are not found.
+    """
+    rheobase = calculate_rheobase(cell=cell, section=injecting_section, segx=injecting_segment)
+
+    list_amp = np.linspace(rheobase, max_current, nb_bins)  # [nA]
+    steps = []
+    spikes = []
+    # inject step current and record spike response
+    stim_factory = StimulusFactory(dt=0.1)
+    for amp in list_amp:
+        step_stimulus = stim_factory.step(pre_delay=stim_start, duration=duration, post_delay=post_delay, amplitude=amp)
+        recording = run_stimulus(cell.template_params,
+                                 step_stimulus,
+                                 section=injecting_section,
+                                 segment=injecting_segment,
+                                 recording_section=recording_section,
+                                 recording_segment=recording_segment,
+                                 enable_spike_detection=True)
+        steps.append(step_stimulus)
+        spikes.append(recording.spike)
+
+    spike_count = [len(spike) for spike in spikes]
+
+    plot_fi_curve(list_amp,
+                  spike_count,
+                  injecting_section=injecting_section,
+                  injecting_segment=injecting_segment,
+                  recording_section=recording_section,
+                  recording_segment=recording_segment)
+
+    return np.array(list_amp), np.array(spike_count)

bluecellulab/analysis/inject_sequence.py

@@ -1,7 +1,7 @@
 """Module for injecting a sequence of protocols to the cell."""
 from __future__ import annotations
 from enum import Enum, auto
-from typing import NamedTuple, Sequence, Dict
+from typing import NamedTuple, Sequence, Dict, Optional
 
 import neuron
 import numpy as np
@@ -11,6 +11,7 @@
 from bluecellulab.simulation.simulation import Simulation
 from bluecellulab.stimulus.circuit_stimulus_definitions import Hyperpolarizing
 from bluecellulab.stimulus.factory import Stimulus, StimulusFactory
+from bluecellulab.tools import validate_section_and_segment
 
 
 class StimulusName(Enum):
@@ -24,10 +25,12 @@ class StimulusName(Enum):
 
 
 class Recording(NamedTuple):
-    """A tuple of the current, voltage and time recordings."""
+    """A tuple of the current, voltage and time recordings with optional spike
+    recordings."""
     current: np.ndarray
     voltage: np.ndarray
     time: np.ndarray
+    spike: np.ndarray | None
 
 
 StimulusRecordings = Dict[str, Recording]
@@ -40,41 +43,91 @@ def run_stimulus(
     segment: float,
     cvode: bool = True,
     add_hypamp: bool = True,
+    recording_section: str = "soma[0]",
+    recording_segment: float = 0.5,
+    enable_spike_detection: bool = False,
+    threshold_spike_detection: float = -30,
 ) -> Recording:
-    """Creates a cell and stimulates it with a given stimulus.
+    """Creates a cell from template parameters, applies a stimulus, and records
+    the response.
+
+    This function simulates the electrical activity of a neuronal cell model by injecting
+    a stimulus into a specified section and segment, recording the voltage response, and
+    optionally detecting spikes.
 
     Args:
-        template_params: The parameters to create the cell from a template.
-        stimulus: The input stimulus to inject into the cell.
-        section: Name of the section of cell where the stimulus is to be injected.
-        segment: The segment of the section where the stimulus is to be injected.
-        cvode: True to use variable time-steps. False for fixed time-steps.
+        template_params (TemplateParams): Parameters required to create the cell from a
+            specified template, including morphology and mechanisms.
+        stimulus (Stimulus): The stimulus waveform to inject, defined by time and current arrays.
+        section (str): Name of the section of the cell where the stimulus is applied.
+            (e.g. soma[0])
+        segment (float): The normalized position (0.0 to 1.0) along the injecting
+            section where the stimulus is applied.
+        cvode (bool, optional): Whether to use variable time-step integration. Defaults to True.
+        add_hypamp (bool, optional): If True, adds a hyperpolarizing stimulus before applying
+            the main stimulus. Defaults to True.
+        recording_section (str): Name of the section of the cell where voltage is recorded.
+        recording_segment (float): The normalized position (0.0 to 1.0) along the recording
+            section where voltage is recorded.
+        enable_spike_detection (bool, optional): If True, enables spike detection at the
+            recording location. Defaults to False.
+        threshold_spike_detection (float, optional): The voltage threshold (mV) for spike detection.
+            Defaults to -30 mV.
 
     Returns:
-        The voltage-time recording at the specified location.
+        Recording: A `Recording` object containing the following:
+            - `current` (np.ndarray): The injected current waveform (nA).
+            - `voltage` (np.ndarray): The recorded membrane potential (mV) over time.
+            - `time` (np.ndarray): The simulation time points (ms).
+            - `spike` (np.ndarray or None): The detected spikes, if spike detection is enabled.
 
     Raises:
-        ValueError: If the time and voltage arrays are not the same length.
+        ValueError: If the time, current, and voltage arrays do not have the same length,
+            or if the specified sections or segments are not found in the cell model.
     """
     cell = Cell.from_template_parameters(template_params)
-    neuron_section = cell.sections[section]
+
+    validate_section_and_segment(cell, section, segment)
+    validate_section_and_segment(cell, recording_section, recording_segment)
+
     if add_hypamp:
         hyp_stim = Hyperpolarizing(target="", delay=0.0, duration=stimulus.stimulus_time)
         cell.add_replay_hypamp(hyp_stim)
-    cell.add_voltage_recording(neuron_section, segment)
+
+    cell.add_voltage_recording(cell.sections[recording_section], recording_segment)
+
+    # Set up spike detection if enabled
+    spikes: Optional[np.ndarray] = None
+    if enable_spike_detection:
+        recording_location = f"{recording_section}({str(recording_segment)})"
+        cell.start_recording_spikes(None, location=recording_location, threshold=threshold_spike_detection)
+
+    # Inject the stimulus and run the simulation
     iclamp, _ = cell.inject_current_waveform(
-        stimulus.time, stimulus.current, section=neuron_section, segx=segment
+        stimulus.time, stimulus.current, section=cell.sections[section], segx=segment
     )
     current_vector = neuron.h.Vector()
     current_vector.record(iclamp._ref_i)
+
     simulation = Simulation(cell)
     simulation.run(stimulus.stimulus_time, cvode=cvode)
+
+    # Retrieve simulation results
     current = np.array(current_vector.to_python())
-    voltage = cell.get_voltage_recording(neuron_section, segment)
+    voltage = cell.get_voltage_recording(cell.sections[recording_section], recording_segment)
     time = cell.get_time()
+
     if len(time) != len(voltage) or len(time) != len(current):
-        raise ValueError("Time, current and voltage arrays are not the same length")
-    return Recording(current, voltage, time)
+        raise ValueError("Time, current, and voltage arrays are not the same length")
+
+    if enable_spike_detection:
+        results = cell.get_recorded_spikes(location=recording_location, threshold=threshold_spike_detection)
+        if results is not None:
+            spikes = np.array(results)
+        else:
+            spikes = None
+
+    return Recording(current=current, voltage=voltage, time=time, spike=spikes)
 
 
 def apply_multiple_stimuli(

bluecellulab/analysis/plotting.py

@@ -3,23 +3,55 @@
 import matplotlib.pyplot as plt
 
 
-def plot_iv_curve(currents, voltages):
+def plot_iv_curve(currents, voltages, injecting_section, injecting_segment, recording_section, recording_segment):
     """Plots the IV curve.
 
     Args:
-        currents (iterable): The injected current levels (nA).
-        voltages (iterable): The corresponding steady-state voltages (mV).
+        currents (list): The injected current levels (nA).
+        voltages (list): The corresponding steady-state voltages (mV).
+        injecting_section (str): The section in the cell where the current was injected.
+        injecting_segment (float): The segment position (0.0 to 1.0) where the current was injected.
+        recording_section (str): The section in the cell where spikes were recorded.
+        recording_segment (float): The segment position (0.0 to 1.0) where spikes were recorded.
+
     Raises:
         ValueError: If the lengths of currents and voltages do not match.
     """
     if len(currents) != len(voltages):
         raise ValueError("currents and voltages must have the same length")
 
-    # Plotting
     plt.figure(figsize=(10, 6))
-    plt.plot(voltages, currents, marker='o', linestyle='-', color='b')
+    plt.plot(currents, voltages, marker='o', linestyle='-', color='b')
     plt.title("I-V Curve")
-    plt.ylabel("Injected current [nA]")
-    plt.xlabel("Steady state voltage [mV]")
+    plt.xlabel(f"Injected Current [nA] at {injecting_section}({injecting_segment:.2f})")
+    plt.ylabel(f"Steady state voltage [mV] at {recording_section}({recording_segment:.2f})")
+    plt.grid(True)
+    plt.tight_layout()
+    plt.show()
+
+
+def plot_fi_curve(currents, spike_count, injecting_section, injecting_segment, recording_section, recording_segment):
+    """Plots the F-I (Frequency-Current) curve.
+
+    Args:
+        currents (list): The injected current levels (nA).
+        spike_count (list): The number of spikes recorded for each current level.
+        injecting_section (str): The section in the cell where the current was injected.
+        injecting_segment (float): The segment position (0.0 to 1.0) where the current was injected.
+        recording_section (str): The section in the cell where spikes were recorded.
+        recording_segment (float): The segment position (0.0 to 1.0) where spikes were recorded.
+
+    Raises:
+        ValueError: If the lengths of currents and spike counts do not match.
+    """
+    if len(currents) != len(spike_count):
+        raise ValueError("currents and spike count must have the same length")
+
+    plt.figure(figsize=(10, 6))
+    plt.plot(currents, spike_count, marker='o')
+    plt.title("F-I Curve")
+    plt.xlabel(f"Injected Current [nA] at {injecting_section}({injecting_segment:.2f})")
+    plt.ylabel(f"Spike Count recorded at {recording_section}({recording_segment:.2f})")
+    plt.grid(True)
     plt.tight_layout()
     plt.show()