diff --git a/bluecellulab/cell/injector.py b/bluecellulab/cell/injector.py index d0384113..c425ab82 100644 --- a/bluecellulab/cell/injector.py +++ b/bluecellulab/cell/injector.py @@ -360,17 +360,15 @@ def add_replay_relative_shotnoise( stimulus: RelativeShotNoise, shotnoise_stim_count=0): """Add a replay relative shot noise stimulus.""" - cv_square = stimulus.amp_cv**2 stim_mode = stimulus.mode rel_prop = self.relativity_proportion(stim_mode) mean = stimulus.mean_percent / 100 * rel_prop sd = stimulus.sd_percent / 100 * rel_prop - var = sd * sd rate, amp_mean, amp_var = get_relative_shotnoise_params( - mean, var, stimulus.decay_time, stimulus.rise_time, cv_square) + mean, sd, stimulus.decay_time, stimulus.rise_time, stimulus.relative_skew) rng = self._get_shotnoise_step_rand(shotnoise_stim_count, stimulus.seed) tvec, svec = gen_shotnoise_signal(stimulus.decay_time, stimulus.rise_time, rate, amp_mean, diff --git a/bluecellulab/cell/stimuli_generator.py b/bluecellulab/cell/stimuli_generator.py index ba556060..e4badaf2 100644 --- a/bluecellulab/cell/stimuli_generator.py +++ b/bluecellulab/cell/stimuli_generator.py @@ -114,17 +114,35 @@ def gen_shotnoise_signal(tau_D, tau_R, rate, amp_mean, amp_var, return tvec, P -def get_relative_shotnoise_params(mean, var, tau_D, tau_R, cv_square): +def get_relative_shotnoise_params(mean, sd, tau_D, tau_R, relative_skew): """Returns Rate, amp_mean and amp_var parameters.""" # bi-exponential time to peak [ms] t_peak = math.log(tau_D / tau_R) / (1 / tau_R - 1 / tau_D) # bi-exponential peak height [1] - x_peak = math.exp(-t_peak / tau_D) - math.exp(-t_peak / tau_R) - - rate_ms = (1 + cv_square) / 2 * (mean ** 2 / var) / (tau_D + tau_R) - rate = rate_ms * 1000 # rate in 1 / s [Hz] - amp_mean = mean * x_peak / rate_ms / (tau_D - tau_R) - amp_var = cv_square * amp_mean ** 2 + F_peak = math.exp(-t_peak / tau_D) - math.exp(-t_peak / tau_R) + + # utility constants + Xi = (tau_D - tau_R) / F_peak + A = 1 / (tau_D + tau_R) + B = 1 / ((tau_D + 2 * tau_R) * (2 * tau_D + tau_R)) + + # skewness + skew_bnd_min = (8 / 3) * (B / A ** 2) * (sd / mean) + skew = (1 + relative_skew) * skew_bnd_min + if skew < skew_bnd_min or skew > 2 * skew_bnd_min: + raise ValueError("skewness out of bounds") + + # cumulants + lambda2_1 = sd ** 2 / mean # lambda2 over lambda1 + lambda3_2 = sd * skew # lambda3 over lambda2 + theta1pk = 2 / (A * Xi) * lambda2_1 # = (1 + k) * theta + theta2pk = (3 * A) / (4 * B * Xi) * lambda3_2 # = (2 + k) * theta + + # derived parameters + amp_mean = 2 * theta1pk - theta2pk # mean amplitude [nA or uS] + amp_var = amp_mean * (theta2pk - theta1pk) # variance of amplitude [nA^2 or uS^2] + rate_ms = mean / (amp_mean * Xi) # event rate in 1 / ms + rate = rate_ms * 1000 # event rate in 1 / s [Hz] return rate, amp_mean, amp_var diff --git a/bluecellulab/stimulus/circuit_stimulus_definitions.py b/bluecellulab/stimulus/circuit_stimulus_definitions.py index 42ed6bab..badde88a 100644 --- a/bluecellulab/stimulus/circuit_stimulus_definitions.py +++ b/bluecellulab/stimulus/circuit_stimulus_definitions.py @@ -170,7 +170,7 @@ def from_blueconfig(cls, stimulus_entry: dict) -> Optional[Stimulus]: decay_time=stimulus_entry["DecayTime"], mean_percent=stimulus_entry["MeanPercent"], sd_percent=stimulus_entry["SDPercent"], - amp_cv=stimulus_entry["AmpCV"], + relative_skew=stimulus_entry.get("RelativeSkew", 0.5), seed=stimulus_entry.get("Seed", None), mode=mode, reversal=stimulus_entry.get("Reversal", 0.0) @@ -269,7 +269,7 @@ def from_sonata(cls, stimulus_entry: dict) -> Optional[Stimulus]: decay_time=stimulus_entry["decay_time"], mean_percent=stimulus_entry["mean_percent"], sd_percent=stimulus_entry["sd_percent"], - amp_cv=stimulus_entry["amp_cv"], + relative_skew=stimulus_entry.get("RelativeSkew", 0.5), seed=stimulus_entry.get("random_seed", None), mode=ClampMode(stimulus_entry.get("input_type", "current_clamp").lower()), reversal=stimulus_entry.get("reversal", 0.0) @@ -365,7 +365,7 @@ class RelativeShotNoise(Stimulus): decay_time: float mean_percent: float sd_percent: float - amp_cv: float + relative_skew: float = 0.5 dt: float = 0.25 seed: Optional[int] = None mode: ClampMode = ClampMode.CURRENT @@ -378,6 +378,13 @@ def decay_time_gt_rise_time(cls, v, values): raise ValueError("decay_time must be greater than rise_time") return v + @field_validator("relative_skew") + @classmethod + def relative_skew_in_range(cls, v): + if v < 0.0 or v > 1.0: + raise ValueError("relative skewness must be in [0,1]") + return v + @dataclass(frozen=True, config=dict(extra="forbid")) class OrnsteinUhlenbeck(Stimulus): diff --git a/tests/test_cell/test_injector.py b/tests/test_cell/test_injector.py index fa2ff866..b01c454e 100644 --- a/tests/test_cell/test_injector.py +++ b/tests/test_cell/test_injector.py @@ -351,26 +351,34 @@ def test_add_replay_relative_shotnoise(self): """Unit test for add_replay_relative_shotnoise.""" stimulus = RelativeShotNoise( target="single-cell", delay=0, duration=2, - rise_time=0.4, decay_time=4, mean_percent=70, sd_percent=40, amp_cv=0.63, + rise_time=0.4, decay_time=4, mean_percent=70, sd_percent=40, seed=12, ) self.cell.threshold = 0.184062 soma = self.cell.soma segx = 0.5 tvec, svec = self.cell.add_replay_relative_shotnoise(soma, segx, stimulus) - assert svec.to_python() == approx([0., 0., 0., 0., 0.0204470197, 0.0301526984, - 0.0341840080, 0.0352485557, 0.0347913472, 0.]) + assert svec.to_python() == approx([0., 0., 0., 0.0261496222807232, 0.07428180924441198, + 0.09639263342799839, 0.10479651831172629, 0.10607149236338352, + 0.10381817238692868, 0.0]) assert tvec.to_python() == approx([0.0, 0.25, 0.5, 0.75, 1.0, 1.25, 1.5, 1.75, 2.0, 2.0]) with raises(ValidationError): invalid_stim = RelativeShotNoise( target="single-cell", delay=0, duration=2, - rise_time=4, decay_time=4, mean_percent=70, sd_percent=40, amp_cv=0.63, + rise_time=4, decay_time=4, mean_percent=70, sd_percent=40, seed=12, ) self.cell.add_replay_relative_shotnoise(soma, segx, invalid_stim) + with raises(ValidationError): + invalid_rel_skew = RelativeShotNoise( + target="single-cell", delay=0, duration=2, + rise_time=0.4, decay_time=4, mean_percent=70, sd_percent=40, + seed=12, relative_skew=1.3 + ) + with raises(ZeroDivisionError): self.cell.threshold = 0.0 self.cell.add_replay_relative_shotnoise(soma, segx, stimulus) diff --git a/tests/test_cell/test_stimuli_generator.py b/tests/test_cell/test_stimuli_generator.py index 316b2119..0d6fa270 100644 --- a/tests/test_cell/test_stimuli_generator.py +++ b/tests/test_cell/test_stimuli_generator.py @@ -28,11 +28,11 @@ def test_gen_shotnoise_signal(): def test_get_relative_shotnoise_params(): """Unit test for _get_relative_shotnoise_params.""" rate, amp_mean, amp_var = get_relative_shotnoise_params( - mean=40e-3, var=16e-4, tau_D=4.0, tau_R=0.4, cv_square=0.63**2 + mean=40e-3, sd=0.04, tau_D=4.0, tau_R=0.4, relative_skew=0.5 ) - assert rate == approx(158.73863636363635) - assert amp_mean == approx(0.048776006926722876) - assert amp_var == approx(0.0009442643342459686) + assert rate == approx(227.27272727272705) + assert amp_mean == approx(0.03406760203796963) + assert amp_var == approx(0.0011606015086174707) def test_gen_ornstein_uhlenbeck(): diff --git a/tox.ini b/tox.ini index 2309407a..b279bcea 100644 --- a/tox.ini +++ b/tox.ini @@ -44,7 +44,7 @@ commands = ruff check . --select F541,F401 --per-file-ignores="__init__.py:F401" pycodestyle {[base]name} --ignore=E501,W504,W503 pycodestyle tests --ignore=E501,W504,W503,E741 - mypy . --ignore-missing-imports --disable-error-code=call-overload # remove once pandas-stubs makes a release after 1.5.3.230321 + mypy . --ignore-missing-imports --disable-error-code=call-overload # remove once pandas-stubs makes a release after 1.5.3.230321 docformatter --check bluecellulab -r [testenv:examples]