updated regimes fiugre

ampullaryexamplecell.py

@@ -232,7 +232,7 @@ if __name__ == '__main__':
     axs[0, 0].text(0, 1.6, 'Ampullary cells:', transform=axs[0, 0].transAxes,
                    color=s.cell_color1,
             fontsize='large')
-    axs[0, -1].text(1, 0, '5\\,ms', ha='right',
+    axs[0, -1].text(0.97, -0.45, '5\\,ms', ha='right',
                     transform=axs[0, -1].transAxes)
     plot_colorbar(axs[1, -1], pc, 0.4)
     fig.common_yticks(axs[1, :])

nonlinearbaseline.tex

@@ -113,6 +113,8 @@
 \newcolumntype{C}[1]{>{\centering\let\newline\\\arraybackslash\hspace{0pt}}m{#1}}
 \newcolumntype{R}[1]{>{\raggedleft\let\newline\\\arraybackslash\hspace{0pt}}m{#1}}
 
+\graphicspath{{../susceptibility1}}
+
 %%%%% figures %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 % captions:
@@ -431,7 +433,7 @@ Here we search for such weakly nonlinear responses in electroreceptors of the tw
 \section{Results}
 
 \begin{figure*}[t]
-  %\includegraphics[width=\columnwidth]{motivation.pdf}
+  \includegraphics[width=\columnwidth]{motivation.pdf}
   \caption{\label{fig:motivation} Nonlinearity in an electrophysiologically recorded P-unit of \lepto{} in a three-fish setting (cell identifier ``2021-08-03-ac"). Receiver with EOD frequency $\feod{} =664$\,Hz encounters fish with EOD frequencies $f_{1}=631$\,Hz and $f_{2}=797$\,Hz. Both foreign signals have the same strength relative to the own field amplitude (10\,\% contrast). Top row: Sketch of signal processing in the nonlinear system (black box). Second row: Interference of the receiver EOD with the EODs of other fish, bold line highlights the amplitude modulation. Third row: Respective spike trains of the recorded P-unit. Fourth row: Firing rate, estimated by convolution of the spike trains with a Gaussian kernel ($\sigma = 1$\,ms). Bottom row: Power spectrum of the firing rate. \figitem{A} Baseline condition: The cell is driven by the self-generated field alone. The baseline firing rate \fbase{} dominates the power spectrum of the firing rate ($f_{base} = 139$\,Hz). \figitem{B} The receiver's EOD and a foreign fish with an EOD frequency $f_{1}=631$\,Hz are present. EOD interference induces an amplitude modulation, referred to as beat. \figitem{C} The receiver and a fish with an EOD frequency $f_{2}=797$\,Hz are present. The resulting beat is faster as the difference between the individual frequencies is larger. \figitem{D} All three fish with the EOD frequencies \feod{}, $f_{1}$ and $f_{2}$ are present, a second-order amplitude modulation occurs, commonly referred to as envelope. Nonlinear peaks occur at the sum and difference of the two beat frequencies in the power spectrum of the firing rate.
   }
 \end{figure*}
@@ -580,7 +582,7 @@ The population of ampullary cells is generally more homogeneous, with lower base
 
 
 \begin{figure*}[t]
-  %\includegraphics[width=\columnwidth]{model_full.pdf}
+  \includegraphics[width=\columnwidth]{model_full.pdf}
   \caption{\label{fig:model_full} Using second-order susceptibility to predict responses to sine-wave stimuli. \figitem[]{A} Absolute value of the second-order susceptibility,  \eqnref{eq:susceptibility}, for both positive and negative frequencies. \susceptf{} was estimated from $N=10^6$ segments of model simulations in the noise-split condition (cell 2013-01-08-aa, see table~\ref{modelparams} for model parameters). White lines indicate zero frequencies. Nonlinear responses at \fsum{} are quantified in the upper right and lower left quadrants.  Nonlinear responses at \fdiff{} are quantified in the upper left and lower right quadrants. The baseline firing rate of this cell was at $\fbase=120$\,Hz. The position of the orange/red letters corresponds to the beat frequencies used for the stimulation with pure sine waves in the subsequent panels and indicates the sum/difference of those beat frequencies.  \figitem{B--E} Black line -- power spectral density of model simulations in response to stimulation with two pure sine waves, \fone{} and \ftwo, in addition to the receiving fish's own EOD (three-fish scenario). The contrast of beat beats is 0.02.  Colored circles highlight the height of selected peaks in the power spectrum. Grey line -- power spectral density of model in the baseline condition. \figitem{B} The sum of the two beat frequencies match \fbase{}. \figitem{C} The difference of \fone{} and \ftwo{} match \fbase{}. \figitem{D} Only the first beat frequency matches \fbase{}. \figitem{E} None of the two beat frequencies matches \fbase{}.}
 \end{figure*}
 
@@ -760,7 +762,7 @@ Using \eqnref{eq:nli_equation} instead of \eqnref{eq:nli_equation2} for the esti
 
 
 \begin{figure*}[t]
-  %\includegraphics[width=\columnwidth]{flowchart.pdf}
+  \includegraphics[width=\columnwidth]{flowchart.pdf}
   \caption{\label{flowchart}
     Architecture of the P-unit model. Each row illustrates subsequent processing steps for three different stimulation regimes: (i) baseline activity without external stimulus, only the fish's self-generated EOD (the carrier, \eqnref{eq:eod}) is present. (ii) RAM stimulation, \eqnref{eq:ram_equation}. The amplitude of the EOD carrier is modulated with a weak (2\,\% contrast) band-limited white-noise stimulus. (iii) Noise split, \eqnsref{eq:ram_split}--\eqref{eq:Noise_split_intrinsic}, where 90\,\% of the intrinsic noise is replaced by a RAM stimulus, whose amplitude is scaled to maintain the mean firing rate and the CV of the ISIs of the model's baseline activity. As an example, simulations of the model for cell ``2012-07-03-ak'' are shown (table~\ref{modelparams}). \figitem{A} The stimuli are thresholded, \eqnref{eq:threshold2}, by setting all negative values to zero. \figitem{B} Subsequent low-pass filtering, \eqnref{eq:dendrite}, attenuates the carrier and carves out the AM signal. \figitem{C} Intrinsic Gaussian white-noise is added to the signals shown in \panel{B}. Note the reduced internal noise amplitude in the noise split (iii) condition. \figitem{D} Spiking output of the LIF model, \eqnsref{eq:LIF}--\eqref{spikethresh}, in response to the sum of \panel{B} and \panel{C}. \figitem{E} Power spectra of the LIF neuron's spiking activity. Both, baseline activity (\panel[i]{E}) and noise split (\panel[iii]{E}), have the same peaks in the response spectrum at $\fbase$, $f_{EOD} - \fbase$, $f_{EOD}$, and $f_{EOD} + \fbase$. With RAM stimulation (\panel[ii]{E}), the peak at the baseline firing rate, $\fbase$, is washed out.}
 \end{figure*}
@@ -864,26 +866,4 @@ Both, the reduced intrinsic noise and the RAM stimulus, need to replace the orig
 %\bibliography{journalsabbrv,references}
 
 
-\setcounter{figure}{0}
-\renewcommand{\thefigure}{S\arabic{figure}}
-
-\newpage
-
-\section{Supporting information}
-
-\subsection{S1 Second-order susceptibility of high-CV P-unit}
-CVs in P-units can range up to 1.5 \citep{Grewe2017, Hladnik2023}. We show the same analysis as in \figrefb{fig:punit} for an example higher-CV P-unit. Similar to the low-CV cell, high-CV P-units fire at multiples of the EOD period (\subfigrefb{fig:punithighcv}{A}). In contrast to low-CV P-units, however, the higher CV characterizes the noisier baseline firing pattern and the peak at \fbase{} is less pronounced in the power spectrum of the baseline activity (\subfigrefb{fig:punithighcv}{B}). High-CV P-units do not exhibit a clear nonlinear structure related to \fbase{} neither in the second-order susceptibility matrices (\subfigrefb{fig:punithighcv}{E--F}), nor in the projected diagonals (\subfigrefb{fig:punithighcv}{G}). The overall level of nonlinearity, however, shows the same dependence on the stimulus contrast. It is much reduced for high-contrast stimuli that drive the neuron much stronger (\subfigrefb{fig:punithighcv}{F}).
-
-\label{S1:highcvpunit}
-\begin{figure*}[!ht]
-  %\includegraphics[width=\columnwidth]{cells_suscept_high_CV.pdf}  
-  \caption{\label{fig:punithighcv} Response of a noisy P-units (cell ``2018-08-24-af") with a relatively high baseline CV of 0.34 to RAM stimuli with two different contrasts. \figitem{A} Interspike intervals (ISI) distribution during baseline. \figitem{B} Baseline power spectrum. \figitem{C} Top: EOD carrier (gray) with RAM (red). Center: Spike trains in response to the 5\,\% RAM contrast. Bottom: Spike trains in response to the 10\,\% RAM contrast. \figitem{D} First-order susceptibility (\eqnref{linearencoding_methods}). \figitem{E} Absolute value $|\chi_2(f_1, f_2)|$ of the second-order susceptibility, \eqnref{eq:susceptibility}, for the 5\,\% RAM contrast.  Pink lines -- edges of the structure when \fone, \ftwo{} or \fsum{} are equal to \fbase{}. \figitem{F} $|\chi_2(f_1, f_2)|$ for the 10\,\% RAM contrast. \figitem{G} Projected diagonals, calculated as the mean of the anti-diagonals of the matrices in \panel{E--F}. Gray dots: \fbase{}. Dashed lines: Medians of the projected diagonals.}
-\end{figure*}
-
-
-\begin{figure*}[t]
-  \includegraphics[width=\columnwidth]{modelsusceptovern.pdf}
-  \caption{\label{fig:trialnr} Dependence of the estimate of the second-order susceptibility on the number of FFT segments $\n$. \figitem[i]{A} -- \figitem[iv]{A} Second-order susceptibilities estimated from an increasing number of segments $N=10, 100, 1000$ and $10^6$ as indicated. Shown are simulations of the P-unit model cell ``2017-07-18-ai'' (table~\ref{modelparams}) in the noise-split configuration (limit to vanishing external RAM signal, see \figrefb{fig:noisesplit}). \figitem[v]{A} Distribution of second-order susceptibility values as a function of simulated segments. See \panel[iv]{E} for a legend of the percentiles. Here, the estimates converge from about $N=10^4$ segments on. \figitem{B--D} Same for model cells ``2012-12-13-ao'', ``2012-12-20-ac'', and ``2013-01-08-ab'', respectively, same models as in \figrefb{fig:modelsusceptcontrasts}. \figitem{E} Dependence of the estimate of the second-order susceptibility on the number segments for four different contrasts of external RAM stimuli as indicated. Same model cell as in \panel{A}. \figitem[i]{E} Noise-split configuration without additional external RAM stimulus ($c=0$\,\%). Here, the RAM stimulus that was used to estimate the second-order susceptibility replaced $90$\,\% of the intrinsic noise and had a contrast of $10.6$\,\%. Same as \panel[v]{A}, but with a different range of the $y$-axis. \figitem[ii]{E} For a low stimulus contrast of $c=1\,\%$, the estimate does not converge yet even for $N=10^7$ segments. \figitem[iii]{E} For $c=3$\,\% the estimate just starts to converge at $N=10^6$ segments. \figitem[iv]{E} At $c=10$\,\% the situation is similar to the noise split, because similarly strong RAM stimuli were used for estimating the second-order susceptibility.}
-\end{figure*}
-
 \end{document}

plotstyle.py

@@ -151,5 +151,5 @@ def plot_style():
     pt.text_params(font_size=7, font_family='sans-serif', latex=True,
                    preamble=('p:sfmath', 'p:marvosym', 'p:xfrac',
                              'p:SIunits'))
-    pt.ticks_params(xtick_dir='out', xtick_size=3)
+    pt.ticks_params(xtick_dir='out', xtick_size=2.5)
     return ns

punitexamplecell.py

@@ -284,7 +284,7 @@ if __name__ == '__main__':
     axs[0, 0].text(0, 1.6, 'P-units:', transform=axs[0, 0].transAxes,
                    color=s.cell_color1,
             fontsize='large')
-    axs[0, -1].text(1, 0, '5\\,ms', ha='right',
+    axs[0, -1].text(0.97, -0.45, '5\\,ms', ha='right',
                     transform=axs[0, -1].transAxes)
     plot_colorbar(axs[1, -1], pc)
     fig.common_yticks(axs[1, :])

regimes.py

@@ -11,6 +11,21 @@ data_path = Path('data')
 sims_path = data_path / 'simulations'
 
 
+def load_data(file_path):
+    data = np.load(file_path)
+    ratebase = float(data['ratebase'])
+    cvbase = float(data['cvbase'])
+    beatf1 = float(data['beatf1'])
+    beatf2 = float(data['beatf2'])
+    contrasts = data['contrasts']
+    powerf1 = data['powerf1']
+    powerf2 = data['powerf2']
+    powerfsum = data['powerfsum']
+    powerfdiff = data['powerfdiff']
+    return (ratebase, cvbase, beatf1, beatf2,
+            contrasts, powerf1, powerf2, powerfsum, powerfdiff)
+
+
 def load_models(file):
     """ Load model parameter from csv file.
 
@@ -25,7 +40,7 @@ def load_models(file):
         For each cell a dictionary with model parameters.
     """
     parameters = []
-    with file.open('r') as file:
+    with open(file, 'r') as file:
         header_line = file.readline()
         header_parts = header_line.strip().split(",")
         keys = header_parts
@@ -118,7 +133,8 @@ def plot_am(ax, s, alpha, beatf1, beatf2, tmax):
     ax.show_spines('l')
     ax.plot(1000*time, -100*am, **s.lsAM)
     ax.set_xlim(0, 1000*tmax)
-    ax.set_ylim(-50, 50)
+    ax.set_ylim(-13, 13)
+    ax.set_yticks_delta(10)
     #ax.set_xlabel('Time', 'ms')
     ax.set_ylabel('AM', r'\%')
     ax.text(1, 1.2, f'Contrast = {100*alpha:g}\\,\\%',
@@ -134,59 +150,64 @@ def plot_raster(ax, s, spikes, tmax):
     #ax.set_ylabel('Trials')
 
 
-def compute_power(spikes, nfft, dt):
+def compute_power(path, contrast, spikes, nfft, dt):
+    if not path.exists():
+        print(f'    Compute power spectrum for contrast = {100*contrast:4.1f}%')
         psds = []
         time = np.arange(nfft + 1)*dt
         tmax = nfft*dt
-    rates = []
-    cvs = []
         for s in spikes:
-        rates.append(len(s)/tmax)
-        isis = np.diff(s)
-        cvs.append(np.std(isis)/np.mean(isis))
             b, _ = np.histogram(s, time)
+            b = b / dt
             fourier = np.fft.rfft(b - np.mean(b))
             psds.append(np.abs(fourier)**2)
-        #psds.append(fourier)
             freqs = np.fft.rfftfreq(nfft, dt)
-    #print('mean rate', np.mean(rates))
+            prr = np.mean(psds, 0)*dt/nfft
-    #print('CV', np.mean(cvs))
+        np.savez(path, nfft=nfft, deltat=dt, nsegs=len(spikes),
-    return freqs, np.mean(psds, 0)
+                 freqs=freqs, prr=prr)
-    #return freqs, np.abs(np.mean(psds, 0))**2/dt
+    else:
+        print(f'    Load power spectrum for contrast = {100*contrast:4.1f}%')
+        data = np.load(path)
+        freqs = data['freqs']
+        prr = data['prr']
+    return freqs, prr
 
 
 def decibel(x):
     return 10*np.log10(x/1e8)
 
-def plot_psd(ax, s, spikes, nfft, dt, beatf1, beatf2):
+
-    offs = 3
+def plot_psd(ax, s, path, contrast, spikes, nfft, dt, beatf1, beatf2):
-    freqs, psd = compute_power(spikes, nfft, dt)
+    offs = 4
+    freqs, psd = compute_power(path, contrast, spikes, nfft, dt)
     psd /= freqs[1]
     ax.plot(freqs, decibel(psd), **s.lsPower)
     ax.plot(beatf2, decibel(peak_ampl(freqs, psd, beatf2)) + offs,
             label=r'$r$', clip_on=False, **s.psF0)
     ax.plot(beatf1, decibel(peak_ampl(freqs, psd, beatf1)) + offs,
             label=r'$\Delta f_1$', clip_on=False, **s.psF01)
-    ax.plot(beatf2, decibel(peak_ampl(freqs, psd, beatf2)) + offs + 4.5,
+    ax.plot(beatf2, decibel(peak_ampl(freqs, psd, beatf2)) + offs + 5.5,
             label=r'$\Delta f_2$', clip_on=False, **s.psF02)
     ax.plot(beatf2 - beatf1, decibel(peak_ampl(freqs, psd, beatf2 - beatf1)) + offs,
             label=r'$\Delta f_2 - \Delta f_1$', clip_on=False, **s.psF01_2)
     ax.plot(beatf1 + beatf2, decibel(peak_ampl(freqs, psd, beatf1 + beatf2)) + offs,
             label=r'$\Delta f_1 + \Delta f_2$', clip_on=False, **s.psF012)
     ax.set_xlim(0, 300)
-    ax.set_ylim(-40, 0)
+    ax.set_ylim(-60, 0)
     ax.set_xlabel('Frequency', 'Hz')
     ax.set_ylabel('Power [dB]')
 
 
-def plot_example(axs, axr, axp, s, cell, alpha, beatf1, beatf2, nfft, trials):
+def plot_example(axs, axr, axp, s, path, cell, alpha, beatf1, beatf2,
+                 nfft, trials):
+    sim_path = path / f'{cell_name}-contrastspectrum-{1000*alpha:03.0f}.npz'
     dt = 0.0001
     tmax = nfft*dt
     t1 = 0.1
     spikes = punit_spikes(cell, alpha, beatf1, beatf2, tmax, trials)
     plot_am(axs, s, alpha, beatf1, beatf2, t1)
     plot_raster(axr, s, spikes, t1)
-    plot_psd(axp, s, spikes, nfft, dt, beatf1, beatf2)
+    plot_psd(axp, s, sim_path, alpha, spikes, nfft, dt, beatf1, beatf2)
 
 
 def peak_ampl(freqs, psd, f):
@@ -195,39 +216,6 @@ def peak_ampl(freqs, psd, f):
     return np.max(psd_snippet)
 
 
-def compute_peaks(name, cell, alpha_max, beatf1, beatf2, nfft, trials):
-    data_file = sims_path / f'{name}-contrastpeaks.csv'
-    data = TableData(data_file)
-    return data
-    """
-    if data_file.exists():
-        data = TableData(data_file)
-        return data
-    dt = 0.0001
-    tmax = nfft*dt
-    alphas = np.linspace(0, alpha_max, 200)
-    ampl_f1 = np.zeros(len(alphas))
-    ampl_f2 = np.zeros(len(alphas))
-    ampl_sum = np.zeros(len(alphas))
-    ampl_diff = np.zeros(len(alphas))
-    for k, alpha in enumerate(alphas):
-        print(alpha)
-        spikes =  punit_spikes(cell, alpha, beatf1, beatf2, tmax, trials)
-        freqs, psd = compute_power(spikes, nfft, dt)
-        ampl_f1[k] = peak_ampl(freqs, psd, beatf1)
-        ampl_f2[k] = peak_ampl(freqs, psd, beatf2)
-        ampl_sum[k] = peak_ampl(freqs, psd, beatf1 + beatf2)
-        ampl_diff[k] = peak_ampl(freqs, psd, beatf2 - beatf1)
-    data = TableData()
-    data.append('contrast', '%', '%.1f', 100*alphas)
-    data.append('f1', 'Hz', '%g', ampl_f1)
-    data.append('f2', 'Hz', '%g', ampl_f2)
-    data.append('f1+f2', 'Hz', '%g', ampl_sum)
-    data.append('f2-f1', 'Hz', '%g', ampl_diff)
-    data.write(data_file)
-    return data
-    """
-
 def amplitude(power):
     power -= power[0]
     power[power<0] = 0
@@ -254,73 +242,83 @@ def amplitude_squarefit(contrast, power, max_contrast):
     return (r.intercept + r.slope*contrast)**2
 
 
-def plot_peaks(ax, s, data, alphas):
+def plot_peaks(ax, s, alphas, contrasts, powerf1, powerf2, powerfsum,
-    contrast = data[:, 'contrast']
+               powerfdiff):
-    ax.plot(contrast, amplitude_linearfit(contrast, data[:, 'f1'], 4), **s.lsF01m)
+    cmax = 10
-    ax.plot(contrast, amplitude_linearfit(contrast, data[:, 'f2'], 2), **s.lsF02m)
+    contrasts *= 100
-    ax.plot(contrast, amplitude_squarefit(contrast, data[:, 'f1+f2'], 4), **s.lsF012m)
+    ax.plot(contrasts, amplitude_linearfit(contrasts, powerf1, 4),
-    ax.plot(contrast, amplitude_squarefit(contrast, data[:, 'f2-f1'], 4), **s.lsF01_2m)
+            **s.lsF01m)
-    ax.plot(contrast, amplitude(data[:, 'f1']), **s.lsF01)
+    ax.plot(contrasts, amplitude_linearfit(contrasts, powerf2, 2),
-    ax.plot(contrast, amplitude(data[:, 'f2']), **s.lsF02)
+            **s.lsF02m)
-    ax.plot(contrast, amplitude(data[:, 'f1+f2']), **s.lsF012)
+    ax.plot(contrasts, amplitude_squarefit(contrasts, powerfsum, 4),
-    ax.plot(contrast, amplitude(data[:, 'f2-f1']), **s.lsF01_2)
+            **s.lsF012m)
+    ax.plot(contrasts, amplitude_squarefit(contrasts, powerfdiff, 4),
+            **s.lsF01_2m)
+    ax.plot(contrasts, amplitude(powerf1), **s.lsF01)
+    ax.plot(contrasts, amplitude(powerf2), **s.lsF02)
+    mask = contrasts < cmax
+    ax.plot(contrasts[mask], amplitude(powerfsum)[mask],
+            clip_on=False, **s.lsF012)
+    ax.plot(contrasts[mask], amplitude(powerfdiff)[mask],
+            clip_on=False, **s.lsF01_2)
+    ymax = 60
     for alpha, tag in zip(alphas, ['A', 'B', 'C', 'D']):
-        contrast = 100*alpha
+        ax.plot(100*alpha, ymax*0.95, 'vk', ms=4, clip_on=False)
-        ax.plot(contrast, 630, 'vk', ms=4, clip_on=False)
+        ax.text(100*alpha, ymax, tag, ha='center')
-        ax.text(contrast, 660, tag, ha='center')
         #ax.axvline(contrast, **s.lsGrid)
         #ax.text(contrast, 630, tag, ha='center')
-    ax.axvline(1.5, **s.lsLine)
+    ax.axvline(1.2, **s.lsLine)
-    ax.axvline(4, **s.lsLine)
+    ax.axvline(3.5, **s.lsLine)
-    yoffs = 340
+    yoffs = 35
-    ax.text(1.5/2, yoffs, 'linear\nregime',
+    ax.text(1.2/2, yoffs, 'linear\nregime',
             ha='center', va='center')
-    ax.text((1.5 + 4)/2, yoffs, 'weakly\nnonlinear\nregime',
+    ax.text((1.2 + 3.5)/2, yoffs, 'weakly\nnonlinear\nregime',
             ha='center', va='center')
-    ax.text(10, yoffs, 'strongly\nnonlinear\nregime',
+    ax.text(5.5, yoffs, 'strongly\nnonlinear\nregime',
             ha='center', va='center')
-    ax.set_xlim(0, 16.5)
+    ax.set_xlim(0, cmax)
-    ax.set_ylim(0, 600)
+    ax.set_ylim(0, ymax)
-    ax.set_xticks_delta(5)
+    ax.set_xticks_delta(2)
-    ax.set_yticks_delta(300)
+    ax.set_yticks_delta(20)
     ax.set_xlabel('Contrast', r'\%')
     ax.set_ylabel('Amplitude', 'Hz')
 
     
 if __name__ == '__main__':
-    parameters = load_models(data_path / 'punitmodels.csv')
+    cell_name = '2018-05-08-ad-invivo-1'      # 228Hz, CV=0.67
-    cell_name = '2013-01-08-aa-invivo-1'     # 132Hz, CV=0.16: perfect!
+    ratebase, cvbase, beatf1, beatf2, \
-    beatf1 = 40
+        contrasts, powerf1, powerf2, powerfsum, powerfdiff = \
-    beatf2 = 132
+            load_data(sims_path / f'{cell_name}-contrastpeaks.npz')
-    # cell_name = '2012-07-03-ak-invivo-1'    # 128Hz, CV=0.24
+    alphas = [0.002, 0.01, 0.03, 0.06]
-    # cell_name = '2018-05-08-ae-invivo-1'    # 142Hz, CV=0.48
 
+    parameters = load_models(data_path / 'punitmodels.csv')
     cell = cell_parameters(parameters, cell_name)
+    nfft = 2**18
     
-    s = plot_style()
+    print(f'Loaded data for cell {cell_name}: '
+          f'baseline rate = {ratebase:.0f}Hz, CV = {cvbase:.2f}')
 
-    nfft = 2**18
+    s = plot_style()
-    fig, axs = plt.subplots(5, 4, cmsize=(s.plot_width, 0.8*s.plot_width),
+    fig, (axes, axa) = plt.subplots(2, 1, height_ratios=[4, 3],
-                            height_ratios=[1, 1.5, 2, 1.5, 4])
+                                    cmsize=(s.plot_width, 0.6*s.plot_width))
-    fig.subplots_adjust(leftm=8, rightm=2, topm=2, bottomm=3.5,
+    fig.subplots_adjust(leftm=8, rightm=2, topm=2, bottomm=3.5, hspace=0.6)
-                        wspace=0.3, hspace=0.3)
+    axe = axes.subplots(3, 4, wspace=0.4, hspace=0.2,
-    ax0 = fig.merge(axs[3, :])
+                        height_ratios=[1, 2, 3])
-    ax0.set_visible(False)
-    axa = fig.merge(axs[4, :])
     fig.show_spines('lb')
-    alphas = [0.01, 0.03, 0.05, 0.16]
+    
-    #alphas = [0.002, 0.01, 0.05, 0.1]
+    # example power spectra:
     for c, alpha in enumerate(alphas):
-        plot_example(axs[0, c], axs[1, c], axs[2, c], s, cell,
+        plot_example(axe[0, c], axe[1, c], axe[2, c], s, sims_path,
-                     alpha, beatf1, beatf2, nfft, 100)
+                     cell, alpha, beatf1, beatf2, nfft, 100)
-    axs[1, 0].xscalebar(1, -0.1, 30, 'ms', ha='right')
+    axe[1, 0].xscalebar(1, -0.1, 20, 'ms', ha='right')
-    axs[2, 0].legend(loc='center left', bbox_to_anchor=(0, -0.8),
+    axe[2, 0].legend(loc='center left', bbox_to_anchor=(0, -0.8),
                      ncol=5, columnspacing=2)
-    data = compute_peaks(cell_name, cell, 0.2, beatf1, beatf2, nfft, 1000)
+    fig.common_yspines(axe[0, :])
-    plot_peaks(axa, s, data, alphas)
+    fig.common_yticks(axe[2, :])
-    fig.common_yspines(axs[0, :])
+    fig.tag(axe[0, :], xoffs=-3, yoffs=1.6)
-    fig.common_yticks(axs[2, :])
+    
-    #fig.common_xlabels(axs[2, :])
+    # contrast dependence:
-    fig.tag(axs[0, :], xoffs=-2, yoffs=1.6)
+    plot_peaks(axa, s, alphas, contrasts, powerf1, powerf2,
-    fig.tag(axa)
+               powerfsum, powerfdiff)
+    fig.tag(axa, yoffs=2)
     fig.savefig()