Note. The following notebook contains code in addition to text and figures. By default, the code has been hidden. You can click the icon that looks like an eye in the toolbar above to show the code. To run the code, click the cell menu, then "run all".



In [ ]:

    
# Import packages, set preferences, etc.
%matplotlib inline
from brian2 import *
import ipywidgets as ipw
from numpy.random import poisson
from scipy.integrate import quad
from scipy.special import erf
import warnings
warnings.filterwarnings("ignore")

prefs.codegen.target = 'cython'
defaultclock.dt = 0.05*ms



In [ ]:

    
%%html
<!-- hack to improve styling of ipywidgets sliders -->
<style type="text/css">
.widget-label {
    min-width: 35ex;
    max-width: 35ex;
}
.widget-hslider {
    width: 100%;
}
</style>

This notebook demonstrates the analytical solution to the diffusion approximation equations in the basic model. These equations are from Brunel 2000 "Dynamics of sparsely connected networks of excitatory and inhibitory spiking neurons", appendix A, which cites Tuckwell 1988 "Introduction to Theoretical Neurobiology".

Without refractoriness, the mean interspike interval is

$$m=\tau\sqrt{\pi}\int_{-\mu/\sigma}^{(1-\mu)/\sigma}e^{x^2}(1+\mathrm{erf}(x))\,\mathrm{d}x$$

so the firing rate is $1/m$. The CV is

$$CV^2 = 2\pi\tau^2/m^2\int_{-\mu/\sigma}^{(1-\mu)/\sigma}e^{x^2}\int_{-\infty}^x e^{y^2}(1+\mathrm{erf}(y))^2\,\mathrm{d}y\,\mathrm{d}x$$

With refractoriness, the mean interspike interval is

$$\hat{m} = m+t_\mathrm{ref}$$

and the CV is

$$\hat{CV}=CV\;\hat{m}\,/\,m$$

The accuracy of this analytical formulation is demonstrated in the interactive figure below.



In [ ]:

    
def analytical_fr_cv(mu, sigma, tau, refrac):
    ytheta = (1-mu)/sigma
    yr = -mu/sigma
    r0 = 1/(tau*sqrt(pi)*quad(lambda x: exp(x*x)*(1+erf(x)), yr, ytheta)[0])
    c = quad(lambda x: exp(x*x)*quad(lambda y: exp(y*y)*(1+erf(y))**2, -20, x)[0], yr, ytheta)[0]
    cv2 = 2*pi*tau**2*r0**2*c
    cv = sqrt(cv2)
    rate_ref = 1/(1/r0+refrac)
    cv_ref = cv*rate_ref/r0
    return rate_ref, cv_ref

def reduced_model(mu=1.5, sigma=0.5, tau_ms=10, t_ref_ms=0.1):
    # Set parameters
    repeats = 1000
    duration = 1000*ms
    tau = tau_ms*ms
    t_ref = t_ref_ms*ms
    # Define and run the model
    eqs = '''
    dv/dt = (mu-v)/tau+sigma*xi*tau**-0.5 : 1 (unless refractory)
    '''
    G = NeuronGroup(repeats, eqs, threshold='v>1', reset='v=0',
                    refractory=t_ref, method='euler')
    spikemon = SpikeMonitor(G)
    statemon = StateMonitor(G, 'v', record=[0])
    run(duration)
    # Compute ISI histograms
    isi = []
    for train in spikemon.spike_trains().values():
        train.sort()
        isi.append(diff(train))
    isi = hstack(isi)
    cv = std(isi)/mean(isi)
    # Plot results
    figure(figsize=(10, 2.5))
    subplot(131)
    plot(spikemon.t/ms, spikemon.i, ',k')
    xlabel('Time (ms)')
    ylabel('Repeat number')
    title('Spike raster plot')
    xlim(0, duration/ms)
    ylim(0, repeats)
    
    subplot(132)
    plot(statemon.t[:1000]/ms, statemon.v.T[:1000], '-k')
    xlabel('Time (ms)')
    ylabel('v')
    title('Membrane potential trace')
    #xlim(0, duration/ms)
    ylim(-0.2, 1.2)
    axhline(0, ls=':', c='r')
    axhline(1, ls=':', c='g')
    
    subplot(133)
    hist(isi/ms, fc='k', bins=arange(60)*0.5)
    yticks([])
    ylabel('Frequency')
    xlabel('ISI (ms)')
    title('Interspike interval histogram')
    #title('CV = %.2f' % cv)
    text(0.95, 0.9, 'CV = %.2f' % cv, ha='right', va='top',
         bbox=dict(facecolor='white'),
         transform=gca().transAxes)
    tight_layout()
    
    sim_fr = spikemon.num_spikes/(duration*repeats)
    sim_cv = cv
    an_fr, an_cv = analytical_fr_cv(mu, sigma, tau, t_ref)
    print 'Firing rate: simulated=%d sp/s, analytical=%d sp/s' % (sim_fr, an_fr)
    print 'CV: simulated=%.2f, analytical=%.2f' % (sim_cv, an_cv)

display(ipw.interact(reduced_model,
             tau_ms=ipw.FloatSlider(
                 min=0.1, max=20.0, step=0.1, value=10.0,
                 continuous_update=False,
                 description=r"Membrane time constant $\tau$ (ms)"),
             t_ref_ms=ipw.FloatSlider(
                 min=0, max=5, step=0.05, value=0.1,
                 continuous_update=False,
                 description=r"Refractory period $t_\mathrm{ref}$ (ms)"),
             mu=ipw.FloatSlider(
                 min=0, max=5, step=0.05, value=1.5,
                 continuous_update=False,
                 description=r"Mean current $\mu$"),
             sigma=ipw.FloatSlider(
                 min=0, max=5, step=0.05, value=0.5,
                 continuous_update=False,
                 description=r"Standard deviation of current $\sigma$"),
             ));