started noisesplit figure

ampullaryexamplecell.py

@@ -94,7 +94,7 @@ def plot_chi2(ax, s, contrast, freqs, chi2, fcutoff, vmax):
     if vmax is None:
         vmax = np.quantile(1e-3*chi2, 0.99)
     pc = ax.pcolormesh(freqs, freqs, 1e-3*chi2, vmin=0, vmax=vmax,
-                       cmap='viridis', rasterized=True, zorder=10)
+                       rasterized=True, zorder=10)
     ax.set_xlim(0, fcutoff)
     ax.set_ylim(0, fcutoff)
     df = 100 if fcutoff == 300 else 50
@@ -220,3 +220,4 @@ if __name__ == '__main__':
     fig.tag(axs, xoffs=-3, yoffs=2)
     
     fig.savefig()
+    print()

dataoverview.py

@@ -274,3 +274,4 @@ if __name__ == '__main__':
     fig.common_yticks(axs[2, :])
     fig.tag(xoffs=-3.5, yoffs=2)
     fig.savefig()
+    print()

lifsuscept.py

@@ -0,0 +1,144 @@
+import numpy as np
+import mpmath as mp
+import matplotlib.pyplot as plt
+from pathlib import Path
+from plotstyle import plot_style
+
+
+sims_path = Path('data') / 'simulations'
+
+
+"""
+LIF code from Maria Schlungbaum, Lidner lab, 2024
+
+LIF model in dimensionless units: dv/dt = -v + mu + sqrt(2D)*xi
+v: membrane voltage
+mu: mean input voltage
+D: noise intensity
+xi: white Gaussian noise
+tau_mem = 1 (membrane time constant, skipped here)
+tau_ref: refractory period
+vT: threshold voltage
+vR: reset voltage
+"""
+
+
+def firingrate(mu, D, tau_ref, vR, vT):
+    x_start = (mu - vT)/mp.sqrt(2.0*D)
+    x_end = (mu - vR)/mp.sqrt(2.0*D)
+    dx = 0.0001
+    r = 0.0
+    for i in np.arange(x_start, x_end, dx):
+        integrand = mp.exp(i**2) * mp.erfc(i)
+        r += integrand*dx
+    r0 = 1.0/(tau_ref + mp.sqrt(mp.pi)*r)
+    return float(r0)
+
+
+def susceptibility1(omega, r0, mu, D, tau_ref, vR, vT):
+    delta = (vR**2 - vT**2 + 2.0*mu*(vT - vR))/(4.0*D)
+    a = (r0 * omega*1.0j)/(mp.sqrt(D) * (omega*1.0j - 1.0))
+    b = mp.pcfd(omega*1.0j - 1.0, (mu - vT)/mp.sqrt(D)) - mp.exp(delta) * mp.pcfd(omega*1.0j - 1.0, (mu - vR)/mp.sqrt(D))
+    c = mp.pcfd(omega*1.0j, (mu - vT)/mp.sqrt(D)) - mp.exp(delta) * mp.exp(omega*1.0j*tau_ref) * mp.pcfd(omega*1.0j, (mu - vR)/mp.sqrt(D))
+    return a * b/c
+
+
+def susceptibility2(omega1, omega2, chi1_1, chi1_2, r0, mu, D, tau_ref, vR, vT):
+    delta = (vR**2 - vT**2 + 2.0*mu*(vT - vR))/(4.0*D)
+    a1 = r0*(1.0 - omega1*1.0j - omega2*1.0j)*(omega1*1.0j + omega2*1.0j)/(2.0*D*(omega1*1.0j - 1.0)*(omega2*1.0j - 1.0))
+    a2 = (omega1*1.0j + omega2*1.0j)/(2.0*mp.sqrt(D))
+    a3 = chi1_1/(omega2*1.0j - 1.0) + chi1_2/(omega1*1.0j - 1.0)
+    b1 = mp.pcfd(omega1*1.0j + omega2*1.0j - 2.0, (mu - vT)/mp.sqrt(D)) - mp.exp(delta) * mp.pcfd(omega1*1.0j + omega2*1.0j - 2.0, (mu - vR)/mp.sqrt(D))
+    b2 = mp.pcfd(omega1*1.0j + omega2*1.0j - 1.0, (mu - vT)/mp.sqrt(D))
+    b3 = mp.exp(delta) * mp.pcfd(omega1*1.0j + omega2*1.0j - 1.0, (mu - vR)/mp.sqrt(D))
+    c = mp.pcfd(omega1*1.0j + omega2*1.0j, (mu - vT)/mp.sqrt(D)) - mp.exp(delta) * mp.exp(1.0j*(omega1 + omega2)*tau_ref) * mp.pcfd(omega1*1.0j + omega2*1.0j, (mu - vR)/mp.sqrt(D))
+    return a1 * b1/c + a2*a3*b2/c - a2*a3*b3/c
+
+
+def susceptibilities(frange, mu, D, tau_ref, vR, vT):
+    print('compute LIF susceptibilites:')
+    r0 = firingrate(mu, D, tau_ref, vR, vT)
+    chi1_data = np.zeros(len(frange), dtype=complex)
+    chi2_data = np.zeros((len(frange), len(frange)), dtype=complex)
+    for f2 in range(len(frange)):
+        print(f'    step {f2 + 1:4d} of {len(frange):4d}')
+        omega2 = 2.0*np.pi*frange[f2]
+        chi1_2 = susceptibility1(omega2, r0, mu, D, tau_ref, vR, vT)
+        chi1_data[f2] = chi1_2
+        for f1 in range(len(frange)):
+            omega1 = 2.0*np.pi*frange[f1]
+            chi1_1 = susceptibility1(omega1, r0, mu, D, tau_ref, vR, vT)
+            chi2 = susceptibility2(omega1, omega2, chi1_1, chi1_2, r0, mu, D, tau_ref, vR, vT)
+            chi2_data[f2, f1] = chi2
+    return r0, chi1_data, chi2_data
+
+        
+def load_lifdata(mu, D, vT=1, vR=0, tau_ref=0):
+    file_path = sims_path / f'lif-mu{10*mu:03.0f}-D{10000*D:04.0f}-chi2.npz'
+    if not file_path.exists():
+        freqs = np.linspace(0.01, 1.0, 200)
+        r0, chi1, chi2 = susceptibilities(freqs, mu, D, tau_ref, vR, vT)
+        np.savez(file_path, mu=mu, D=D, vT=vT, vR=vR, tau_mem=1, tau_ref=tau_ref,
+                 r0=r0, freqs=freqs, chi1=chi1, chi2=chi2)
+    data = np.load(file_path)
+    r0 = float(data['r0'])
+    freqs = data['freqs']
+    chi1 = data['chi1']
+    chi2 = data['chi2']
+    return r0, freqs, chi1, chi2
+
+
+def plot_gain(ax, s, r0, freqs, chi1):
+    ax.plot(freqs, np.abs(chi1), **s.lsM1)
+    ax.set_xlabel('$f$')
+    ax.set_ylabel('$|\\chi_1(f)|$', labelpad=6)
+    ax.set_xlim(0, 1)
+    ax.set_ylim(0, 14)
+    ax.set_xticks_delta(0.2)
+    ax.set_yticks_delta(3)
+
+
+def plot_chi2(ax, s, r0, freqs, chi2):
+    chi2 = np.abs(chi2)
+    vmax = np.quantile(chi2, 0.996)
+    ten = 10**np.floor(np.log10(vmax))
+    for fac, delta in zip([1, 2, 3, 4, 6, 8, 10],
+                          [0.5, 1, 1, 2, 3, 4, 5]):
+        if fac*ten >= vmax:
+            vmax = fac*ten
+            ten *= delta
+            break
+    pc = ax.pcolormesh(freqs, freqs, chi2, vmin=0, vmax=vmax,
+                       rasterized=True)
+    ax.set_aspect('equal')
+    ax.set_xlabel('$f_1$')
+    ax.set_ylabel('$f_2$', labelpad=6)
+    ax.set_xlim(0, 1)
+    ax.set_ylim(0, 1)
+    ax.set_xticks_delta(0.2)
+    ax.set_yticks_delta(0.2)
+    cax = ax.inset_axes([1.04, 0, 0.05, 1])
+    cax.set_spines_outward('lrbt', 0)
+    cb = fig.colorbar(pc, cax=cax)
+    cb.outline.set_color('none')
+    cb.outline.set_linewidth(0)
+    cax.set_ylabel('$|\\chi_2(f_1, f_2)|$')
+    cax.set_yticks_delta(ten)
+
+    
+if __name__ == '__main__':
+    mu = 1.1
+    D = 0.001
+
+    r0, freqs, chi1, chi2 = load_lifdata(mu, D)
+
+    s = plot_style()
+    plt.rcParams['axes.labelpad'] = 2
+    fig, (axg, axc) = plt.subplots(1, 2, cmsize=(s.plot_width, 0.38*s.plot_width))
+    fig.subplots_adjust(leftm=8, rightm=8.5, topm=1, bottomm=3.5, wspace=0.4)
+    fig.set_align(autox=False)
+    plot_gain(axg, s, r0, freqs, chi1)
+    plot_chi2(axc, s, r0, freqs, chi2)
+    fig.tag()
+    fig.savefig()
+    print()

modelsusceptcontrasts.py

@@ -47,7 +47,7 @@ def plot_chi2(ax, s, data_file):
             ten *= delta
             break
     pc = ax.pcolormesh(freqs, freqs, prss, vmin=0, vmax=vmax,
-                       cmap='viridis', rasterized=True)
+                       rasterized=True)
     if 'noise_frac' in data:
         ax.set_title('$c$=0\\,\\%', fontsize='medium')
     else:
@@ -182,3 +182,4 @@ if __name__ == '__main__':
     fig.common_yticks(axs[5, :])
     fig.tag(axs, xoffs=-4.5, yoffs=1.8)
     fig.savefig()
+    print()

modelsusceptlown.py

@@ -47,7 +47,7 @@ def plot_chi2(ax, s, data_file):
             ten *= delta
             break
     pc = ax.pcolormesh(freqs, freqs, prss, vmin=0, vmax=vmax,
-                       cmap='viridis', rasterized=True)
+                       rasterized=True)
     ns = f'$N={n}$' if n <= 100 else f'$N=10^{np.log10(n):.0f}$'
     if 'noise_frac' in data:
         ax.set_title(f'$c$=0\\,\\%, {ns}', fontsize='medium')
@@ -197,3 +197,4 @@ if __name__ == '__main__':
     fig.tag(axs, xoffs=-4.5, yoffs=1.8)
     axs[1, 0].set_visible(False)
     fig.savefig()
+    print()

modelsusceptovern.py

@@ -46,7 +46,7 @@ def plot_chi2(ax, s, data_file):
             ten *= delta
             break
     pc = ax.pcolormesh(freqs, freqs, prss, vmin=0, vmax=vmax,
-                       cmap='viridis', rasterized=True)
+                       rasterized=True)
     ax.set_title(f'$N=10^{np.log10(n):.0f}$', fontsize='medium')
     ax.set_xlim(0, 300)
     ax.set_ylim(0, 300)
@@ -164,3 +164,5 @@ if __name__ == '__main__':
     fig.common_yticks(axs[-1, :4])
     fig.tag(axs, xoffs=-2.5, yoffs=1.8)
     fig.savefig()
+    print()
+    

noisesplit.py

@@ -0,0 +1,150 @@
+import numpy as np
+import matplotlib.pyplot as plt
+from pathlib import Path
+from plotstyle import plot_style
+
+
+base_path = Path('data')
+data_path = base_path / 'cells'
+sims_path = base_path / 'simulations'
+
+
+def sort_files(cell_name, all_files, n):
+    files = [fn for fn in all_files if '-'.join(fn.stem.split('-')[2:-n]) == cell_name]
+    if len(files) == 0:
+        return None, 0
+    nums = [int(fn.stem.split('-')[-1]) for fn in files]
+    idxs = np.argsort(nums)
+    files = [files[i] for i in idxs]
+    nums = [nums[i] for i in idxs]
+    return files, nums
+
+
+def plot_chi2(ax, s, freqs, chi2, nsegs):
+    ax.set_aspect('equal')
+    i0 = np.argmin(freqs < 0)
+    i1 = np.argmax(freqs > 300)
+    if i1 == 0:
+        i1 = len(freqs)
+    freqs = freqs[i0:i1]
+    chi2 = chi2[i0:i1, i0:i1]
+    vmax = np.quantile(chi2, 0.996)
+    ten = 10**np.floor(np.log10(vmax))
+    for fac, delta in zip([1, 2, 3, 4, 6, 8, 10],
+                          [0.5, 1, 1, 2, 3, 4, 5]):
+        if fac*ten >= vmax:
+            vmax = fac*ten
+            ten *= delta
+            break
+    pc = ax.pcolormesh(freqs, freqs, chi2, vmin=0, vmax=vmax,
+                       rasterized=True)
+    ax.set_xlim(0, 300)
+    ax.set_ylim(0, 300)
+    ax.set_xticks_delta(100)
+    ax.set_yticks_delta(100)
+    ax.set_xlabel('$f_1$', 'Hz')
+    ax.set_ylabel('$f_2$', 'Hz')
+    ax.text(1, 1.1, f'$N=10^{{{np.log10(nsegs):.0f}}}$',
+            ha='right', transform=ax.transAxes)
+    
+    cax = ax.inset_axes([1.04, 0, 0.05, 1])
+    cax.set_spines_outward('lrbt', 0)
+    cb = fig.colorbar(pc, cax=cax)
+    cb.outline.set_color('none')
+    cb.outline.set_linewidth(0)
+    cax.set_ylabel(r'$|\chi_2|$ [Hz]')
+    cax.set_yticks_delta(ten)
+
+
+def plot_chi2_contrast(ax1, ax2, s, cell_name, contrast, nsmall, nlarge):
+    data_files = sims_path.glob(f'chi2-noisen-{cell_name}-{1000*contrast:03.0f}-*.npz')
+    files, nums = sort_files(cell_name, data_files, 2)
+    for ax, n in zip([ax1, ax2], [nsmall, nlarge]):
+        i = nums.index(n)
+        data = np.load(files[i])
+        n = data['n']
+        alpha = data['alpha']
+        freqs = data['freqs']
+        pss = data['pss']
+        chi2 = np.abs(data['prss'])*0.5/np.sqrt(pss.reshape(1, -1)*pss.reshape(-1, 1))
+        plot_chi2(ax, s, freqs, chi2, n)
+
+
+def plot_chi2_split(ax1, ax2, s, cell_name, nsmall, nlarge):
+    data_files = sims_path.glob(f'chi2-split-{cell_name}-*.npz')
+    files, nums = sort_files(cell_name, data_files, 1)
+    for ax, n in zip([ax1, ax2], [nsmall, nlarge]):
+        i = nums.index(n)
+        data = np.load(files[i])
+        n = data['n']
+        alpha = data['alpha']
+        noise_frac = data['noise_frac']
+        freqs = data['freqs']
+        pss = data['pss']
+        chi2 = np.abs(data['prss'])*0.5/np.sqrt(pss.reshape(1, -1)*pss.reshape(-1, 1))
+        plot_chi2(ax, s, freqs, chi2, n)
+    return alpha, noise_frac
+
+        
+def plot_chi2_data(ax, s, cell_name, run):
+    data_file = data_path / f'{cell_name}-spectral-s{run:02d}.npz'
+    data = np.load(data_file)
+    n = data['n']
+    alpha = data['alpha']
+    freqs = data['freqs']
+    pss = data['pss']
+    chi2 = np.abs(data['prss'])*0.5/np.sqrt(pss.reshape(1, -1)*pss.reshape(-1, 1))
+    print(f'Measured cell {data_file.name} at {100*alpha:.1f}% contrast')
+    plot_chi2(ax, s, freqs, chi2, n)
+    return alpha
+
+    
+def plot_noise_split(ax, contrast, noise_contrast, noise_frac):
+    axr, axs, axn = ax.subplots(3, 1, hspace=0.1)
+    tmax = 50
+    
+    axr.show_spines('l')
+    axr.set_xlim(0, tmax)
+    axr.set_ylim(-8, 8)
+    axr.set_yticks_delta(6)
+    axr.set_ylabel('\\%')
+    
+    axs.show_spines('l')
+    axs.set_xlim(0, tmax)
+    axs.set_ylim(-8, 8)
+    axs.set_yticks_delta(6)
+    axs.set_ylabel('\\%')
+    
+    axn.set_ylim(-6, 6)
+    axn.set_xlim(0, tmax)
+    axn.set_yticks_delta(6)
+    axn.set_yticks_blank()
+    axn.set_xticks_delta(25)
+    axn.set_xlabel('Time', 'ms')
+
+    
+if __name__ == '__main__':
+    cell_name = '2012-07-03-ak-invivo-1'
+    nsmall = 100
+    nlarge = 1000000
+    contrast = 0.03
+                 
+    s = plot_style()
+    fig, axs = plt.subplots(2, 4, cmsize=(s.plot_width, 0.4*s.plot_width),
+                            width_ratios=[1, 0, 1, 1, 1])
+    fig.subplots_adjust(leftm=7, rightm=8, topm=2, bottomm=3.5,
+                        wspace=0.4, hspace=0.6)
+    axs[1, 0].set_visible(False)
+    data_contrast = plot_chi2_data(axs[0, 0], s, cell_name[:13], 0)
+    plot_noise_split(axs[0, 1], data_contrast, 0, 1)
+    plot_chi2_contrast(axs[0, 2], axs[0, 3], s, cell_name, contrast, nsmall, nlarge)
+    noise_contrast, noise_frac = plot_chi2_split(axs[1, 2], axs[1, 3], s,
+                                                 cell_name, nsmall, nlarge)
+    plot_noise_split(axs[1, 1], contrast, noise_contrast, noise_frac)
+    fig.common_xticks(axs[:, 2])
+    fig.common_xticks(axs[:, 3])
+    fig.common_yticks(axs[0, 2:])
+    fig.common_yticks(axs[1, 2:])
+    #fig.tag(axs, xoffs=-4.5, yoffs=1.8)
+    fig.savefig()
+    print()

nonlinearbaseline.tex

@@ -417,7 +417,7 @@ We here analyze nonlinear responses in two types of primary electroreceptor affe
 \section{Introduction}
 
 \begin{figure*}[t]
-  %\includegraphics[width=\columnwidth]{plot_chi2.pdf}
+  \includegraphics[width=\columnwidth]{lifsuscept}
   \caption{\label{fig:lifresponse} First- (linear) and second-order response functions of the leaky integrate-and-fire model. \figitem{A} Magnitude of the first-order response function $|\chi_1(f)|$, also known as the ``gain'' function, quantifies the response amplitude relative to the stimulus amplitude, both measured at the stimulus frequency. \figitem{B} Magnitude of the second-order response function $|\chi_2(f_1, f_2)|$ quantifies the response at the sum of two stimulus frequencies. For linear systems, the second-order response function is zero, because linear systems do not create new frequencies and thus there is no response at the sum of the two frequencies. The plots show the analytical solutions from \citep{Lindner2001} and \citep{Voronenko2017} with $\mu = 1.1$ and $D = 0.001$. Note that the leaky integrate-and-fire model is formulated without dimensions, frequencies are given in multiples of the inverse membrane time constant.}
 \end{figure*}
 We like to think about signal encoding in terms of linear relations with unique mapping of a given input value to a certain output of the system under consideration. Indeed, such linear methods, for example the transfer function or first-oder susceptibility shown in figure~\ref{fig:lifresponse}, have been widely and successfully applied to describe and predict neuronal responses and are an invaluable tool to characterize neural systems \citep{Borst1999}. Nonlinear mechanisms, on the other hand, are key on different levels of neural processing. Deciding for one action over another is a nonlinear process on the systemic level. On the cellular level, spiking neurons are inherently nonlinear. Whether an action potential is elicited depends on the membrane potential to exceed a threshold \citep{Hodgkin1952, Koch1995}. Because of such nonlinearities, understanding and predicting neuronal responses to sensory stimuli is in general a difficult task.
@@ -495,7 +495,7 @@ Electric fish possess an additional electrosensory system, the passive or ampull
 In the example recordings shown above (\figsrefb{fig:punit} and \fref{fig:ampullary}), we only observe nonlinear responses on the anti-diagonal of the second-order susceptibility, where the sum of the two stimulus frequencies matches the neuron's baseline firing rate, which is in line with theoretical expectations \citep{Voronenko2017,Franzen2023}. However, a pronounced nonlinear response at frequencies \foneb{} or \ftwob{}, although predicted by theory, cannot be observed. In the following, we investigate how these discrepancies can be understood.
 
 \begin{figure*}[t]
-  %\includegraphics[width=\columnwidth]{model_and_data.pdf}
+  \includegraphics[width=\columnwidth]{noisesplit}
   \notejb{This model in the next figure shows a triangle for 3\% contrast ...}
   \notejb{We cannot really make up this twist with the 3\% contrast not converging into a triangle.}
   \caption{\label{fig:noisesplit} Estimation of second-order susceptibilities in the limit of weak stimuli. \figitem{A} \suscept{} estimated from $N=11$ 0.5\,s long segments of an electrophysiological recording of another low-CV P-unit (cell 2012-07-03-ak, $\fbase=120$\,Hz, CV=0.20) driven with a weak RAM stimulus with contrast 2.5\,\%.  Pink edges mark the baseline firing rate where enhanced nonlinear responses are expected. \figitem[i]{B} \textit{Standard condition} of model simulations with intrinsic noise (bottom) and a RAM stimulus (top). \figitem[ii]{B} \suscept{} estimated from simulations of the cell's LIF model counterpart (cell 2012-07-03-ak, table~\ref{modelparams}) based on the same number of FFT segments $N=11$ as in the electrophysiological recording. \figitem[iii]{B} Same as \panel[ii]{B} but using $10^6$ segments. \figitem[i-iii]{C} Same as in \panel[i-iii]{B} but in the \textit{noise split} condition: there is no external RAM signal driving the model. Instead, a large part (90\,\%) of the total intrinsic noise is treated as a signal and is presented as an equivalent amplitude modulation (\signalnoise, center), while the intrinsic noise is reduced to 10\,\% of its original strength (see methods for details). Simulating one million segments, this reveals the full expected trangular structure of the second-order susceptibility.}

plotstyle.py

@@ -126,7 +126,7 @@ def plot_style():
     pt.colormap('YR', cmcolors, cmvalues)
     cycle_colors = ['blue', 'red', 'orange', 'lightgreen', 'magenta',
                     'yellow', 'cyan', 'pink']
-    pt.colors_params(palette, cycle_colors, 'RdYlBu')
+    pt.colors_params(palette, cycle_colors, 'viridis')
     pt.figure_params(palette['white'], format='pdf', compression=6,
                      fonttype=3, stripfonts=True)
     pt.labels_params('{label} [{unit}]', labelsize='medium', labelpad=6)

punitexamplecell.py

@@ -160,7 +160,7 @@ def plot_chi2(ax, s, contrast, freqs, chi2, fcutoff, vmax):
     if vmax is None:
         vmax = np.quantile(1e-3*chi2, 0.99)
     pc = ax.pcolormesh(freqs, freqs, 1e-3*chi2, vmin=0, vmax=vmax,
-                       cmap='viridis', rasterized=True, zorder=10)
+                       rasterized=True, zorder=10)
     ax.set_xlim(0, fcutoff)
     ax.set_ylim(0, fcutoff)
     df = 100 if fcutoff == 300 else 50
@@ -272,3 +272,4 @@ if __name__ == '__main__':
     fig.tag(axs, xoffs=-3, yoffs=2)
     
     fig.savefig()
+    print()