MLI Current parameters tuning

Goal: choose parameters for a random current injected into the MLI every time step such that the model exhibits ISI distribution statistics that match those described in Hausser and Clark (1997) in the presence of GABA blockers.

  • Hausser and Clark (1997) Figure 1F show that under control conditions, there is a distribution of mean firing rates among MLIs
    • We refer to this as a 'prior distribution' and make no attempt to model it. Instead we model the parameters for a single neuron -- the average neuron -- having mean firing rate described below.
  • Hausser and Clark (1997) report mean firing rates of 19.1 Hz (= 12.3 * 1.55) and CV of .24
  • Hausser and Clark (1997) Figure 3 shows an example ISI histogram of an MLI that is symmetric (Gaussian looking) with mean firing rate of 30 Hz, CV of .14.
  • We model random currents injected every time step drawn from a gamma distribution with parameters k,theta chosen to match 30Hz / .14
    • We choose parameters and generate a representative example for display at 30 Hz.
  • We use a grid search over a range of parameters of the probability distribution and measure their effect (mean firing rate and CV) of the neuron's ISI distribution.
    • This grid search is performed in a separate python script to take advantage of parallel processing.
    • The results are loaded in and analyzed using Pandas

In [2]:
%load_ext autoreload
%autoreload 2
import time
import sys
sys.path.append('../../')
from brian import *
from MLI_PKJ_net import *
from statsmodels.tsa.stattools import acf
from pandas import *
set_global_preferences(useweave=True, usenewpropagate=True, usecodegen=True, usecodegenweave=True)
defaultclock.dt = .25*ms


/home/bill/anaconda/lib/python2.7/site-packages/brian/utils/sparse_patch/__init__.py:39: UserWarning: Couldn't find matching sparse matrix patch for scipy version 0.13.3, but in most cases this shouldn't be a problem.
  warnings.warn("Couldn't find matching sparse matrix patch for scipy version %s, but in most cases this shouldn't be a problem." % scipy.__version__)
/home/bill/anaconda/lib/python2.7/site-packages/brian/synapses/spikequeue.py:490: UserWarning: Using C++ SpikeQueue
  warnings.warn('Using C++ SpikeQueue')

Data analysis from current parameter sweep


In [1]:
# This is the output file from the MLI gamma distribution parameter sweep python script
data_dir = '~/data/neuron_models/molecular_layer/mli_gamma_current_param_sweep/2013-10-17T11:10:42.982338/results.txt'

In [2]:
with open(data_dir) as f:
    results = read_csv(f,delimiter='\t')
target = array([30.,.14])
results['error'] = ((results.ix[:,-2:] - target).abs()/target).sum(axis=1)
results['mean'] = results.ix[:,0] + results.ix[:,1]
results.sort(columns='error')[:5]


Out[2]:
k theta mli_mean_firing_rate mli_cv error mean
2991 3.966333 0.006653 29.666667 0.139010 0.018180 3.972985
3099 4.293988 0.006148 29.500000 0.138495 0.027417 4.300136
2863 3.616834 0.007157 29.333333 0.139248 0.027592 3.623990
2707 3.223647 0.008165 29.333333 0.140755 0.027612 3.231812
3085 4.250301 0.006148 29.000000 0.140900 0.039764 4.256449

Run network with parameters


In [2]:
def run_net(T,N_MLI,k,theta):
    reinit()
    reinit_default_clock()
    clear(True)
    
    MLI = MLIGroup(N_MLI)
    
    @network_operation(Clock(dt=defaultclock.dt))
    def random_current():
        MLI.I = gamma(k,theta,size=len(MLI)) * nA
        
    # Monitor
    MS_MLI = SpikeMonitor(MLI)
    MR_MLI = PopulationRateMonitor(MLI,bin=1*ms)
    MISI_MLI = ISIHistogramMonitor(MLI,bins=arange(0,162,2)*ms)
    MV_MLI = StateMonitor(MLI,'V',record=0)
    
    run(T*msecond)
    
    return MS_MLI, MR_MLI, MISI_MLI, MV_MLI

Case #1


In [3]:
T = 300000
N_MLI = 1

# chosen parameters
k, theta = 3.966333,0.006653

# run simulation
MS_MLI, MR_MLI, MISI_MLI, MV_MLI = run_net(T,N_MLI,k,theta)


brian.stateupdater: WARNING  Using codegen CStateUpdater

In [10]:
# save monitors for reproducibility
monitors = {'MS_MLI_iso':MS_MLI, 'MR_MLI_iso':MR_MLI, 'MISI_MLI_iso':MISI_MLI, 'MV_MLI_iso':MV_MLI}
out_dir = '~/'
for name, mon in monitors.iteritems():
    cPickle.dump(mon, open(out_dir+'%s.pkl'%name,'w'))

In [7]:
load_monitors = False
if load_monitors:
    in_dir = '~/'
    MS_MLI = cPickle.load(open(in_dir+'MS_MLI_iso.pkl'))
    MV_MLI = cPickle.load(open(in_dir+'MV_MLI_iso.pkl'))
    MISI_MLI = cPickle.load(open(in_dir+'MISI_MLI_iso.pkl'))
    MR_MLI = cPickle.load(open(in_dir+'MR_MLI_iso.pkl'))

In [9]:
fig = figure(figsize=(25,5), dpi=600)
# Plot trace
ax = fig.add_subplot(131)
MV_MLI.insert_spikes(MS_MLI)
ax.plot(MV_MLI.times[:6000],(MV_MLI.values[0,:6000]), color='#8C2318')
xlim([-.1,1.5])
add_scalebar(ax,labelx='This')
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
ax.spines['left'].set_visible(False)
ax.spines['bottom'].set_visible(False)

# Plot ISI
ax = fig.add_subplot(132)
plot_ISI_histogram(MISI_MLI, MS_MLI, MR_MLI, (20,0),(95,MISI_MLI.count.max()*.7), color='#8C2318', edgecolor='w')
simpleaxis(ax)
xlim([0,200])
tick_params(labelsize=18)

# Plot spike autocorrelation
ax = fig.add_subplot(133)
plot_spike_correlogram(MS_MLI.spiketimes[0],MS_MLI.spiketimes[0], bin=1*ms, width=200*ms, ax=ax, color='#8C2318')
simpleaxis(ax)
tick_params(labelsize=18)

tight_layout()
fig.subplots_adjust(wspace=.3)
fig.savefig('MLI_iso_color.tiff',dpi=600*6.93/25)
fig.savefig('MLI_iso_color.png',dpi=600*6.93/25)



In [17]:
# statistical tests
from scipy.stats import normaltest, shapiro
isis = diff(MS_MLI.spiketimes[0])
shapiro(isis[:5000])


Out[17]:
(0.9462147355079651, 2.3953936079014855e-39)