Replaces SEClam and IClamp with new ConductanceSource and MembraneCur…

…rentSource point processes

docs/online-lfp.rst

@@ -140,6 +140,8 @@ 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, it is neccessary that total current over the neuron sums to zero in order to produce valid extracellular recording results. For IClamp electrodes, this is always the case. However, the Neuron SEClamp class does not fulfill this criterion due to numerical issues. We have created a new point process, `new_conductance_source`, available in `neurodamus-neocortex`, which does fulfill the criterion. Therefore, If an SEClamp source is present in the simulation config file, and `new_conductance_source` is compiled, this will be used instead of the SEClamp mechanism. If `new_conductance_source` is not available, SEClamp will be used, and extracellular recording results should not be trusted.
+
 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('represents_physical_electrode', 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