Solving the integral equation ($N\to\infty$)

Microscopic model

We consider here an infinitely large population of unconnected LIF neurons with escape noise. Between spikes the membrane potential of neuron $i$ obeys: $\tau_m\frac{du_i}{dt}=-u_i+I(t).$ Given $u_i(t)$ at time $t$, neuron $i$ fires a spike stochastically with conditional intensity $\rho_i(t)=c\exp\left(\frac{u_i(t)-\vartheta}{\Delta_u}\right).$ Immediately after a spike a neuron is reset to $u_r=0$ and clamped at this potential for an absolute refractory period $t_{ref}$.

%matplotlib inline
from pylab import *

tref=0.001 #in sec
taum=0.02
mu= 20.
c=10.
du=1.
vth=10.

def hazard(u):
    return c*exp((u-vth)/du)

Macroscopic model

$A(t)=\int_{-\infty}^t\rho(t|\hat{t})S(t|\hat{t})A(\hat{t})\,d\hat{t},\qquad \text{$A(t)=\delta(t)$ for $t\le 0$}$ The integral is discretized and solved forward in time with time step $dt$.

def solve(tmax,dt,I):
    Nt=int(tmax/dt+0.5)
    T=5*taum  #max ISI, in sec
    kmax=int(T/dt+0.5)
    tauv=dt*np.arange(kmax+1)
    rho=zeros(kmax+1)
    n=zeros(Nt)   #A(t)* dt

    m=zeros(kmax+1) #refractory states, age vector
    m[1]=1.

    ufree=mu
    mfree=0.
    rhofree=hazard(ufree)
    u=zeros(kmax+1)
    kref=int(tref/dt+0.5)
    for i in range(Nt):
        mfree += m[kmax]     
        for k in range(kmax)[::-1]:
            if tauv[k]>=tref:
                u[k+1]=u[k]+dt*(-u[k]+I[i])/taum
                rho_new=hazard(u[k+1])
                Pfire = 1-exp(-0.5*(rho[k]+rho_new)*dt)
                rho[k+1]=rho_new
                dm = Pfire * m[k] #spike contribution from bin k
                m[k+1] = m[k]-dm
                n[i] += dm
            else:
                m[k+1]=m[k]

        ufree += (-ufree + I[i])*dt/taum
        rhofree_new=hazard(ufree)
        dmfree= (1-exp(-0.5*(rhofree+rhofree_new)*dt)) * mfree
        rhofree = rhofree_new
        mfree -= dmfree
        n[i] += dmfree
        m[0] = n[i]

    return n/dt

We can now solve the model:

dt=0.001
tmax=0.5 #simulation time
Nt=int(tmax/dt+0.5)
tv=dt* np.arange(Nt)
I = mu*np.ones(Nt)+10*(tv>0.3)-20*(tv>0.4)
A = solve(tmax, dt,I)

and plot the population activity

figure(1)
clf()
plot(tv,A,label='population act.')
plot(tv,I,'r',label='Stimulus [mV]')
xlabel('t [s]')
ylabel('A [Hz]')
legend(loc=0)

<matplotlib.legend.Legend at 0x7efc6d987410>

png

Problem: Compute analytical prediction of stationary firing rate (e.g., in order to check the simulation)

def get_singleneuron_rate_theory(mu,V_th,V_reset,t_ref,tau_m, delta_u,rho_0):
    tmax=0.
    dt=0.0001
    S_end=1
    while (S_end>0.001 and tmax<10.):
        tmax=tmax+1  #max ISI, in sec
        t=np.arange(0,tmax,dt)
        K=len(t)
        S=np.ones(K)
        v=mu+(V_reset-mu) * np.exp(-(t-t_ref)/tau_m)
        rho=rho_0*np.exp((v-V_th)/delta_u)*(t>t_ref)
        S=np.exp(-np.cumsum(rho)*dt)
        S_end=S[-1]
    return 1./(np.sum(S)*dt)

print(get_singleneuron_rate_theory(mu,vth,0,tref,taum, du,c))    

44.49762905220395

Problem: How does the time-dependent response to the stimulus change when the noise is increased?