[BBPBGLIB-1137] Account for nonzero total current when using SEClamp; allow nonzero total current when simulating physical electrodes

docs/online-lfp.rst

@@ -130,6 +130,12 @@ Subsequently, an ERROR will be encountered when instantiating the LFP report:
 
 To ensure accurate and valid LFP reports, make sure that the electrodes file corresponds to the circuit being used in your simulation.
 
+- **Stimulus Electrode Compatibility**: A common use case is that current will be injected into a population to account for synaptic inputs from neural populations that are not modeled. In this case,the injected current should be considered a membrane current rather than an electrode current, and it is neccessary that total current over the neuron sums to zero in order to produce valid extracellular recording results. The Neuron SEClamp class does not fulfill these criteria due to numerical issues. We have created a new point process, `ConductanceSource`, available in `neurodamus-neocortex`, which does fulfill the criteria. If an conductance source stimulus is present in the simulation config file, `ConductanceSource` will be used by default instead of the SEClamp mechanism. The injected current will be reported as part of the `i_membrane` variable, rather than as an electrode current. 
+
+However, it may be the case that the user wishes to model a physical electrode, rather than missing synaptic input, using the conductance source mechanism. In this case, the total current over the neuron is nonzero, and the injected current should not be considered a membrane current. For this reason, we have added the key `represents_phsyical_electrode` to the stimulus block. With the key-value pair `represents_physical_electrode:true`, SEClamp will be used rather than ConductanceSource.
+
+Similarly, current sources may also be used to model the effects of missing synaptic inputs. We have created a new point process, `MembraneCurrentSource`, which is used instead of IClamp if the key `represents_phsyical_electrode` is set to false or is not set. `MembraneCurrentSource` behaves identically to IClamp, but is considered a membrane current, and is therefore accounted for in the calculation of the extracellular signal. It is not reported on as an electrode current. Setting `represents_physical_electrode:true` will result in using IClamp instead of `MembraneCurrentSource` 
+
 By keeping these considerations in mind, you can ensure a smooth and successful usage of the online LFP calculation feature.
 
 Conclusion

neurodamus/core/stimuli.py

@@ -10,19 +10,22 @@
 
 class SignalSource:
 
-    def __init__(self, base_amp=0.0, *, delay=0, rng=None):
+    def __init__(self, base_amp=0.0, *, delay=0, rng=None, represents_physical_electrode=False):
         """
         Creates a new signal source, which can create composed signals
         Args:
             base_amp: The base (resting) amplitude of the signal (Default: 0)
             rng: The Random Number Generator. Used in the Noise functions
+            represents_physical_electrode: Whether the source represents a phsyical
+            electrode or missing synaptic input
         """
         h = Neuron.h
         self.stim_vec = h.Vector()
         self.time_vec = h.Vector()
         self._cur_t = 0
         self._base_amp = base_amp
         self._rng = rng
+        self._represents_physical_electrode = represents_physical_electrode
         if delay > .0:
             self._add_point(base_amp)
             self._cur_t = delay
@@ -356,7 +359,7 @@ def shot_noise(cls, tau_D, tau_R, rate, amp_mean, var, duration, dt=0.25, base_a
 
     @classmethod
     def ornstein_uhlenbeck(cls, tau, sigma, mean, duration, dt=0.25, base_amp=.0, **kw):
-        return cls(base_amp, **kw).add_ornstein_uhlenbeck(tau, sigma, mean,duration, dt)
+        return cls(base_amp, **kw).add_ornstein_uhlenbeck(tau, sigma, mean, duration, dt)
 
     # Operations
     def __add__(self, other):
@@ -367,19 +370,27 @@ def __add__(self, other):
 class CurrentSource(SignalSource):
     _all_sources = []
 
-    def __init__(self, base_amp=0.0, *, delay=0, rng=None):
+    def __init__(self, base_amp=0.0, *, delay=0, rng=None, physical_electrode=False):
         """
         Creates a new current source that injects a signal under IClamp
         """
-        super().__init__(base_amp, delay=delay, rng=rng)
+        super().__init__(base_amp, delay=delay, rng=rng,
+                         represents_physical_electrode=physical_electrode)
         self._clamps = set()
         self._all_sources.append(self)
 
     class _Clamp:
         def __init__(self, cell_section, position=0.5, clamp_container=None,
                      stim_vec_mode=True, time_vec=None, stim_vec=None,
                      **clamp_params):
-            self.clamp = Neuron.h.IClamp(position, sec=cell_section)
+            represents_physical_electrode = clamp_params.get('represents_physical_electrode', False)
+            # Checks if new MembraneCurrentSource mechanism is available and if source does not
+            # represent physical electrode, otherwise fall back to IClamp.
+            if not represents_physical_electrode and hasattr(Neuron.h, "MembraneCurrentSource"):
+                self.clamp = Neuron.h.MembraneCurrentSource(position, sec=cell_section)
+            else:
+                self.clamp = Neuron.h.IClamp(position, sec=cell_section)
+
             if stim_vec_mode:
                 assert time_vec is not None and stim_vec is not None
                 self.clamp.dur = time_vec[-1]
@@ -397,8 +408,9 @@ def detach(self):
             del self.clamp  # Force del on the clamp (there might be references to self)
 
     def attach_to(self, section, position=0.5):
-        return CurrentSource._Clamp(section, position, self._clamps, True,
-                                    self.time_vec, self.stim_vec)
+        return CurrentSource._Clamp(
+            section, position, self._clamps, True, self.time_vec, self.stim_vec,
+            represents_physical_electrode=self._represents_physical_electrode)
 
     # Constant has a special attach_to and doesnt share any composing method
     class Constant:
@@ -418,14 +430,16 @@ def attach_to(self, section, position=0.5):
 class ConductanceSource(SignalSource):
     _all_sources = []
 
-    def __init__(self, reversal=0.0, *, delay=.0, rng=None):
+    def __init__(self, reversal=0.0, *, delay=.0, rng=None, physical_electrode=False):
         """
         Creates a new conductance source that injects a conductance by driving
         the rs of an SEClamp at a given reversal potential.
 
         reversal: reversal potential of conductance (mV)
         """
-        super().__init__(0.0, delay=delay, rng=rng)  # set SignalSource's base_amp to zero
+        # set SignalSource's base_amp to zero
+        super().__init__(reversal, delay=delay, rng=rng,
+                         represents_physical_electrode=physical_electrode)
         self._reversal = reversal   # set reversal from base_amp parameter in classmethods
         self._clamps = set()
         self._all_sources.append(self)
@@ -434,7 +448,14 @@ class _DynamicClamp:
         def __init__(self, cell_section, position=0.5, clamp_container=None,
                      stim_vec_mode=True, time_vec=None, stim_vec=None,
                      reversal=0.0, **clamp_params):
-            self.clamp = Neuron.h.SEClamp(position, sec=cell_section)
+            represents_physical_electrode = clamp_params.get('represents_physical_electrode', False)
+            # Checks if new conductanceSource mechanism is available and if source does not
+            # represent physical electrode, otherwise fall back to SEClamp.
+            if not represents_physical_electrode and hasattr(Neuron.h, "ConductanceSource"):
+                self.clamp = Neuron.h.ConductanceSource(position, sec=cell_section)
+            else:
+                self.clamp = Neuron.h.SEClamp(position, sec=cell_section)
+
             if stim_vec_mode:
                 assert time_vec is not None and stim_vec is not None
                 self.clamp.dur1 = time_vec[-1]
@@ -460,8 +481,9 @@ def detach(self):
             del self.clamp  # Force del on the clamp (there might be references to self)
 
     def attach_to(self, section, position=0.5):
-        return ConductanceSource._DynamicClamp(section, position, self._clamps, True,
-                                               self.time_vec, self.stim_vec, self._reversal)
+        return ConductanceSource._DynamicClamp(
+            section, position, self._clamps, True, self.time_vec, self.stim_vec,
+            self._reversal, represents_physical_electrode=self._represents_physical_electrode)
 
 
 # EStim class is a derivative of TStim for stimuli with an extracelular electrode. The main

neurodamus/stimulus_manager.py

@@ -85,6 +85,7 @@ class BaseStim:
     def __init__(self, _target, stim_info: dict, _cell_manager):
         self.duration = float(stim_info["Duration"])  # duration [ms]
         self.delay = float(stim_info["Delay"])        # start time [ms]
+        self.represents_physical_electrode = stim_info.get('RepresentsPhysicalElectrode', False)
 
 
 @StimulusManager.register_type
@@ -123,7 +124,10 @@ def __init__(self, target, stim_info: dict, cell_manager):
 
                 rng = random.Random123(seed1, seed2, seed3(gid))  # setup RNG
                 ou_args = (self.tau, self.sigma, self.mean, self.duration)
-                ou_kwargs = {'dt': self.dt, 'delay': self.delay, 'rng': rng}
+                ou_kwargs = {
+                    'dt': self.dt, 'delay': self.delay, 'rng': rng,
+                    'physical_electrode': self.represents_physical_electrode
+                }
                 # inject Ornstein-Uhlenbeck signal
                 if stim_info["Mode"] == "Conductance":
                     cs = ConductanceSource.ornstein_uhlenbeck(*ou_args, **ou_kwargs,
@@ -252,7 +256,10 @@ def __init__(self, target, stim_info: dict, cell_manager):
                 rng = random.Random123(seed1, seed2, seed3(gid))  # setup RNG
                 shotnoise_args = (self.tau_D, self.tau_R, self.rate,
                                   self.amp_mean, self.amp_var, self.duration)
-                shotnoise_kwargs = {'dt': self.dt, 'delay': self.delay, 'rng': rng}
+                shotnoise_kwargs = {
+                    'dt': self.dt, 'delay': self.delay, 'rng': rng,
+                    'physical_electrode': self.represents_physical_electrode
+                }
                 # generate shot noise current source
                 if stim_info["Mode"] == "Conductance":
                     cs = ConductanceSource.shot_noise(*shotnoise_args, **shotnoise_kwargs,
@@ -454,8 +461,11 @@ def __init__(self, target, stim_info: dict, cell_manager):
                     continue
 
                 # generate ramp current source
-                cs = CurrentSource.ramp(self.amp_start, self.amp_end, self.duration,
-                                        delay=self.delay)
+                cs = CurrentSource.ramp(
+                    self.amp_start, self.amp_end, self.duration,
+                    delay=self.delay,
+                    physical_electrode=self.represents_physical_electrode
+                )
                 # attach current source to section
                 cs.attach_to(sc.sec, tpoint_list.x[sec_id])
                 self.stimList.append(cs)  # save CurrentSource
@@ -580,8 +590,10 @@ def __init__(self, target, stim_info: dict, cell_manager):
                     continue
 
                 # generate noise current source
-                cs = CurrentSource.noise(self.mean, self.var, self.duration,
-                                         dt=self.dt, delay=self.delay, rng=rng)
+                cs = CurrentSource.noise(
+                    self.mean, self.var, self.duration, dt=self.dt, delay=self.delay, rng=rng,
+                    physical_electrode=self.represents_physical_electrode
+                )
                 # attach current source to section
                 cs.attach_to(sc.sec, tpoint_list.x[sec_id])
                 self.stimList.append(cs)  # save CurrentSource
@@ -660,8 +672,11 @@ def __init__(self, target, stim_info: dict, cell_manager):
                     continue
 
                 # generate pulse train current source
-                cs = CurrentSource.train(self.amp, self.freq, self.width,
-                                         self.duration, delay=self.delay)
+                cs = CurrentSource.train(
+                    self.amp, self.freq, self.width, self.duration,
+                    delay=self.delay,
+                    physical_electrode=self.represents_physical_electrode
+                )
                 # attach current source to section
                 cs.attach_to(sc.sec, tpoint_list.x[sec_id])
                 self.stimList.append(cs)  # save CurrentSource
@@ -696,8 +711,10 @@ def __init__(self, target, stim_info: dict, cell_manager):
                     continue
 
                 # generate sinusoidal current source
-                cs = CurrentSource.sin(self.amp, self.duration, self.freq,
-                                       step=self.dt, delay=self.delay)
+                cs = CurrentSource.sin(
+                    self.amp, self.duration, self.freq, step=self.dt, delay=self.delay,
+                    physical_electrode=self.represents_physical_electrode
+                )
                 # attach current source to section
                 cs.attach_to(sc.sec, tpoint_list.x[sec_id])
                 self.stimList.append(cs)  # save CurrentSource
@@ -735,6 +752,7 @@ def __init__(self, target, stim_info: dict, cell_manager):
 
                 # create single electrode voltage clamp at location
                 seclamp = Nd.h.SEClamp(tpoint_list.x[sec_id], sec=sc.sec)
+
                 seclamp.rs = self.rs
                 seclamp.dur1 = self.duration
                 seclamp.amp1 = self.vhold

tests/integration-e2e/test_dry_run_workflow.py

@@ -24,7 +24,7 @@ def test_dry_run_workflow(USECASE3):
 
     assert 20.0 <= nd._dry_run_stats.cell_memory_total <= 30.0
     assert 0.0 <= nd._dry_run_stats.synapse_memory_total <= 1.0
-    assert 80.0 <= nd._dry_run_stats.base_memory <= 120.0
+    assert 70.0 <= nd._dry_run_stats.base_memory <= 120.0
     expected_items = {
         'L4_PC-dSTUT': 2,
         'L4_MC-dSTUT': 1,

tests/scientific/test_current_injection.py

@@ -0,0 +1,109 @@
+import json
+import numpy
+import os
+import pytest
+from tempfile import NamedTemporaryFile
+
+
+@pytest.fixture
+def sonata_config_files(sonata_config, input_type):
+    config_files = []
+    for represents_physical_electrode in [True, False]:
+        # Create a deep copy of sonata_config for each configuration to avoid conflicts
+        config_copy = json.loads(json.dumps(sonata_config))
+
+        stimulus_config = {
+            "input_type": input_type,
+            "delay": 5,
+            "duration": 2100,
+            "node_set": "l4pc",
+            "represents_physical_electrode": represents_physical_electrode
+        }
+
+        if input_type == "current_clamp":
+            stimulus_config.update({
+                "module": "noise",
+                "mean": 0.05,
+                "variance": 0.01
+            })
+        elif input_type == "conductance":
+            stimulus_config.update({
+                "module": "ornstein_uhlenbeck",
+                "mean": 0.05,
+                "sigma": 0.01,
+                "tau": 0.1
+            })
+
+        config_copy["inputs"] = {"Stimulus": stimulus_config}
+        config_copy["reports"] = {
+            "current": {
+                "type": "summation",
+                "cells": "l4pc",
+                "variable_name": "i_membrane",
+                "unit": "nA",
+                "dt": 0.1,
+                "start_time": 0.0,
+                "end_time": 50.0
+            },
+            "voltage": {
+                "type": "compartment",
+                "cells": "l4pc",
+                "variable_name": "v",
+                "unit": "mV",
+                "dt": 0.1,
+                "start_time": 0.0,
+                "end_time": 50.0
+            }
+        }
+
+        with NamedTemporaryFile("w", suffix='.json', delete=False) as config_file:
+            json.dump(config_copy, config_file)
+            config_files.append(config_file.name)
+
+    yield tuple(config_files)
+
+    # Cleanup
+    for config_file in config_files:
+        os.unlink(config_file)
+
+
+def _read_sonata_soma_report(report_name):
+    import libsonata
+    report = libsonata.SomaReportReader(report_name)
+    pop_name = report.get_population_names()[0]
+    ids = report[pop_name].get_node_ids()
+    data = report[pop_name].get(node_ids=[ids[0]])
+    return numpy.array(data.data).flatten()
+
+
+def _run_simulation(config_file):
+    import subprocess
+    output_dir = "output_current_conductance"
+    command = [
+        "neurodamus",
+        config_file,
+        f"--output-path={output_dir}"
+    ]
+    config_dir = os.path.dirname(config_file)
+    subprocess.run(command, cwd=config_dir, check=True)
+    soma_report_path = os.path.join(config_dir, output_dir, "voltage.h5")
+    return _read_sonata_soma_report(soma_report_path)
+
+
+@pytest.mark.parametrize("input_type", [
+    "current_clamp",
+    "conductance",
+])
+def test_current_conductance_injection(sonata_config_files):
+    """
+    Test the consistency of voltage traces between original and new configurations
+    (set by 'represents_physical_electrode': true/false)
+    under different types of input (current clamp and conductance).
+    """
+    import numpy.testing as npt
+    config_file_original, config_file_new = sonata_config_files
+
+    voltage_vec_original = _run_simulation(config_file_original)
+    voltage_vec_new = _run_simulation(config_file_new)
+
+    npt.assert_equal(voltage_vec_original, voltage_vec_new)

tests/simulations/v5_sonata/simulation_config.json

@@ -92,7 +92,8 @@
             "variance": 0.001,
             "delay": 0.0,
             "duration": 30000,
-            "node_set": "Excitatory"
+            "node_set": "Excitatory",
+            "represents_physical_electrode":true
         },
         "ThresholdInh": {
             "module": "noise",
@@ -101,14 +102,16 @@
             "variance": 0.001,
             "delay": 0.0,
             "duration": 30000,
-            "node_set": "Inhibitory"
+            "node_set": "Inhibitory",
+            "represents_physical_electrode":true
         },
         "hypamp_mosaic": {
             "module": "hyperpolarizing",
             "input_type": "current_clamp",
             "delay": 0.0,
             "duration": 30000,
-            "node_set": "Mosaic"
+            "node_set": "Mosaic",
+            "represents_physical_electrode":true
         }
     },
     "reports": {