everything for 1ms kernels

dataoverview.py

@@ -6,9 +6,6 @@ from scipy.stats import mannwhitneyu
 from thunderlab.tabledata import TableData
 from plotstyle import plot_style, lighter, significance_str
 
-
-data_path = Path('data')
-
 from noisesplit import model_cell as model_split_example
 from modelsusceptcontrasts import model_cells as model_contrast_examples
 from modelsusceptlown import model_cell as model_lown_example
@@ -18,6 +15,9 @@ from noisesplit import example_cell as punit_split_example
 from ampullaryexamplecell import example_cell as ampul_example
 from ampullaryexamplecell import example_cells as ampul_examples
 
+
+data_path = Path('data')
+
 model_examples = ([[model_lown_example, 0.01],
                    [model_lown_example, 0.03],
                    [model_lown_example, 0.1]],
@@ -26,6 +26,8 @@ model_examples = ([[model_lown_example, 0.01],
 punit_examples = (punit_example, [punit_split_example], punit_examples)
 ampul_examples = (ampul_example, [], ampul_examples)
 
+respmod = 'respmod1'
+
 
 def plot_corr(ax, data, xcol, ycol, zcol, zmin, zmax, xpdfmax, cmap, color,
               si_thresh, example=[], split_example=[], examples=[]):
@@ -266,7 +268,7 @@ def si_stats(title, data, sicol, si_thresh, nsegscol):
     print(f'    total duration:     {np.min(duration):.1f}s - {np.max(duration):.1f}s, median={np.median(duration):.1f}s, mean={np.mean(duration):.1f}s, std={np.std(duration):.1f}s')
     duration = data['durationbase']
     print(f'    baseline duration:  {np.min(duration):.1f}s - {np.max(duration):.1f}s, median={np.median(duration):.1f}s, mean={np.mean(duration):.1f}s, std={np.std(duration):.1f}s')
-    cols = ['cvbase', 'respmod2', 'ratebase', 'vsbase', 'serialcorr1', 'burstfrac', 'ratestim', 'cvstim']
+    cols = ['cvbase', respmod, 'ratebase', 'vsbase', 'serialcorr1', 'burstfrac', 'ratestim', 'cvstim']
     for i in range(len(cols)):
         for j in range(i + 1, len(cols)):
             xcol = cols[i]
@@ -285,7 +287,7 @@ def plot_cvbase_si_punit(ax, data, ycol, si_thresh, color):
     ax.set_ylim(0, 6.5)
     ax.set_yticks_delta(2)
     examples = punit_examples if 'stimindex' in data else model_examples
-    #cax = plot_corr(ax, data, 'cvbase', ycol, 'respmod2', 0, 250, 3,
+    #cax = plot_corr(ax, data, 'cvbase', ycol, respmod, 0, 250, 3,
     #                'coolwarm', color, si_thresh, *examples)
     #cax.set_ylabel('Response mod.', 'Hz')
     plot_corr_contrast(ax, data, 'cvbase', ycol, 3, color,
@@ -300,7 +302,7 @@ def plot_cvstim_si_punit(ax, data, ycol, si_thresh, color):
     ax.set_ylim(0, 6.5)
     ax.set_yticks_delta(2)
     examples = punit_examples if 'stimindex' in data else model_examples
-    #cax = plot_corr(ax, data, 'cvstim', ycol, 'respmod2', 0, 250, 2,
+    #cax = plot_corr(ax, data, 'cvstim', ycol, respmod, 0, 250, 2,
     #                'coolwarm', color, si_thresh, *examples)
     #cax.set_ylabel('Response mod.', 'Hz')
     #cax = plot_corr(ax, data, 'cvstim', ycol, 'cvbase', 0, 1.5, 2,
@@ -318,16 +320,16 @@ def plot_cvstim_si_punit(ax, data, ycol, si_thresh, color):
     
 def plot_rmod_si_punit(ax, data, ycol, si_thresh, color):
     ax.set_xlabel('Response modulation', 'Hz')
-    ax.set_xlim(0, 250)
+    ax.set_xlim(0, 350)
     ax.set_xticks_delta(100)
     ax.set_ylabel('SI($r$)')
     ax.set_ylim(0, 6.5)
     ax.set_yticks_delta(2)
     examples = punit_examples if 'stimindex' in data else model_examples
-    #cax = plot_corr(ax, data, 'respmod2', ycol, 'cvbase', 0, 1.5, 0.016,
+    #cax = plot_corr(ax, data, respmod, ycol, 'cvbase', 0, 1.5, 0.016,
     #                'coolwarm', color, si_thresh, *examples)
     #cax.set_ylabel('CV$_{\\rm base}$')
-    plot_corr_contrast(ax, data, 'respmod2', ycol, 0.016, color,
+    plot_corr_contrast(ax, data, respmod, ycol, 0.016, color,
                        si_thresh, *examples)
 
     
@@ -355,7 +357,7 @@ def plot_cvbase_si_ampul(ax, data, ycol, si_thresh, color):
     ax.set_ylabel('SI($r$)')
     ax.set_ylim(0, 10)
     ax.set_yticks_delta(2)
-    #cax = plot_corr(ax, data, 'cvbase', ycol, 'respmod2', 0, 80, 20,
+    #cax = plot_corr(ax, data, 'cvbase', ycol, respmod, 0, 80, 20,
     #                'coolwarm', color, si_thresh, *ampul_examples)
     #cax.set_ylabel('Response mod.', 'Hz')
     plot_corr_contrast(ax, data, 'cvbase', ycol, 20, color,
@@ -369,7 +371,7 @@ def plot_cvstim_si_ampul(ax, data, ycol, si_thresh, color):
     ax.set_ylabel('SI($r$)')
     ax.set_ylim(0, 10)
     ax.set_yticks_delta(2)
-    #cax = plot_corr(ax, data, 'cvstim', ycol, 'respmod2', 0, 80, 6,
+    #cax = plot_corr(ax, data, 'cvstim', ycol, respmod, 0, 80, 6,
     #                'coolwarm', color, si_thresh, *ampul_examples)
     #cax.set_ylabel('Response mod.', 'Hz')
     #cax = plot_corr(ax, data, 'cvstim', ycol, 'cvbase', 0, 0.2, 6,
@@ -386,16 +388,16 @@ def plot_cvstim_si_ampul(ax, data, ycol, si_thresh, color):
     
 def plot_rmod_si_ampul(ax, data, ycol, si_thresh, color):
     ax.set_xlabel('Response modulation', 'Hz')
-    ax.set_xlim(0, 80)
-    ax.set_xticks_delta(20)
+    ax.set_xlim(0, 100)
+    ax.set_xticks_delta(30)
     ax.set_ylabel('SI($r$)')
     ax.set_ylim(0, 10)
     ax.set_yticks_delta(2)
-    #cax = plot_corr(ax, data, 'respmod2', ycol, 'cvbase', 0, 0.2, 0.06,
+    #cax = plot_corr(ax, data, respmod, ycol, 'cvbase', 0, 0.2, 0.06,
     #                'coolwarm', color, si_thresh, *ampul_examples)
     #cax.set_ylabel('CV$_{\\rm base}$')
     #cax.set_yticks_delta(0.1)
-    plot_corr_contrast(ax, data, 'respmod2', ycol, 0.06, color,
+    plot_corr_contrast(ax, data, respmod, ycol, 0.06, color,
                        si_thresh, *ampul_examples)
 
         
@@ -412,7 +414,7 @@ def plot_rate_si_ampul(ax, data, ycol, si_thresh, color):
     #cax.set_yticks_delta(0.1)
     plot_corr_contrast(ax, data, 'ratebase', ycol, 0.06, color,
                        si_thresh, *ampul_examples)
-    ax.legend(title='contrast', loc='upper right', bbox_to_anchor=(1.95, 1),
+    ax.legend(title='contrast', loc='upper right', bbox_to_anchor=(1.9, 1),
               markerfirst=False, handletextpad=0.5)
 
     
@@ -453,8 +455,8 @@ if __name__ == '__main__':
     print(f'    baseline rate model: min={np.min(ratemodel):3.0f}Hz max={np.max(ratemodel):3.0f}Hz median={np.median(ratemodel):3.0f}Hz')
     print(f'    baseline rate data:  min={np.min(ratedata):3.0f}Hz max={np.max(ratedata):3.0f}Hz median={np.median(ratedata):3.0f}Hz')
     print()
-    rmmodel = punit_model['respmod2']
-    rmdata = punit_data['respmod2']
+    rmmodel = punit_model[respmod]
+    rmdata = punit_data[respmod]
     u, p = mannwhitneyu(rmmodel, rmdata)
     print('Response modulation differs between P-unit models and data:')
     print(f'    U={u:g}, p={p:.2g}')
@@ -489,7 +491,7 @@ if __name__ == '__main__':
     s = plot_style()
     fig, axs = plt.subplots(3, 3, cmsize=(s.plot_width, 0.75*s.plot_width),
                             height_ratios=[1, 0, 1, 0.3, 1])
-    fig.subplots_adjust(leftm=6.5, rightm=18, topm=4.5, bottomm=4,
+    fig.subplots_adjust(leftm=6.5, rightm=17.5, topm=4.5, bottomm=4,
                         wspace=0.6, hspace=0.6)
     
     si_stats('P-unit model:', punit_model, 'dsinorm100', si_thresh,

nonlinearbaseline.tex

@@ -278,7 +278,7 @@ Supported by SPP 2205 ``Evolutionary optimisation of neuronal processing'' by th
 % 250 words
 \section{Abstract}
 Spiking thresholds in neurons or rectification at synapses are essential for neuronal computations rendering neuronal processing inherently nonlinear. Nevertheless, linear response theory has been instrumental for understanding, for example, the impact of noise or neuronal synchrony on signal transmission, or the emergence of oscillatory activity, but is valid only at low stimulus amplitudes or large levels of intrinsic noise. At higher signal-to-noise ratios, however, nonlinear response components become relevant. Theoretical results for leaky integrate-and-fire neurons in the weakly nonlinear regime suggest strong responses at the sum of two input frequencies if one of these frequencies or their sum match the neuron's baseline firing rate.
-We here analyze nonlinear responses in two types of primary electroreceptor afferents, the P-units of the active and the ampullary cells of the passive electrosensory system of the wave-type electric fish \textit{Apteronotus leptorhynchus}. In our combined experimental and modeling approach we identify these predicted nonlinear responses primarily in low-noise P-units and in more than every second ampullary cell. Our results provide experimental evidence for nonlinear responses of spike generators in the weakly nonlinear regime. We conclude that such nonlinear responses occur in any sensory neuron that operates in similar regimes particularly at near-threshold stimulus conditions.
+We here analyze nonlinear responses in two types of primary electroreceptor afferents, the P-units of the active and the ampullary cells of the passive electrosensory system of the wave-type electric fish \textit{Apteronotus leptorhynchus}. In our combined experimental and modeling approach we identify these predicted nonlinear responses primarily in a few low-noise P-units and in more than every second ampullary cell. Our results provide experimental evidence for nonlinear responses of spike generators in the weakly nonlinear regime. We conclude that such nonlinear responses occur in any sensory neuron that operates in similar regimes particularly at near-threshold stimulus conditions.
 
 % max 120 words
 \section{Significance statement}
@@ -333,7 +333,7 @@ The P-unit model parameters and spectral analysis algorithms are available at \u
 The baseline firing rate $r$ was calculated as the number of spikes divided by the duration of the baseline recording (median 32\,s). The coefficient of variation (CV) of the interspike intervals (ISI) is their standard deviation relative to their mean: $\rm{CV}_{\rm base} = \sqrt{\langle (ISI- \langle ISI \rangle) ^2 \rangle} / \langle ISI \rangle$. If the baseline was recorded several times in a recording, the measures from the longest recording were taken.
 
 \paragraph{White noise analysis} \label{response_modulation}
-When stimulated with band-limited white noise stimuli, neuronal activity is modulated around the average firing rate that is similar to the baseline firing rate and in that way encodes the time-course of the stimulus. For an estimate of the time-dependent firing rate $r(t)$ we convolved each spike train with normalized Gaussian kernels with standard deviation of 2\,ms and averaged the resulting single-trail firing rates over trials. The response modulation quantifies the variation of $r(t)$ computed as the standard deviation in time $\sigma_{s} = \sqrt{\langle (r(t)-\langle r(t) \rangle_t )^2\rangle_t}$, where $\langle \cdot \rangle_t$ denotes averaging over time.
+When stimulated with band-limited white noise stimuli, neuronal activity is modulated around the average firing rate that is similar to the baseline firing rate and in that way encodes the time-course of the stimulus. For an estimate of the time-dependent firing rate $r(t)$ we convolved each spike train with normalized Gaussian kernels with standard deviation of 1\,ms and averaged the resulting single-trail firing rates over trials. The response modulation quantifies the variation of $r(t)$ computed as the standard deviation in time $\sigma_{s} = \sqrt{\langle (r(t)-\langle r(t) \rangle_t )^2\rangle_t}$, where $\langle \cdot \rangle_t$ denotes averaging over time.
 
 \paragraph{Spectral analysis}\label{susceptibility_methods}
 To characterize the relation between the spiking response evoked by white-noise  stimuli, we estimated the first- and second-order susceptibilities in the frequency domain. For this we converted spike times into binary vectors $x_k$ with $\Delta t = 0.5$\,ms wide bins that are set to 2\,kHz where a spike occurred and zero otherwise. Fast Fourier transforms (FFT) $S(\omega)$ and $X(\omega)$ of the stimulus $s_k$ (also down-sampled to a sampling rate of 2\,kHz) and $x_k$, respectively, were computed numerically according to
@@ -650,7 +650,7 @@ In the P-unit models, each model cell contributed with three RAM stimulus presen
 
 The effective stimulus strength also plays a role in predicting the SI($r$) values. We quantify the effect of stimulus strength on a cell's response by the response modulation --- the standard deviation of a cell's firing rate in response to a RAM stimulus. The lower the response modulation, i.e. the weaker the effective stimulus, the higher the S($r$) (\figrefb{fig:dataoverview}\,\panel[ii]{A}). Although there is a tendency of low stimulus contrasts to evoke lower response modulations, response modulations evoked by each of the three contrasts overlap substantially, emphasizing the strong heterogeneity of the P-units' sensitivity \citep{Grewe2017}. Cells with high SI($r$) values are the ones with baseline firing rate below 200\,Hz (\figrefb{fig:dataoverview}\,\panel[iii]{A}).
 
-In comparison to the experimentally measured P-unit recordings, the model cells are skewed to lower baseline CVs (Mann-Whitney $U=13986$, $p=3\times 10^{-9}$), because the models are not able to reproduce bursting, which we observe in many P-units and which leads to high CVs. Also the response modulation of the models is skewed to lower values (Mann-Whitney $U=15312$, $p=7\times 10^{-7}$), because in the measured cells, response modulation is positively correlated with baseline CV (Pearson $R=0.45$, $p=1\times 10^{-19}$), i.e. bursting cells are more sensitive. Median baseline firing rate in the models is by 53\,Hz smaller than in the experimental data (Mann-Whitney $U=17034$, $p=0.0002$).
+In comparison to the experimentally measured P-unit recordings, the model cells are skewed to lower baseline CVs (Mann-Whitney $U=13986$, $p=3\times 10^{-9}$), because the models are not able to reproduce bursting, which we observe in many P-units and which leads to high CVs. Also the response modulation of the models is skewed to lower values (Mann-Whitney $U=14051$, $p=4\times 10^{-9}$), because in the measured cells, response modulation is positively correlated with baseline CV (Pearson $R=0.34$, $p=1\times 10^{-4}$), i.e. bursting cells are more sensitive. Median baseline firing rate in the models is by 53\,Hz smaller than in the experimental data (Mann-Whitney $U=17034$, $p=0.0002$).
 
 In the experimentally measured P-units, each of the $172$ cells contributes on average with two RAM stimulus presentations, presented at contrasts ranging from 1 to 20\,\% to the 376 samples. Despite the mentioned differences between the P-unit models and the measured data, the SI($r$) values do not differ between models and data (median of 1.3, Mann-Whitney $U=19702$, $p=0.09$) and also 16\,\% of the samples from all presented stimulus contrasts exceed the threshold of 1.8. The SI($r$) values of the P-unit population correlate weakly with the CV of the baseline ISIs that range from 0.18 to 1.35 (median 0.49). Cells with lower baseline CVs tend to have more pronounced ridges in their second-order susceptibilities than those with higher baseline CVs (\figrefb{fig:dataoverview}\,\panel[i]{B}).
 

rebuttal2.tex

@@ -229,13 +229,15 @@
   P-units, but "much stronger" does not clearly convey this,
   especially in the abstract.}
 
-\response{We changed the sentence to ``...  we identify these predicted nonlinear responses primarily in low-noise P-units and in more than every second ampullary cell.''}
+\response{We changed the sentence to ``...  we identify these
+  predicted nonlinear responses primarily in a few low-noise P-units
+  and in more than every second ampullary cell.''}
 
 \issue{(3) Figure 1A. "r" needs to be clearly defined here. Based on
   the text, it seems to be the baseline firing rate of the neuron, but
   this needs to be made clear in the figure legend.}
 
-\response{We added a brief sentence to the caption.}
+\response{We added a brief explanatory sentence to the caption.}
 
 \issue{(4) Figure 1B. "Because frequencies can also be negative..."
   This is unclear and needs more explanation, especially because there
@@ -251,15 +253,15 @@
   clipped in these two figures.}
 
 \response{You are right. In figure 4 we show now the spectrum up to
-  750Hz, such that fEOD and its interactions with df2 and harmonics
-  are included. We labeled the additonal peaks accordingly. In figure
-  3 we stay with the small range, because we have so little data (only
-  three trials of 500ms duration) for this special setting where one
-  of the beat frequencies approximately matches the P-units baseline
-  firing rate. This is why the power spectra are very noisy. Also, for
-  an introductory figure we prefer to only show the few peaks that are
-  relevant for the rest of the manuscript, such that the reader does
-  not get overwhelmed.  }
+  750\,Hz, such that $f_{EOD}§ and its interactions with $\Delta f_2$
+  and harmonics are included. We labeled the additonal peaks
+  accordingly. In figure 3 we stay with the small range, because we
+  have so little data for this special setting where one of the beat
+  frequencies approximately matches the P-units baseline firing rate
+  (only three trials of 500ms duration). This is why the power
+  spectra are very noisy. Also, for an introductory figure we prefer
+  to only show the few peaks that are relevant for the rest of the
+  manuscript, such that the reader does not get overwhelmed.}
 
 \issue{(6) Figure 3. Why are these example firing rates based on
   convolution with a 1 ms Gaussian kernel if the analyses were based
@@ -268,15 +270,34 @@
   actually analyzed. More fundamentally, why would a 2-fold difference
   in kernel width be appropriate for presentation vs. analysis?}
 
-\response{This was for historical reasons. We updated figure 3 to also
-  use the 2ms kernel. Now all firing rates in the manuscript are based
-  on the 2ms kernel.}
+\response{This was for historical reasons. We now decided to use the
+  1\,ms kernel for all figures and analysis. In doing so we also added
+  panels showing firing rates in addition to the response spectra in
+  figure 4. Using the more narrow kernel better reveals the details of
+  the time course of the firing rate and this way improves the
+  connection between the firing rate and the response spectra. In
+  figure 10, middel column, the range of possible values of the
+  response modulations is a bit enlarged by using the 1\,ms kernel,
+  but the correlations and their significance did not change a lot
+  either.}
 
 \issue{(7) Figure 3D legend. The relationship between 2nd order AM
   (envelope) and the two nonlinear peaks should be made clear. I
   believe the envelope is represented by both peaks, correct?}
 
-\response{ADD SENTENCE.}
+\response{No, what is shown is the power spectrum of the spike
+  response, not the one of the amplitude modulation or envelope of the
+  stimulus. We added a sentence to the end of the figure caption to
+  make this clear.}
+
+\response{If it were the power spectrum of the signal after it passed
+  a non-linearity (rectification or threhsolding at zero), then there
+  could be also peaks at the sum and difference of the beat
+  frequencies. However, since they are close to the higher one of the
+  two beat frequencies they do not show up in the AM as obviously as
+  for the settings used in the social envelope papers by Eric Fortune
+  and Andre Longtin and colleges (I guess this is what you have in
+  mind).}
 
 \issue{(8) Line 302. "not-small amplitude" is arbitrary and
   vague. Please be clearer and more precise.}

regimes.py

@@ -16,22 +16,22 @@ amax = 60
 cthresh1 = 1.2
 cthresh2 = 3.5
 
-model_cell = '2018-05-08-ab-invivo-1'      # 116, CV=0.68
-alphas = [0.002, 0.008, 0.025, 0.05]
-rmax = 400
-amax = 50
-cthresh1 = 1.2
-cthresh2 = 3.0
+#model_cell = '2018-05-08-ab-invivo-1'      # 116, CV=0.68
+#alphas = [0.002, 0.008, 0.025, 0.05]
+#rmax = 400
+#amax = 50
+#cthresh1 = 1.2
+#cthresh2 = 3.0
 
 data_path = Path('data')
 sims_path = data_path / 'simulations'
 
 trials = 1000
-spec_trials = 100   # set to zero to only recompute firng rates
+spec_trials = 1000   # set to zero to only recompute firng rates
 sigma = 0.001
 nfft = 2**18
 
-recompute = False
+recompute = True
 
 
 def load_data(file_path):
@@ -268,18 +268,18 @@ def plot_psd(ax, s, path, contrast, spikes, nfft, dt, beatf1, beatf2, eodf):
         ax.plot(beatf2 + 2*beatf1, decibel(peak_ampl(freqs, psd, beatf2 + 2*beatf1)) + offsm,
                 label=r'$\Delta f_2 \pm k\Delta f_1$', zorder=40, **s.psF012m)
         ax.plot(beatf2 + 3*beatf1, decibel(peak_ampl(freqs, psd, beatf2 + 3*beatf1)) + offsm, zorder=40, **s.psF012m)
-        #ax.plot(beatf2 - 2*beatf1, decibel(peak_ampl(freqs, psd, beatf2 - 2*beatf1)) + offsm, zorder=40, **s.psF012m)
-        #ax.plot(beatf2 - 3*beatf1, decibel(peak_ampl(freqs, psd, beatf2 - 3*beatf1)) + offsm, zorder=40, **s.psF012m)
+        ax.plot(beatf2 - 2*beatf1, decibel(peak_ampl(freqs, psd, beatf2 - 2*beatf1)) + offsm, zorder=40, **s.psF012m)
+        ax.plot(beatf2 - 3*beatf1, decibel(peak_ampl(freqs, psd, beatf2 - 3*beatf1)) + offsm, zorder=40, **s.psF012m)
         ax.plot(3*beatf1, decibel(peak_ampl(freqs, psd, 3*beatf1)) + offsm, zorder=40, **s.psF01m)
         ax.plot(4*beatf1, decibel(peak_ampl(freqs, psd, 4*beatf1)) + offsm, zorder=40, **s.psF01m)
         ax.plot(eodf - 2*beatf1, decibel(peak_ampl(freqs, psd, eodf - 2*beatf1)) + offsm, zorder=40, **s.psFEODm)
         ax.plot(eodf - 3*beatf1, decibel(peak_ampl(freqs, psd, eodf - 3*beatf1)) + offsm, zorder=40, **s.psFEODm)
-        #ax.plot(eodf - 4*beatf1, decibel(peak_ampl(freqs, psd, eodf - 4*beatf1)) + offsm, zorder=40, **s.psFEODm)
+        ax.plot(eodf - 4*beatf1, decibel(peak_ampl(freqs, psd, eodf - 4*beatf1)) + offsm, zorder=40, **s.psFEODm)
         ax.plot(eodf - 2*beatf2, decibel(peak_ampl(freqs, psd, eodf - 2*beatf2)) + offsm, zorder=40, **s.psF0m)
-        #ax.plot(eodf - beatf2 + 2*beatf1, decibel(peak_ampl(freqs, psd, eodf - beatf2 + 2*beatf1)) + offsm,
-        #        zorder=40, **s.psF02m)
-        #ax.plot(eodf - beatf2 + 3*beatf1, decibel(peak_ampl(freqs, psd, eodf - beatf2 + 2*beatf1)) + offsm,
-        #        zorder=40, **s.psF02m)
+        ax.plot(eodf - beatf2 + 2*beatf1, decibel(peak_ampl(freqs, psd, eodf - beatf2 + 2*beatf1)) + offsm,
+                zorder=40, **s.psF02m)
+        ax.plot(eodf - beatf2 + 3*beatf1, decibel(peak_ampl(freqs, psd, eodf - beatf2 + 2*beatf1)) + offsm,
+                zorder=40, **s.psF02m)
     ax.set_xlim(0, 750)
     ax.set_ylim(-60, 0)
     ax.set_xticks_delta(200)
@@ -399,8 +399,8 @@ if __name__ == '__main__':
     fig, (axes, axa) = plt.subplots(2, 1, height_ratios=[5, 3],
                                     cmsize=(s.plot_width, 0.7*s.plot_width))
     fig.subplots_adjust(leftm=8, rightm=2, topm=2, bottomm=3.5, hspace=0.6)
-    axe = axes.subplots(4, 4, wspace=0.4, hspace=0.2,
-                        height_ratios=[1, 2, 2, 0.6, 3])
+    axe = axes.subplots(4, 4, wspace=0.4, hspace=0.3,
+                        height_ratios=[1, 2, 2, 0.2, 3])
     fig.show_spines('lb')
     
     # example power spectra:

twobeats.py

@@ -12,7 +12,7 @@ cell = '2021-08-03-ac-invivo-1'
 
 data_path = Path('data')
 
-sigma = 0.002
+sigma = 0.001
 
 
 def load_spikes(cell_path, f1=797, f2=631):