Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

amplitude based stimuli in apply_multiple_stimuli #182

Merged
merged 4 commits into from
May 8, 2024
Merged
Show file tree
Hide file tree
Changes from 2 commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
34 changes: 28 additions & 6 deletions bluecellulab/analysis/inject_sequence.py
Original file line number Diff line number Diff line change
Expand Up @@ -74,6 +74,7 @@ def apply_multiple_stimuli(
cell: Cell,
stimulus_name: StimulusName,
amplitudes: Sequence[float],
threshold_based: bool = True,
section_name: str | None = None,
segment: float = 0.5,
n_processes: int | None = None,
Expand All @@ -84,6 +85,8 @@ def apply_multiple_stimuli(
cell: The cell to which the stimuli are applied.
stimulus_name: The name of the stimulus to apply.
amplitudes: The amplitudes of the stimuli to apply.
threshold_based: Whether to consider amplitudes to be
threshold percentages or to be raw amplitudes.
section_name: Section name of the cell where the stimuli are applied.
If None, the stimuli are applied at the soma[0] of the cell.
segment: The segment of the section where the stimuli are applied.
Expand All @@ -103,18 +106,37 @@ def apply_multiple_stimuli(

# Prepare arguments for each stimulus
for amplitude in amplitudes:
if threshold_based:
thres_perc = amplitude
amp = None
else:
thres_perc = None
amp = amplitude

if stimulus_name == StimulusName.AP_WAVEFORM:
stimulus = stim_factory.ap_waveform(threshold_current=cell.threshold, threshold_percentage=amplitude)
stimulus = stim_factory.ap_waveform(
threshold_current=cell.threshold, threshold_percentage=thres_perc, amplitude=amp
)
elif stimulus_name == StimulusName.IDREST:
stimulus = stim_factory.idrest(threshold_current=cell.threshold, threshold_percentage=amplitude)
stimulus = stim_factory.idrest(
threshold_current=cell.threshold, threshold_percentage=thres_perc, amplitude=amp
)
elif stimulus_name == StimulusName.IV:
stimulus = stim_factory.iv(threshold_current=cell.threshold, threshold_percentage=amplitude)
stimulus = stim_factory.iv(
threshold_current=cell.threshold, threshold_percentage=thres_perc, amplitude=amp
)
elif stimulus_name == StimulusName.FIRE_PATTERN:
stimulus = stim_factory.fire_pattern(threshold_current=cell.threshold, threshold_percentage=amplitude)
stimulus = stim_factory.fire_pattern(
threshold_current=cell.threshold, threshold_percentage=thres_perc, amplitude=amp
)
elif stimulus_name == StimulusName.POS_CHEOPS:
stimulus = stim_factory.pos_cheops(threshold_current=cell.threshold, threshold_percentage=amplitude)
stimulus = stim_factory.pos_cheops(
threshold_current=cell.threshold, threshold_percentage=thres_perc, amplitude=amp
)
elif stimulus_name == StimulusName.NEG_CHEOPS:
stimulus = stim_factory.neg_cheops(threshold_current=cell.threshold, threshold_percentage=amplitude)
stimulus = stim_factory.neg_cheops(
threshold_current=cell.threshold, threshold_percentage=thres_perc, amplitude=amp
)
else:
raise ValueError("Unknown stimulus name.")

Expand Down
160 changes: 115 additions & 45 deletions bluecellulab/stimulus/factory.py
Original file line number Diff line number Diff line change
@@ -1,5 +1,6 @@
from __future__ import annotations
from abc import ABC, abstractmethod
from typing import Optional
import matplotlib.pyplot as plt
import numpy as np

Expand Down Expand Up @@ -299,106 +300,167 @@ def ramp(
)

def ap_waveform(
self, threshold_current: float, threshold_percentage: float = 220.0
self,
threshold_current: Optional[float] = None,
threshold_percentage: Optional[float] = 220.0,
amplitude: Optional[float] = None,
) -> Stimulus:
"""Returns the APWaveform Stimulus object, a type of Step stimulus.

Args:
threshold_current: The threshold current of the Cell.
threshold_percentage: Percentage of desired threshold_current amplification.
amplitude: Raw amplitude of input current.
"""
pre_delay = 250.0
duration = 50.0
post_delay = 250.0
return Step.threshold_based(
self.dt,
pre_delay=pre_delay,
duration=duration,
post_delay=post_delay,
threshold_current=threshold_current,
threshold_percentage=threshold_percentage,
)

if threshold_current is not None and threshold_current != 0 and threshold_percentage is not None:
return Step.threshold_based(
self.dt,
pre_delay=pre_delay,
duration=duration,
post_delay=post_delay,
threshold_current=threshold_current,
threshold_percentage=threshold_percentage,
)

if amplitude is not None:
return Step.amplitude_based(
self.dt,
pre_delay=pre_delay,
duration=duration,
post_delay=post_delay,
amplitude=amplitude,
)

raise TypeError("You have to give either threshold_current or amplitude")
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I like to catch these edge cases in the tests. It may look trivial but there are sometimes issues due to those.


def idrest(
self,
threshold_current: float,
threshold_percentage: float = 200.0,
threshold_current: Optional[float] = None,
threshold_percentage: Optional[float] = 200.0,
amplitude: Optional[float] = None,
Comment on lines +350 to +352
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

What if the user mistakenly provides all the 3 of threshold_current, threshold_percentage and amplitude?
Or what if both the amplitude and threshold_current are provided but the threshold_percentage missing?
Can those cases also be covered with logs/exceptions/warnings?

) -> Stimulus:
"""Returns the IDRest Stimulus object, a type of Step stimulus.

Args:
threshold_current: The threshold current of the Cell.
threshold_percentage: Percentage of desired threshold_current amplification.
amplitude: Raw amplitude of input current.
"""
pre_delay = 250.0
duration = 1350.0
post_delay = 250.0
return Step.threshold_based(
self.dt,
pre_delay=pre_delay,
duration=duration,
post_delay=post_delay,
threshold_current=threshold_current,
threshold_percentage=threshold_percentage,
)
if threshold_current is not None and threshold_current != 0 and threshold_percentage is not None:
return Step.threshold_based(
self.dt,
pre_delay=pre_delay,
duration=duration,
post_delay=post_delay,
threshold_current=threshold_current,
threshold_percentage=threshold_percentage,
)

if amplitude is not None:
return Step.amplitude_based(
self.dt,
pre_delay=pre_delay,
duration=duration,
post_delay=post_delay,
amplitude=amplitude,
)

raise TypeError("You have to give either threshold_current or amplitude")

def iv(
self,
threshold_current: float,
threshold_percentage: float = -40.0,
threshold_current: Optional[float] = None,
threshold_percentage: Optional[float] = -40.0,
amplitude: Optional[float] = None,
) -> Stimulus:
"""Returns the IV Stimulus object, a type of Step stimulus.

Args:
threshold_current: The threshold current of the Cell.
threshold_percentage: Percentage of desired threshold_current amplification.
amplitude: Raw amplitude of input current.
"""
pre_delay = 250.0
duration = 3000.0
post_delay = 250.0
return Step.threshold_based(
self.dt,
pre_delay=pre_delay,
duration=duration,
post_delay=post_delay,
threshold_current=threshold_current,
threshold_percentage=threshold_percentage,
)
if threshold_current is not None and threshold_current != 0 and threshold_percentage is not None:
return Step.threshold_based(
self.dt,
pre_delay=pre_delay,
duration=duration,
post_delay=post_delay,
threshold_current=threshold_current,
threshold_percentage=threshold_percentage,
)

if amplitude is not None:
return Step.amplitude_based(
self.dt,
pre_delay=pre_delay,
duration=duration,
post_delay=post_delay,
amplitude=amplitude,
)

raise TypeError("You have to give either threshold_current or amplitude")

def fire_pattern(
self,
threshold_current: float,
threshold_percentage: float = 200.0,
threshold_current: Optional[float] = None,
threshold_percentage: Optional[float] = 200.0,
amplitude: Optional[float] = None,
) -> Stimulus:
"""Returns the FirePattern Stimulus object, a type of Step stimulus.

Args:
threshold_current: The threshold current of the Cell.
threshold_percentage: Percentage of desired threshold_current amplification.
amplitude: Raw amplitude of input current.
"""
pre_delay = 250.0
duration = 3600.0
post_delay = 250.0
return Step.threshold_based(
self.dt,
pre_delay=pre_delay,
duration=duration,
post_delay=post_delay,
threshold_current=threshold_current,
threshold_percentage=threshold_percentage,
)
if threshold_current is not None and threshold_current != 0 and threshold_percentage is not None:
return Step.threshold_based(
self.dt,
pre_delay=pre_delay,
duration=duration,
post_delay=post_delay,
threshold_current=threshold_current,
threshold_percentage=threshold_percentage,
)

if amplitude is not None:
return Step.amplitude_based(
self.dt,
pre_delay=pre_delay,
duration=duration,
post_delay=post_delay,
amplitude=amplitude,
)

raise TypeError("You have to give either threshold_current or amplitude")

def pos_cheops(
self,
threshold_current: float,
threshold_percentage: float = 300.0,
threshold_current: Optional[float] = None,
threshold_percentage: Optional[float] = 300.0,
amplitude: Optional[float] = None,
) -> Stimulus:
"""A combination of pyramid shaped Ramp stimuli with a positive
amplitude.

Args:
threshold_current: The threshold current of the Cell.
threshold_percentage: Percentage of desired threshold_current amplification.
amplitude: Raw amplitude of input current.
"""
delay = 250.0
ramp1_duration = 4000.0
Expand All @@ -407,7 +469,10 @@ def pos_cheops(
inter_delay = 2000.0
post_delay = 250.0

amplitude = threshold_current * threshold_percentage / 100
if amplitude is None:
if threshold_current is None or threshold_current == 0 or threshold_percentage is None:
raise TypeError("You have to give either threshold_current or amplitude")
amplitude = threshold_current * threshold_percentage / 100
result = (
Empty(self.dt, duration=delay)
+ Slope(self.dt, duration=ramp1_duration, amplitude_start=0.0, amplitude_end=amplitude)
Expand All @@ -424,15 +489,17 @@ def pos_cheops(

def neg_cheops(
self,
threshold_current: float,
threshold_percentage: float = 300.0,
threshold_current: Optional[float] = None,
threshold_percentage: Optional[float] = 300.0,
amplitude: Optional[float] = None,
) -> Stimulus:
"""A combination of pyramid shaped Ramp stimuli with a negative
amplitude.

Args:
threshold_current: The threshold current of the Cell.
threshold_percentage: Percentage of desired threshold_current amplification.
amplitude: Raw amplitude of input current.
"""
delay = 1750.0
ramp1_duration = 3333.0
Expand All @@ -441,7 +508,10 @@ def neg_cheops(
inter_delay = 2000.0
post_delay = 250.0

amplitude = - threshold_current * threshold_percentage / 100
if amplitude is None:
if threshold_current is None or threshold_current == 0 or threshold_percentage is None:
raise TypeError("You have to give either threshold_current or amplitude")
amplitude = - threshold_current * threshold_percentage / 100
result = (
Empty(self.dt, duration=delay)
+ Slope(self.dt, duration=ramp1_duration, amplitude_start=0.0, amplitude_end=amplitude)
Expand Down
12 changes: 9 additions & 3 deletions tests/test_analysis/test_inject_sequence.py
Original file line number Diff line number Diff line change
Expand Up @@ -32,16 +32,19 @@ def mock_run_stimulus():
def test_apply_multiple_step_stimuli(mock_run_stimulus):
"""Do not run the code in parallel, mock the return value via MockRecording."""
amplitudes = [80, 100, 120, 140]
thres_perc = [0.08, 0.1, 0.12, 0.14]
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

For test time maybe instead of 4, just 1 or 2 values can be kept here

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

In the future more protocols and therefore more tests will be added, so the tests will take even longer. Maybe after running the protocols once, the other tests can just mock the expensive neuron simulation run part to retrieve precomputed (or even constant) voltage results to us. That way the newly added code still can be tested.

cell = create_ball_stick()

with patch('bluecellulab.analysis.inject_sequence.IsolatedProcess') as mock_isolated_process, \
patch('bluecellulab.analysis.inject_sequence.run_stimulus', mock_run_stimulus):
# the mock process pool to return a list of MockRecordings
mock_isolated_process.return_value.__enter__.return_value.starmap.return_value = [MockRecording() for _ in amplitudes]

recordings = apply_multiple_stimuli(cell, StimulusName.FIRE_PATTERN, amplitudes, n_processes=4)
recordings = apply_multiple_stimuli(cell, StimulusName.FIRE_PATTERN, amplitudes, threshold_based=False, n_processes=4)
recordings_thres = apply_multiple_stimuli(cell, StimulusName.FIRE_PATTERN, thres_perc, n_processes=4)
assert len(recordings) == len(amplitudes)
for recording in recordings.values():
assert len(recordings_thres) == len(thres_perc)
for recording in list(recordings.values()) + list(recordings_thres.values()):
assert len(recording.time) > 0
assert len(recording.time) == len(recording.voltage)
assert len(recording.time) == len(recording.current)
Expand All @@ -52,7 +55,10 @@ def test_apply_multiple_step_stimuli(mock_run_stimulus):
assert "Unknown stimulus name" in str(exc_info.value)

short_amplitudes = [80]
short_thres = [0.08]
other_stim = [StimulusName.AP_WAVEFORM, StimulusName.IV, StimulusName.IDREST, StimulusName.POS_CHEOPS, StimulusName.NEG_CHEOPS]
for stim in other_stim:
res = apply_multiple_stimuli(cell, stim, short_amplitudes, n_processes=1)
res = apply_multiple_stimuli(cell, stim, short_amplitudes, threshold_based=False, n_processes=1)
res_thres = apply_multiple_stimuli(cell, stim, short_thres, n_processes=1)
assert len(res) == len(short_amplitudes)
assert len(res_thres) == len(short_thres)