Add subthreshold stimulus support

bluecellulab/cell/injector.py

@@ -38,6 +38,7 @@
     ShotNoise,
     RelativeOrnsteinUhlenbeck,
     RelativeShotNoise,
+    SubThreshold,
 )
 from bluecellulab.type_aliases import NeuronSection, TStim
 
@@ -331,6 +332,20 @@ def add_replay_relativelinear(self, stimulus: RelativeLinear, section=None, segx
 
         return tstim
 
+    def add_replay_subthreshold(self, stimulus: SubThreshold, section=None, segx=0.5):
+        """Inject a current step at some percent below the cell's threshold.
+
+        The injected amplitude is: threshold * (100 - percent_less) / 100
+        This matches the Neurodamus SubThreshold implementation.
+        """
+        if section is None:
+            section = self.soma  # type: ignore
+        amp = self.threshold * (100 - stimulus.percent_less) / 100  # type: ignore
+        tstim = neuron.h.TStim(segx, sec=section)  # type: ignore
+        tstim.pulse(stimulus.delay, stimulus.duration, amp)
+        self.persistent.append(tstim)  # type: ignore
+        return tstim
+
     def _get_ornstein_uhlenbeck_rand(self, stim_count, seed):
         """Return rng for ornstein_uhlenbeck simulation."""
         rng_settings = RNGSettings.get_instance()

bluecellulab/circuit_simulation.py

@@ -180,6 +180,7 @@ def instantiate_gids(
         add_sinusoidal_stimuli: bool = False,
         add_linear_stimuli: bool = False,
         add_seclamp_stimuli: bool = False,
+        add_subthreshold_stimuli: bool = False,
     ):
         """Instantiate a list of cells.
 
@@ -268,6 +269,11 @@ def instantiate_gids(
                                 Setting add_stimuli=True,
                                 will automatically set this option to
                                 True.
+        add_subthreshold_stimuli : Process the 'subthreshold' stimuli
+                                blocks of the simulation config.
+                                Setting add_stimuli=True,
+                                will automatically set this option to
+                                True.
         """
         if not isinstance(cells, list):
             cells = [cells]
@@ -354,6 +360,7 @@ def instantiate_gids(
             add_ornstein_uhlenbeck_stimuli = True
             add_linear_stimuli = True
             add_seclamp_stimuli = True
+            add_subthreshold_stimuli = True
 
         if (
             add_noise_stimuli
@@ -365,6 +372,7 @@ def instantiate_gids(
             or add_sinusoidal_stimuli
             or add_linear_stimuli
             or add_seclamp_stimuli
+            or add_subthreshold_stimuli
         ):
             self._add_stimuli(
                 add_noise_stimuli=add_noise_stimuli,
@@ -376,6 +384,7 @@ def instantiate_gids(
                 add_sinusoidal_stimuli=add_sinusoidal_stimuli,
                 add_linear_stimuli=add_linear_stimuli,
                 add_seclamp_stimuli=add_seclamp_stimuli,
+                add_subthreshold_stimuli=add_subthreshold_stimuli,
             )
 
         self.recording_index, self.sites_index = prepare_recordings_for_reports(
@@ -403,6 +412,7 @@ def _add_stimuli(
         add_sinusoidal_stimuli=False,
         add_linear_stimuli=False,
         add_seclamp_stimuli=False,
+        add_subthreshold_stimuli=False,
     ) -> None:
         """Instantiate all the stimuli."""
         stimuli_entries = self.circuit_access.config.get_all_stimuli_entries()
@@ -515,6 +525,11 @@ def _add_stimuli(
                 elif isinstance(stimulus, circuit_stimulus_definitions.SEClamp):  # sonata only
                     if add_seclamp_stimuli:
                         self.cells[cell_id].add_seclamp(stimulus, section=sec, segx=segx)
+                elif isinstance(stimulus, circuit_stimulus_definitions.SubThreshold):
+                    if add_subthreshold_stimuli:
+                        self.cells[cell_id].add_replay_subthreshold(
+                            stimulus, section=sec, segx=segx
+                        )
                 elif isinstance(
                     stimulus, circuit_stimulus_definitions.SynapseReplay
                 ):  # sonata only

bluecellulab/stimulus/circuit_stimulus_definitions.py

@@ -51,6 +51,7 @@ class Pattern(Enum):
     RELATIVE_ORNSTEIN_UHLENBECK = "relative_ornstein_uhlenbeck"
     SINUSOIDAL = "sinusoidal"
     SECLAMP = "seclamp"
+    SUBTHRESHOLD = "subthreshold"
 
     @classmethod
     def from_blueconfig(cls, pattern: str) -> Pattern:
@@ -72,6 +73,8 @@ def from_blueconfig(cls, pattern: str) -> Pattern:
             return Pattern.ORNSTEIN_UHLENBECK
         elif pattern == "RelativeOrnsteinUhlenbeck":
             return Pattern.RELATIVE_ORNSTEIN_UHLENBECK
+        elif pattern == "SubThreshold":
+            return Pattern.SUBTHRESHOLD
         else:
             raise ValueError(f"Unknown pattern {pattern}")
 
@@ -101,6 +104,8 @@ def from_sonata(cls, pattern: str) -> Pattern:
             return Pattern.SINUSOIDAL
         elif pattern == "seclamp":
             return Pattern.SECLAMP
+        elif pattern == "subthreshold":
+            return Pattern.SUBTHRESHOLD
         else:
             raise ValueError(f"Unknown pattern {pattern}")
 
@@ -230,6 +235,15 @@ def from_blueconfig(cls, stimulus_entry: dict) -> Optional[Stimulus]:
                 node_set=stimulus_entry["Target"],
                 compartment_set=None,
             )
+        elif pattern == Pattern.SUBTHRESHOLD:
+            return SubThreshold(
+                target=stimulus_entry["Target"],
+                delay=stimulus_entry["Delay"],
+                duration=stimulus_entry["Duration"],
+                percent_less=stimulus_entry["PercentLess"],
+                node_set=stimulus_entry["Target"],
+                compartment_set=None,
+            )
         else:
             raise ValueError(f"Unknown pattern {pattern}")
 
@@ -395,6 +409,15 @@ def from_sonata(cls, stimulus_entry: dict, config_dir: Optional[str] = None) ->
                 node_set=node_set,
                 compartment_set=compartment_set,
             )
+        elif pattern == Pattern.SUBTHRESHOLD:
+            return SubThreshold(
+                target=target_name,
+                delay=stimulus_entry["delay"],
+                duration=stimulus_entry["duration"],
+                percent_less=stimulus_entry["percent_less"],
+                node_set=node_set,
+                compartment_set=compartment_set,
+            )
         else:
             raise ValueError(f"Unknown pattern {pattern}")
 
@@ -534,3 +557,9 @@ class SEClamp(Stimulus):
     durations: Optional[list[float]]
     voltages: Optional[list[float]]
     series_resistance: float
+
+
+@dataclass(frozen=True, config=dict(extra="forbid"))
+class SubThreshold(Stimulus):
+    """Injects a current step at some percent below a cell's threshold."""
+    percent_less: float

tests/test_cell/test_injector.py

@@ -286,6 +286,20 @@ def test_add_replay_relativelinear(self):
         assert tstim.stim.to_python() == [0.0, 0.0, percent_60, percent_100, 0.0, 0.0]
         assert tstim.tvec.to_python() == [0.0, 0.0, 0.0, 20.0, 20.0, 20.0]
 
+    def test_add_replay_subthreshold(self):
+        """Unit test for add_replay_subthreshold."""
+        from bluecellulab.stimulus.circuit_stimulus_definitions import SubThreshold
+
+        stimulus = SubThreshold(
+            target="single-cell",
+            delay=5, duration=15, percent_less=20)
+        tstim = self.cell.add_replay_subthreshold(stimulus)
+
+        # Expected amplitude: threshold * (100 - 20) / 100 = threshold * 0.8
+        expected_amp = self.cell.threshold * 0.8
+        assert tstim.stim.to_python() == approx([0.0, expected_amp, expected_amp, 0.0, 0.0])
+        assert tstim.tvec.to_python() == approx([5.0, 5.0, 20.0, 20.0, 20.0])
+
     def test_get_ornstein_uhlenbeck_rand(self):
         """Unit test to check RNG generated for ornstein_uhlenbeck."""
         rng = self.cell._get_ornstein_uhlenbeck_rand(0, 144)

tests/test_stimulus/test_circuit_simulation_stimuli.py

@@ -18,7 +18,7 @@
 import pytest
 
 from bluecellulab.circuit_simulation import CircuitSimulation
-from bluecellulab.stimulus.circuit_stimulus_definitions import Noise, Pulse
+from bluecellulab.stimulus.circuit_stimulus_definitions import Noise, Pulse, SubThreshold
 
 
 class FakeCell:
@@ -42,6 +42,9 @@ def add_replay_noise(self, stimulus, noise_seed=None, noisestim_count=0, section
     def add_pulse(self, stimulus, section=None, segx=0.5):
         self.calls.append(("pulse", section, segx, stimulus))
 
+    def add_replay_subthreshold(self, stimulus, section=None, segx=0.5):
+        self.calls.append(("subthreshold", section, segx, stimulus))
+
 
 class FakeConfig:
     def __init__(self, stimuli, compartment_sets=None):
@@ -446,3 +449,31 @@ def test_compartment_set_applied_only_to_matching_population():
     assert cell_a.calls[0][0] == "noise"
 
     assert cell_b.calls == []
+
+
+def test_subthreshold_stimulus_dispatched_to_cell():
+    """SubThreshold stimulus should be dispatched via add_replay_subthreshold."""
+    subthreshold = SubThreshold(
+        target="node_set_A",
+        delay=0.0,
+        duration=1000.0,
+        percent_less=20.0,
+        node_set="node_set_A",
+    )
+
+    cfg = FakeConfig([subthreshold], compartment_sets={})
+    access = FakeCircuitAccess(cfg, target_map={"node_set_A": [1]})
+
+    cell = FakeCell(1)
+    cells = {1: cell}
+
+    sim = make_dummy_sim(access, cells)
+
+    CircuitSimulation._add_stimuli(sim, add_subthreshold_stimuli=True)
+
+    assert len(cell.calls) == 1
+    call = cell.calls[0]
+    assert call[0] == "subthreshold"
+    assert call[1] == cell.soma
+    assert call[2] == 0.5
+    assert call[3].percent_less == 20.0

tests/test_stimulus/test_circuit_stimulus_definitions.py

@@ -14,7 +14,7 @@
 # limitations under the License.
 
 import pytest
-from bluecellulab.stimulus.circuit_stimulus_definitions import Noise, Pattern, Stimulus
+from bluecellulab.stimulus.circuit_stimulus_definitions import Noise, Pattern, Stimulus, SubThreshold
 
 
 def test_pattern_from_sonata_valid():
@@ -31,6 +31,7 @@ def test_pattern_from_sonata_valid():
         "ornstein_uhlenbeck": Pattern.ORNSTEIN_UHLENBECK,
         "relative_ornstein_uhlenbeck": Pattern.RELATIVE_ORNSTEIN_UHLENBECK,
         "seclamp": Pattern.SECLAMP,
+        "subthreshold": Pattern.SUBTHRESHOLD,
     }
 
     for sonata_pattern, expected_enum in valid_patterns.items():
@@ -73,3 +74,54 @@ def test_from_sonata_noise_requires_one_mean_field():
 
     with pytest.raises(ValueError, match="Noise input must contain exactly one of 'mean' or 'mean_percent'."):
         Stimulus.from_sonata({**base, "mean": 0.01, "mean_percent": 5.0})
+
+
+def test_pattern_from_blueconfig_subthreshold():
+    """Test SubThreshold mapping from BlueConfig."""
+    assert Pattern.from_blueconfig("SubThreshold") == Pattern.SUBTHRESHOLD
+
+
+def test_subthreshold_from_sonata():
+    """Test parsing SubThreshold stimulus from SONATA config."""
+    entry = {
+        "module": "subthreshold",
+        "delay": 10.0,
+        "duration": 500.0,
+        "percent_less": 20.0,
+        "node_set": "Excitatory",
+    }
+    stim = Stimulus.from_sonata(entry)
+    assert isinstance(stim, SubThreshold)
+    assert stim.percent_less == 20.0
+    assert stim.delay == 10.0
+    assert stim.duration == 500.0
+    assert stim.target == "Excitatory"
+
+
+def test_subthreshold_from_blueconfig():
+    """Test parsing SubThreshold stimulus from BlueConfig."""
+    entry = {
+        "Pattern": "SubThreshold",
+        "Target": "Mosaic",
+        "Delay": 0.0,
+        "Duration": 1000.0,
+        "PercentLess": 15.0,
+    }
+    stim = Stimulus.from_blueconfig(entry)
+    assert isinstance(stim, SubThreshold)
+    assert stim.percent_less == 15.0
+    assert stim.target == "Mosaic"
+
+
+def test_subthreshold_dataclass():
+    """Test SubThreshold dataclass creation."""
+    stim = SubThreshold(
+        target="All",
+        delay=0.0,
+        duration=100.0,
+        percent_less=10.0,
+        node_set="All",
+    )
+    assert stim.percent_less == 10.0
+    assert stim.delay == 0.0
+    assert stim.duration == 100.0