Seclamp stimulus

bluecellulab/cell/injector.py

@@ -142,7 +142,7 @@ def add_ramp(
         return tstim
 
     def add_voltage_clamp(
-            self, stop_time, level, rs=None, section=None, segx=0.5,
+            self, stop_time, level, durations=None, levels=None, rs=None, section=None, segx=0.5,
             current_record_name=None, current_record_dt=None):
         """Add a voltage clamp.
 
@@ -153,6 +153,10 @@ def add_voltage_clamp(
             Time at which voltage clamp should stop
         level : float
             Voltage level of the vc (in mV)
+        durations: list of float
+            Durations of each step of the vc (in ms)
+        levels : list of float
+            Voltage levels of the vc (in mV)
         rs: float
             Series resistance of the vc (in MOhm)
         section: NEURON object
@@ -169,6 +173,7 @@ def add_voltage_clamp(
 
         SEClamp (NEURON) object of the created vc
         """
+        from neuron import h  # noqa: PLC0415
 
         if section is None:
             section = self.soma
@@ -177,12 +182,29 @@ def add_voltage_clamp(
         vclamp = neuron.h.SEClamp(segx, sec=section)
         self.persistent.append(vclamp)
 
-        vclamp.amp1 = level
-        vclamp.dur1 = stop_time
-
         if rs is not None:
             vclamp.rs = rs
 
+        vclamp.dur1 = stop_time
+        vclamp.amp1 = level
+
+        if durations is not None and levels is not None:
+            if len(levels) != len(durations) - 1:
+                raise BluecellulabError("Inconsistent durations and levels for seclamp.")
+
+            voltage_vec = h.Vector(levels)
+            time_vec = h.Vector(np.cumsum(durations))
+
+            self.persistent.append(time_vec)
+            self.persistent.append(voltage_vec)
+
+            voltage_vec.play(
+                vclamp._ref_amp1,  # noqa: SLF001
+                time_vec,
+                0,
+                sec=section,
+            )
+
         current = neuron.h.Vector()
         if current_record_dt is None:
             current.record(vclamp._ref_i)
@@ -537,3 +559,15 @@ def add_sinusoidal(self, stimulus) -> TStim:
             stimulus.duration,
             stimulus.frequency,
         )
+
+    def add_seclamp(self, stimulus, section=None, segx=0.5):
+        """Add a SEClamp stimulus."""
+        return self.add_voltage_clamp(
+            stimulus.duration,
+            stimulus.voltage,
+            stimulus.durations,
+            stimulus.voltages,
+            rs=stimulus.series_resistance,
+            section=section,
+            segx=segx,
+        )

bluecellulab/circuit_simulation.py

@@ -142,6 +142,7 @@ def instantiate_gids(
         add_ornstein_uhlenbeck_stimuli: bool = False,
         add_sinusoidal_stimuli: bool = False,
         add_linear_stimuli: bool = False,
+        add_seclamp_stimuli: bool = False,
     ):
         """Instantiate a list of cells.
 
@@ -222,6 +223,11 @@ def instantiate_gids(
                                 Setting add_stimuli=True,
                                 will automatically set this option to
                                 True.
+        add_seclamp_stimuli : Process the 'seclamp' stimuli
+                                blocks of the simulation config.
+                                Setting add_stimuli=True,
+                                will automatically set this option to
+                                True.
         """
         if not isinstance(cells, list):
             cells = [cells]
@@ -282,6 +288,7 @@ def instantiate_gids(
             add_shotnoise_stimuli = True
             add_ornstein_uhlenbeck_stimuli = True
             add_linear_stimuli = True
+            add_seclamp_stimuli = True
 
         if add_noise_stimuli or \
                 add_hyperpolarizing_stimuli or \
@@ -290,7 +297,8 @@ def instantiate_gids(
                 add_shotnoise_stimuli or \
                 add_ornstein_uhlenbeck_stimuli or \
                 add_sinusoidal_stimuli or \
-                add_linear_stimuli:
+                add_linear_stimuli or \
+                add_seclamp_stimuli:
             self._add_stimuli(
                 add_noise_stimuli=add_noise_stimuli,
                 add_hyperpolarizing_stimuli=add_hyperpolarizing_stimuli,
@@ -299,7 +307,8 @@ def instantiate_gids(
                 add_shotnoise_stimuli=add_shotnoise_stimuli,
                 add_ornstein_uhlenbeck_stimuli=add_ornstein_uhlenbeck_stimuli,
                 add_sinusoidal_stimuli=add_sinusoidal_stimuli,
-                add_linear_stimuli=add_linear_stimuli
+                add_linear_stimuli=add_linear_stimuli,
+                add_seclamp_stimuli=add_seclamp_stimuli,
             )
 
         configure_all_reports(
@@ -319,7 +328,8 @@ def _add_stimuli(self, add_noise_stimuli=False,
                      add_shotnoise_stimuli=False,
                      add_ornstein_uhlenbeck_stimuli=False,
                      add_sinusoidal_stimuli=False,
-                     add_linear_stimuli=False
+                     add_linear_stimuli=False,
+                     add_seclamp_stimuli=False,
                      ) -> None:
         """Instantiate all the stimuli."""
         stimuli_entries = self.circuit_access.config.get_all_stimuli_entries()
@@ -389,6 +399,9 @@ def _add_stimuli(self, add_noise_stimuli=False,
                 elif isinstance(stimulus, circuit_stimulus_definitions.Sinusoidal):
                     if add_sinusoidal_stimuli:
                         self.cells[cell_id].add_sinusoidal(stimulus)
+                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.SynapseReplay):  # sonata only
                     if self.circuit_access.target_contains_cell(
                         stimulus.target, cell_id

bluecellulab/stimulus/circuit_stimulus_definitions.py

@@ -50,6 +50,7 @@ class Pattern(Enum):
     ORNSTEIN_UHLENBECK = "ornstein_uhlenbeck"
     RELATIVE_ORNSTEIN_UHLENBECK = "relative_ornstein_uhlenbeck"
     SINUSOIDAL = "sinusoidal"
+    SECLAMP = "seclamp"
 
     @classmethod
     def from_blueconfig(cls, pattern: str) -> Pattern:
@@ -98,6 +99,8 @@ def from_sonata(cls, pattern: str) -> Pattern:
             return Pattern.RELATIVE_ORNSTEIN_UHLENBECK
         elif pattern == "sinusoidal":
             return Pattern.SINUSOIDAL
+        elif pattern == "seclamp":
+            return Pattern.SECLAMP
         else:
             raise ValueError(f"Unknown pattern {pattern}")
 
@@ -374,6 +377,18 @@ def from_sonata(cls, stimulus_entry: dict, config_dir: Optional[str] = None) ->
                 node_set=node_set,
                 compartment_set=compartment_set,
             )
+        elif pattern == Pattern.SECLAMP:
+            return SEClamp(
+                target=target_name,
+                delay=stimulus_entry["delay"],
+                duration=stimulus_entry["duration"],
+                voltage=stimulus_entry["voltage"],
+                durations=stimulus_entry.get("duration_levels", None),
+                voltages=stimulus_entry.get("voltage_levels", None),
+                series_resistance=stimulus_entry.get("series_resistance", 0.01),
+                node_set=node_set,
+                compartment_set=compartment_set,
+            )
         else:
             raise ValueError(f"Unknown pattern {pattern}")
 
@@ -497,3 +512,11 @@ class RelativeOrnsteinUhlenbeck(Stimulus):
 class Sinusoidal(Stimulus):
     amp_start: float
     frequency: float
+
+
+@dataclass(frozen=True, config=dict(extra="forbid"))
+class SEClamp(Stimulus):
+    voltage: float
+    durations: Optional[list[float]]
+    voltages: Optional[list[float]]
+    series_resistance: float

tests/test_cell/test_injector.py

@@ -112,7 +112,7 @@ def test_voltage_clamp(self):
         assert seclamp_obj.rs == rs
 
         self.sim.run(10, dt=1, cvode=False)
-        current = self.cell.get_recording("test_volt_clamp")
+        current = self.cell.get_recording(n_recording)
         assert current == approx(np.array(
             [66.5, 5.39520998, -10.76796553, 20.6887735,
              17.8876999, 15.14995787, 13.47384441, 12.55945316,
@@ -134,10 +134,39 @@ def test_voltage_clamp_dt(self):
         assert seclamp_obj.rs == rs
 
         self.sim.run(10, dt=1, cvode=False)
-        current = self.cell.get_recording("test_volt_clamp_dt")
+        current = self.cell.get_recording(n_recording)
         assert current == approx(np.array(
             [66.5, -10.76796553, 17.8876999, 13.47384441, 12.09052411]), abs=1e-3)
 
+    def test_multilevel_voltage_clamp(self):
+        """Test voltage clamp with multiple voltages."""
+        amp = 10
+        stop_time = 10
+        voltages = [20, 30]
+        durations = [2, 3, 5]
+        rs = 1
+        seclamp_obj = self.cell.add_voltage_clamp(
+            stop_time=stop_time, level=amp, rs=rs, levels=voltages, durations=durations,
+        )
+        assert seclamp_obj.amp1 == amp
+        assert seclamp_obj.dur1 == stop_time
+        assert seclamp_obj.rs == rs
+
+        v_rec = neuron.h.Vector()
+        v_rec.record(seclamp_obj._ref_vc)
+
+        self.sim.run(10, dt=0.2, cvode=False)
+        voltage_rec = np.array(v_rec.to_python())
+        # has to add one 10 at the beginning for t=0
+        assert (voltage_rec == np.array(
+            [
+                10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10,  # 10 + 1 extra for t=0
+                20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20,  # 15
+                30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30,
+                30, 30, 30, 30, 30, 30, 30, 30, 30, 30  # 25
+            ]
+        )).all()
+
     def test_get_noise_step_rand(self):
         """Unit test for _get_noise_step_rand."""
         noisestim_count = 5

tests/test_stimulus/test_circuit_stimulus_definitions.py

@@ -30,6 +30,7 @@ def test_pattern_from_sonata_valid():
         "relative_shot_noise": Pattern.RELATIVE_SHOT_NOISE,
         "ornstein_uhlenbeck": Pattern.ORNSTEIN_UHLENBECK,
         "relative_ornstein_uhlenbeck": Pattern.RELATIVE_ORNSTEIN_UHLENBECK,
+        "seclamp": Pattern.SECLAMP,
     }
 
     for sonata_pattern, expected_enum in valid_patterns.items():