145 lines
5.1 KiB
Python
145 lines
5.1 KiB
Python
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()
|