replace amp_cv with relative_skew

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,

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
 

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):

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)

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():

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]