[BBPBGLIB-1137] Account for nonzero total current when using SEClamp…

…; allow nonzero total current when simulating physical electrodes (#134)

## Context
NEURON reports SEClamp currents at different time points than all other
currents, thus giving inaccurate results in lfp calculations, in cases
where the SEClamp represents synaptic inputs.

This can be fixed by replacing the SEClamp with a new point process
mechanism that is identical, except for producing a NONSPECIFIC_CURRENT
instead of an ELECTRODE_CURRENT.

However, in cases where SEClamp or IClamp are used to represent physical
electrodes, rather than synaptic inputs, the currents should not be
included in the LFP calcualtion

## Scope
When a current_source or conductance_source is used, we check if the
alternative mod file, which is expected to be named
'MembraneCurrentSource.mod' or 'ConductanceSource.mod' is available. If
so, it is used, with nothing else changed, unless the key-value pair
"represents_physical_electrode" is set to True in the Stimulus Config
file. Otherwise the IClamp or SEClamp is used as before.

## Testing
Scientific test added to ensure that new mechanisms produce identical membrane voltages to the old ones

## Review
* [x] PR description is complete
* [x] Coding style (imports, function length, New functions, classes or
files) are good
* [x] Unit/Scientific test added
* [x] Updated Readme, in-code, developer documentation

---------

Co-authored-by: Jorge Blanco Alonso <jorge.blancoalonso@epfl.ch>
Co-authored-by: Jorge Blanco Alonso <41900536+jorblancoa@users.noreply.github.com>

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": {