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