Memodel validation

bluecellulab/analysis/analysis.py

@@ -5,10 +5,10 @@
     efel = None
 import numpy as np
 
-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, plot_fi_curve
+from bluecellulab.stimulus import StimulusFactory
+from bluecellulab.tools import calculate_rheobase
 
 
 def compute_plot_iv_curve(cell,
@@ -19,8 +19,13 @@ def compute_plot_iv_curve(cell,
                           stim_start=100.0,
                           duration=500.0,
                           post_delay=100.0,
-                          threshold_voltage=-30,
-                          nb_bins=11):
+                          threshold_voltage=-20,
+                          nb_bins=11,
+                          rheobase=None,
+                          show_figure=True,
+                          save_figure=False,
+                          output_dir="./",
+                          output_fname="iv_curve.pdf"):
     """Compute and plot the Current-Voltage (I-V) curve for a given cell by
     injecting a range of currents.
 
@@ -46,6 +51,12 @@ def compute_plot_iv_curve(cell,
             response. Default is -30 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
+            be calculated using the `calculate_rheobase` function.
+        show_figure (bool): Whether to display the figure. Default is True.
+        save_figure (bool): Whether to save the figure. Default is False.
+        output_dir (str): The directory to save the figure if save_figure is True. Default is "./".
+        output_fname (str): The filename to save the figure as if save_figure is True. Default is "iv_curve.png".
 
     Returns:
         tuple: A tuple containing:
@@ -57,7 +68,8 @@ def compute_plot_iv_curve(cell,
         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=cell, section=injecting_section, segx=injecting_segment)
+    if rheobase is None:
+        rheobase = calculate_rheobase(cell=cell, section=injecting_section, segx=injecting_segment)
 
     list_amp = np.linspace(rheobase - 2, rheobase - 0.1, nb_bins)  # [nA]
 
@@ -89,15 +101,19 @@ def compute_plot_iv_curve(cell,
             'stim_end': [stim_start + duration]
         }
         features_results = efel.get_feature_values([trace], ['steady_state_voltage_stimend'])
-        steady_state = features_results[0]['steady_state_voltage_stimend']
+        steady_state = features_results[0]['steady_state_voltage_stimend'][0]
         steady_states.append(steady_state)
 
     plot_iv_curve(list_amp,
                   steady_states,
                   injecting_section=injecting_section,
                   injecting_segment=injecting_segment,
                   recording_section=recording_section,
-                  recording_segment=recording_segment)
+                  recording_segment=recording_segment,
+                  show_figure=show_figure,
+                  save_figure=save_figure,
+                  output_dir=output_dir,
+                  output_fname=output_fname)
 
     return np.array(list_amp), np.array(steady_states)
 
@@ -111,7 +127,13 @@ def compute_plot_fi_curve(cell,
                           duration=500.0,
                           post_delay=100.0,
                           max_current=0.8,
-                          nb_bins=11):
+                          threshold_voltage=-20,
+                          nb_bins=11,
+                          rheobase=None,
+                          show_figure=True,
+                          save_figure=False,
+                          output_dir="./",
+                          output_fname="fi_curve.pdf"):
     """Compute and plot the Frequency-Current (F-I) curve for a given cell by
     injecting a range of currents.
 
@@ -135,8 +157,16 @@ def compute_plot_fi_curve(cell,
             (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.
+        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 `max_current`.
             Default is 11.
+        rheobase (float, optional): The rheobase current (in nA) for the cell. If not provided, it will
+            be calculated using the `calculate_rheobase` function.
+        show_figure (bool): Whether to display the figure. Default is True.
+        save_figure (bool): Whether to save the figure. Default is False.
+        output_dir (str): The directory to save the figure if save_figure is True. Default is "./".
+        output_fname (str): The filename to save the figure as if save_figure is True. Default is "iv_curve.png".
 
     Returns:
         tuple: A tuple containing:
@@ -146,7 +176,8 @@ def compute_plot_fi_curve(cell,
     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)
+    if rheobase is None:
+        rheobase = calculate_rheobase(cell=cell, section=injecting_section, segx=injecting_segment)
 
     list_amp = np.linspace(rheobase, max_current, nb_bins)  # [nA]
     steps = []
@@ -161,7 +192,8 @@ def compute_plot_fi_curve(cell,
                                  segment=injecting_segment,
                                  recording_section=recording_section,
                                  recording_segment=recording_segment,
-                                 enable_spike_detection=True)
+                                 enable_spike_detection=True,
+                                 threshold_spike_detection=threshold_voltage)
         steps.append(step_stimulus)
         spikes.append(recording.spike)
 
@@ -172,6 +204,10 @@ def compute_plot_fi_curve(cell,
                   injecting_section=injecting_section,
                   injecting_segment=injecting_segment,
                   recording_section=recording_section,
-                  recording_segment=recording_segment)
+                  recording_segment=recording_segment,
+                  show_figure=show_figure,
+                  save_figure=save_figure,
+                  output_dir=output_dir,
+                  output_fname=output_fname)
 
     return np.array(list_amp), np.array(spike_count)

bluecellulab/analysis/inject_sequence.py

@@ -36,18 +36,17 @@
 StimulusRecordings = Dict[str, Recording]
 
 
-def run_stimulus(
+def run_multirecordings_stimulus(
     template_params: TemplateParams,
     stimulus: Stimulus,
     section: str,
     segment: float,
     cvode: bool = True,
     add_hypamp: bool = True,
-    recording_section: str = "soma[0]",
-    recording_segment: float = 0.5,
+    recording_locations: list[tuple[str, float]] = [("soma[0]", 0.5)],
     enable_spike_detection: bool = False,
-    threshold_spike_detection: float = -30.0,
-) -> Recording:
+    threshold_spike_detection: float = -20.0,
+) -> list[Recording]:
     """Creates a cell from template parameters, applies a stimulus, and records
     the response.
 
@@ -66,16 +65,17 @@
         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
+        recording_location (list): List of tuples containing the name of the section of the cell
+            where voltage is recorded and the normalized position (0.0 to 1.0) along the recording
             section where voltage is recorded.
+            (e.g. [("soma[0]", 0.5), ("dend[0]", 0.5)])
         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.
+            Defaults to -20 mV.
 
     Returns:
-        Recording: A `Recording` object containing the following:
+        list[Recording]: `Recording` objects 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).
@@ -88,19 +88,22 @@
     cell = Cell.from_template_parameters(template_params)
 
     validate_section_and_segment(cell, section, segment)
-    validate_section_and_segment(cell, recording_section, recording_segment)
+    for recording_section, recording_segment in recording_locations:
+        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(cell.sections[recording_section], recording_segment)
+    for recording_section, recording_segment in recording_locations:
+        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)
+        for recording_section, recording_segment in recording_locations:
+            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(
@@ -113,21 +116,90 @@
     simulation.run(stimulus.stimulus_time, cvode=cvode)
 
     # Retrieve simulation results
+    recordings = []
     current = np.array(current_vector.to_python())
-    voltage = cell.get_voltage_recording(cell.sections[recording_section], recording_segment)
-    time = cell.get_time()
+    for recording_section, recording_segment in recording_locations:
+        recording_location = f"{recording_section}({str(recording_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")
+        if len(time) != len(voltage) or len(time) != len(current):
+            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
+        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
+
+        recordings.append(
+            Recording(current=current, voltage=voltage, time=time, spike=spikes)
+        )
 
-    return Recording(current=current, voltage=voltage, time=time, spike=spikes)
+    return recordings
+
+
+def run_stimulus(
+    template_params: TemplateParams,
+    stimulus: Stimulus,
+    section: str,
+    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 = -20.0,
+) -> Recording:
+    """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 (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 -20 mV.
+
+    Returns:
+        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, current, and voltage arrays do not have the same length,
+            or if the specified sections or segments are not found in the cell model.
+    """
+    return run_multirecordings_stimulus(
+        template_params=template_params,
+        stimulus=stimulus,
+        section=section,
+        segment=segment,
+        cvode=cvode,
+        add_hypamp=add_hypamp,
+        recording_locations=[(recording_section, recording_segment)],
+        enable_spike_detection=enable_spike_detection,
+        threshold_spike_detection=threshold_spike_detection,
+    )[0]
 
 
 def apply_multiple_stimuli(

bluecellulab/analysis/plotting.py

@@ -1,9 +1,21 @@
 """Module for plotting analysis results of cell simulations."""
 
 import matplotlib.pyplot as plt
+import pathlib
 
 
-def plot_iv_curve(currents, voltages, injecting_section, injecting_segment, recording_section, recording_segment):
+def plot_iv_curve(
+    currents,
+    voltages,
+    injecting_section,
+    injecting_segment,
+    recording_section,
+    recording_segment,
+    show_figure=True,
+    save_figure=False,
+    output_dir="./",
+    output_fname="iv_curve.pdf",
+):
     """Plots the IV curve.
 
     Args:
@@ -13,6 +25,10 @@
         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.
+        show_figure (bool): Whether to display the figure. Default is True.
+        save_figure (bool): Whether to save the figure. Default is False.
+        output_dir (str): The directory to save the figure if save_figure is True. Default is "./".
+        output_fname (str): The filename to save the figure as if save_figure is True. Default is "iv_curve.pdf".
 
     Raises:
         ValueError: If the lengths of currents and voltages do not match.
@@ -27,10 +43,27 @@
     plt.ylabel(f"Steady state voltage [mV] at {recording_section}({recording_segment:.2f})")
     plt.grid(True)
     plt.tight_layout()
-    plt.show()
+    if show_figure:
+        plt.show()
 
+    if save_figure:
+        pathlib.Path(output_dir).mkdir(parents=True, exist_ok=True)
+        plt.savefig(pathlib.Path(output_dir) / output_fname, format='pdf')
+        plt.close()
 
-def plot_fi_curve(currents, spike_count, injecting_section, injecting_segment, recording_section, recording_segment):
+
+def plot_fi_curve(
+    currents,
+    spike_count,
+    injecting_section,
+    injecting_segment,
+    recording_section,
+    recording_segment,
+    show_figure=True,
+    save_figure=False,
+    output_dir="./",
+    output_fname="fi_curve.pdf",
+):
     """Plots the F-I (Frequency-Current) curve.
 
     Args:
@@ -40,6 +73,10 @@
         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.
+        show_figure (bool): Whether to display the figure. Default is True.
+        save_figure (bool): Whether to save the figure. Default is False.
+        output_dir (str): The directory to save the figure if save_figure is True. Default is "./".
+        output_fname (str): The filename to save the figure as if save_figure is True. Default is "fi_curve.pdf".
 
     Raises:
         ValueError: If the lengths of currents and spike counts do not match.
@@ -54,4 +91,10 @@
     plt.ylabel(f"Spike Count recorded at {recording_section}({recording_segment:.2f})")
     plt.grid(True)
     plt.tight_layout()
-    plt.show()
+    if show_figure:
+        plt.show()
+
+    if save_figure:
+        pathlib.Path(output_dir).mkdir(parents=True, exist_ok=True)
+        plt.savefig(pathlib.Path(output_dir) / output_fname, format='pdf')
+        plt.close()