From 1d2bf99009b5d6d0726c8d5f61b94bd1fdd95277 Mon Sep 17 00:00:00 2001 From: joseph-tharayil Date: Fri, 10 May 2024 09:33:25 +0200 Subject: [PATCH] [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 Co-authored-by: Jorge Blanco Alonso <41900536+jorblancoa@users.noreply.github.com> --- docs/online-lfp.rst | 6 + neurodamus/core/stimuli.py | 46 ++++++-- neurodamus/stimulus_manager.py | 38 ++++-- .../integration-e2e/test_dry_run_workflow.py | 2 +- tests/scientific/test_current_injection.py | 109 ++++++++++++++++++ .../v5_sonata/simulation_config.json | 9 +- 6 files changed, 184 insertions(+), 26 deletions(-) create mode 100644 tests/scientific/test_current_injection.py diff --git a/docs/online-lfp.rst b/docs/online-lfp.rst index 9fc554b2..577ecb7f 100644 --- a/docs/online-lfp.rst +++ b/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 diff --git a/neurodamus/core/stimuli.py b/neurodamus/core/stimuli.py index 92debf55..c9bc0118 100644 --- a/neurodamus/core/stimuli.py +++ b/neurodamus/core/stimuli.py @@ -10,12 +10,14 @@ 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() @@ -23,6 +25,7 @@ def __init__(self, base_amp=0.0, *, delay=0, rng=None): 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,11 +370,12 @@ 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) @@ -379,7 +383,14 @@ 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 diff --git a/neurodamus/stimulus_manager.py b/neurodamus/stimulus_manager.py index 169e7a86..2b40b8a3 100644 --- a/neurodamus/stimulus_manager.py +++ b/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 diff --git a/tests/integration-e2e/test_dry_run_workflow.py b/tests/integration-e2e/test_dry_run_workflow.py index c0335feb..fa247318 100644 --- a/tests/integration-e2e/test_dry_run_workflow.py +++ b/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, diff --git a/tests/scientific/test_current_injection.py b/tests/scientific/test_current_injection.py new file mode 100644 index 00000000..3e220387 --- /dev/null +++ b/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) diff --git a/tests/simulations/v5_sonata/simulation_config.json b/tests/simulations/v5_sonata/simulation_config.json index c7efd394..a8d9d516 100644 --- a/tests/simulations/v5_sonata/simulation_config.json +++ b/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": {