diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..a3066f8 --- /dev/null +++ b/.gitignore @@ -0,0 +1,10 @@ +*.dll +mods/*.c +mods/*.o + +x86_64/* +__pycache__/* + +data/* +figures/* + diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md new file mode 100644 index 0000000..9ebb0ff --- /dev/null +++ b/CONTRIBUTING.md @@ -0,0 +1,6 @@ +# Contributing + +1. **NEVER PUSH TO MASTER:** For each new feature you want to add, create a new branch. +2. When you want to merge a feature into the main branch, submit a pull request. +3. Ask for code review by Abe when you want to merge into the main branch. +4. Add an _Issue_ if you encounter problems or if you want to add new features, todos, etc. \ No newline at end of file diff --git a/FitSimScore_ForallFigs.m b/FitSimScore_ForallFigs.m deleted file mode 100644 index 4b72264..0000000 --- a/FitSimScore_ForallFigs.m +++ /dev/null @@ -1,88 +0,0 @@ -%%% Analysis of DG network data %%% -% This Matlab code fits the data using least-squares methods: -% "Least squares" means that the overall solution minimizes the sum of the squares of the errors made in the results of every single data point -% http://en.wikipedia.org/wiki/Least_squares -% Enter the files to import. - -% ModelDB file along with publication: -% Yim MY, Hanuschkin A, Wolfart J (2015) Hippocampus 25:297-308. -% http://onlinelibrary.wiley.com/doi/10.1002/hipo.22373/abstract -% modified and augmented by -% Man Yi Yim / 2015 -% Alexander Hanuschkin / 2011 - -close all; -set(0,'defaultaxesfontsize',10); - -offset_corrected = 4; % 0= not corrected; 1= corrected 2= corrected 3= two parameters (formula 1) 4= two parameters (formula 2) 5= three parameters -simscore1_corrected = 0; -fig1 = figure(1); - -idname = '-pp10-gaba1-kir1-st0'; -for figure_nr = 1:1 - A = importdata(strcat('sim_score',idname,'.txt')); - - % power law fit (http://en.wikipedia.org/wiki/Power_function) - fun = @(x,xdata)xdata.^x(1); % function to fit - x0 = .1; % initial values for fitting parameters - xdata = A(:,1); % x values data - ydata = A(:,2); % y values data - - switch offset_corrected - case 1 - ydata = ydata - (mean(ydata(find(xdata<0)))); - if (simscore1_corrected==1) - ydata = ydata * 1/ydata(find(A(:,1)==1)); - end - case 2 - % power law fit (http://en.wikipedia.org/wiki/Power_function) - % offset - offset = mean(ydata(find(xdata<0))); - fun = @(x,xdata)offset+(1-offset)*xdata.^x(1); % function to fit - case 3 - fun = @(x,xdata)x(1)*(1-xdata)+xdata.^x(2); - x0 = [0.1, 3]; - case 4 - fun = @(x,xdata)x(1) + (1-x(1))*xdata.^x(2); - x0 = [0.1, 1]; - case 5 - %offset = mean(ydata(find(xdata<0))); - fun = @(x,xdata)x(1)*(1-xdata).^x(2)+xdata.^x(3); - x0 = [0.1, 1, 3]; - end - - - x = lsqcurvefit(fun,x0,xdata,ydata); % - line([0 1],[0 1],'LineStyle','--','Color','blue','LineWidth',3,'MarkerSize',1); hold on; - line([0.6 0.6],[0 1],'LineStyle','-.','Color','green','LineWidth',3,'MarkerSize',1); hold on; - h1 = plot([0:0.005:1],fun(real(x),[0:0.005:1]),'-','Color','red','LineWidth',3,'MarkerSize',1);hold on; - ts = sprintf('fit x^{%.3f}',x); - h2 = plot(xdata,ydata,'o','Color',[1 0.50 0.25],'LineWidth',2,'MarkerSize',10);hold on; - xlabel('Input'); - xlim([-0.05 1.05]); - ylabel('Output'); - ylim([-0.05 1.05]); - title_str = sprintf('Similarity Score Fig. %d',figure_nr); - title(title_str); - - ts = sprintf('%.3f',x); - text(0.2,0.5,ts,'FontSize',16,'Color','red') - ts = sprintf('%.3f',fun(real(x),0.6)); - text(0.2,0.8,ts,'FontSize',16,'Color','green') - axis square -end - -switch offset_corrected - case 0 - text(-0.5,4.35,'(sim_{in})^c','FontSize',22,'Color','black'); - case 1 - text(-0.5,4.35,'(sim_{in})^c - mean subtraced','FontSize',22,'Color','black'); - case 2 - text(-1.3,4.35,'+(1-)(sim_{in})^c','FontSize',22,'Color','black'); - case 3 - text(-1.3,4.35,'a(1-sim_{in})+(sim_{in})^c','FontSize',22,'Color','black'); - case 4 - text(-1.3,4.35,'a+(1-a)(sim_{in})^c','FontSize',22,'Color','black'); - case 5 - text(-1.3,4.35,'a(1-sim_{in})^b+(sim_{in})^c','FontSize',22,'Color','black'); -end \ No newline at end of file diff --git a/GC.hoc b/GC.hoc deleted file mode 100644 index b192ea3..0000000 --- a/GC.hoc +++ /dev/null @@ -1,298 +0,0 @@ -//********************** GRANULE CELL **************************************** - -// extracted from -// Dentate gyrus network model -// Santhakumar V, Aradi I, Soltesz I (2005) J Neurophysiol 93:437-53 -// https://senselab.med.yale.edu/ModelDB/showModel.cshtml?model=51781&file=\dentategyrusnet2005\DG500_M7.hoc - -// ModelDB file along with publication: -// Yim MY, Hanuschkin A, Wolfart J (2015) Hippocampus 25:297-308. -// http://onlinelibrary.wiley.com/doi/10.1002/hipo.22373/abstract - -// modified and augmented by -// Man Yi Yim / 2015 -// Alexander Hanuschkin / 2011 - -objref Gcell[ngcell] - -begintemplate GranuleCell - -external scale_gpas_dg_ -external scale_sk_dg_ -external scale_gabaa_ -external scale_kir_ - -ndend1=4 -ndend2=4 -public pre_list, connect_pre, subsets, is_art, is_connected -public vbc2gc, vmc2gc, vhc2gc, vgc2bc, vbc2bc, vmc2bc, vhc2bc, vgc2mc, vbc2mc, vmc2mc, vhc2mc, vgc2hc, vmc2hc -public soma, gcdend1, gcdend2 -public all, gcldend, pdend, mdend, ddend - -create soma, gcdend1[ndend1], gcdend2[ndend2] -objref syn, pre_list - - -//to include steady state current injection -nst=1 -objectvar stim[nst] -public stim -// double stimdur[nst], stimdel[nst], stimamp[nst] -// public stim, stimdur, stimamp, stimdel - -proc init() { - pre_list = new List() - subsets() - gctemp() - synapse() -} - -objref all, gcldend, pdend, mdend, ddend -proc subsets(){ local i - objref all, gcldend, pdend, mdend, ddend - all = new SectionList() - soma all.append() - for i=0, 3 gcdend1 [i] all.append() - for i=0, 3 gcdend2 [i] all.append() - - gcldend = new SectionList() - gcdend1 [0] gcldend.append() - gcdend2 [0] gcldend.append() - - pdend = new SectionList() - gcdend1 [1] pdend.append() - gcdend2 [1] pdend.append() - - mdend = new SectionList() - gcdend1 [2] mdend.append() - gcdend2 [2] mdend.append() - - ddend = new SectionList() - gcdend1 [3] ddend.append() - gcdend2 [3] ddend.append() -} -proc gctemp() { - -scale_area = 1./1.13 - -// ********** Parameters for reversal potentials (assigned below) ********* -e_gabaa_ = -70. // reversal potential GABAA - -// ***************** Parameters -g_pas_fit_ = 1.44e-05 * (scale_gpas_dg_ / 100.) -gkbar_kir_fit_ = 1.44e-05 * (scale_kir_ / 100.) -ggabaabar_fit_ = 0.722e-05 * (scale_gabaa_ / 100.) - -// *********************** PAS ****************************************** -cm_fit_ = 1. -Ra_fit_ = 184. // fitted - -// *********************** KIR ***************************************** -vhalfl_kir_fit_ = -98.923594 // for Botzman I/V curve, fitted -kl_kir_fit_ = 10.888538 // for Botzman I/V curve, fitted -q10_kir_fit_ = 1. // temperature factor, set to 1 -vhalft_kir_fit_ = 67.0828 // 3 values for tau func from Stegen et al. 2011 -at_kir_fit_ = 0.00610779 -bt_kir_fit_ = 0.0817741 - -// ********************* Neuron Morphology etc *************************** -LJP_ = -10. // Liquid junction potential [mV] -V_rest = -68.16+LJP_ // resting potential [mV] -V_init = -68.16+LJP_ // initial potential [mV] - -// ******************** GABAA ******************** -e_pas_fit_ = -83.8 -e_pas_fit_Dend = -81.74 - - soma {nseg=1 L=16.8*scale_area diam=16.8*scale_area} // changed L & diam - - gcdend1 [0] {nseg=1 L=50*scale_area diam=3*scale_area} - for i = 1, 3 gcdend1 [i] {nseg=1 L=150*scale_area diam=3*scale_area} - - gcdend2 [0] {nseg=1 L=50*scale_area diam=3*scale_area} - for i = 1, 3 gcdend2 [i] {nseg=1 L=150*scale_area diam=3*scale_area} - - forsec all { - insert ccanl - catau_ccanl = 10 - caiinf_ccanl = 5.e-6 - Ra=Ra_fit_ - } - - soma {insert ichan2 - gnatbar_ichan2=0.12 // value Aradi & Holmes 1999 - gkfbar_ichan2=0.016 - gksbar_ichan2=0.006 - gl_ichan2 = g_pas_fit_ - el_ichan2 = e_pas_fit_ // set leak reversal poti to gain Vrest of cell - insert ka - gkabar_ka=0.012 - insert nca - gncabar_nca=0.002 - insert lca - glcabar_lca=0.005 - insert tca - gcatbar_tca=0.000037 - insert sk - gskbar_sk=0.001 * (scale_sk_dg_ / 100.) - insert bk - gkbar_bk=0.0006 - cm=cm_fit_ - } - - forsec gcldend {insert ichan2 - gnatbar_ichan2=0.018 // value Aradi & Holmes 1999 - gkfbar_ichan2=0.004 - gksbar_ichan2=0.006 - gl_ichan2 = g_pas_fit_ - el_ichan2 = e_pas_fit_ // set leak reversal poti to gain Vrest of cell - insert nca // HAV-N- Ca channel - gncabar_nca=0.003 // value Aradi & Holmes 1999 - insert lca - glcabar_lca=0.0075 - insert tca - gcatbar_tca=0.000075 - insert sk - gskbar_sk=0.0004 * (scale_sk_dg_ / 100.) - insert bk - gkbar_bk=0.0006 - cm=cm_fit_ - } - - forsec pdend {insert ichan2 - gnatbar_ichan2=0.013 // value Aradi & Holmes 1999 - gkfbar_ichan2=0.004 - gksbar_ichan2=0.006 - gl_ichan2 = g_pas_fit_ * (0.000063/0.00004) - el_ichan2 = e_pas_fit_Dend // see comment above - insert nca // HAV-N- Ca channel - gncabar_nca=0.001 // value Aradi & Holmes 1999 - insert lca - glcabar_lca=0.0075 - insert tca - gcatbar_tca=0.00025 - insert sk - gskbar_sk=0.0002 * (scale_sk_dg_ / 100.) - insert bk - gkbar_bk=0.001 - cm=cm_fit_*1.6 - } - - forsec mdend {insert ichan2 - gnatbar_ichan2=0.008 // value Aradi & Holmes 1999 - gkfbar_ichan2=0.001 - gksbar_ichan2=0.006 - gl_ichan2 = g_pas_fit_ * (0.000063/0.00004) - el_ichan2 = e_pas_fit_Dend // see comment above - insert nca - gncabar_nca=0.001 // value Aradi & Holmes 1999 - insert lca - glcabar_lca=0.0005 - insert tca - gcatbar_tca=0.0005 - insert sk - gskbar_sk=0.0 * (scale_sk_dg_ / 100.) - insert bk - gkbar_bk=0.0024 - cm=cm_fit_*1.6 - } - - forsec ddend {insert ichan2 - gnatbar_ichan2=0.0 - gkfbar_ichan2=0.001 - gksbar_ichan2=0.008 - gl_ichan2 = g_pas_fit_ * (0.000063/0.00004) - el_ichan2 = e_pas_fit_Dend // see comment above - insert nca - gncabar_nca=0.001 // value Aradi & Holmes 1999 - insert lca - glcabar_lca=0.0 - insert tca - gcatbar_tca=0.001 - insert sk - gskbar_sk=0.0 * (scale_sk_dg_ / 100.) - insert bk - gkbar_bk=0.0024 - cm=cm_fit_*1.6 - } - - - connect gcdend1[0](0), soma(1) - connect gcdend2[0](0), soma(1) - for i=1,3 { - connect gcdend1[i](0), gcdend1[i-1](1) - } - for i=1,3 { - connect gcdend2[i](0), gcdend2[i-1](1) - } - - forsec all { - insert kir // kir conductance added in Yim et al. 2015, note that eK=-90mV is used instead of -105mV as reported in the paper - gkbar_kir = gkbar_kir_fit_ - vhalfl_kir = vhalfl_kir_fit_ - kl_kir = kl_kir_fit_ - vhalft_kir = vhalft_kir_fit_ - at_kir = at_kir_fit_ - bt_kir = bt_kir_fit_ - ggabaa_ichan2 = ggabaabar_fit_ // added GabaA in Yim et al. 2015 - egabaa_ichan2 = e_gabaa_ // reversal potential GABAA added in Yim et al. 2015 - ena = 50 // ena was unified from enat=55 (BC, HIPP, MC) and enat=45 (GC) in Santhakumar et al. (2005) - ek = -90 // simplified ekf=eks=ek=esk; note the eK was erroneously reported as -105mV in the Yim et al. 2015 - cao_ccanl = 2 } -} // end of gctemp() - -// Retrieval of objref arguments uses the syntax: $o1, $o2, ..., $oi. -// http://web.mit.edu/neuron_v7.1/doc/help/neuron/general/ocsyntax.html#arguments -proc connect_pre() { - soma $o2 = new NetCon (&v(1), $o1) -} - - -// Define synapses on to GCs using -//- an Exp2Syn object (parameters tau1 -rise, tau2 -decay, -// time constant [ms] and e - rev potential [mV] -// delay [ms] and weight -variablr betw 0 and 1 [1 corresponding to 1 'S] - -proc synapse() { - gcdend1[3] syn = new Exp2Syn(0.5) // PP syn based on data from Greg Hollrigel and Kevin Staley NOTE: both synapses are identical! - syn.tau1 = 1.5 syn.tau2 = 5.5 syn.e = 0 - pre_list.append(syn) - - gcdend2[3] syn = new Exp2Syn(0.5) // PP syn based on Greg and Staley - syn.tau1 = 1.5 syn.tau2 = 5.5 syn.e = 0 - pre_list.append(syn) - - gcdend1[1] syn = new Exp2Syn(0.5) // MC syn *** Estimated - syn.tau1 = 1.5 syn.tau2 = 5.5 syn.e = 0 - pre_list.append(syn) - - gcdend2[1] syn = new Exp2Syn(0.5) // MC syn *** Estimated - syn.tau1 = 1.5 syn.tau2 = 5.5 syn.e = 0 - pre_list.append(syn) - - gcdend1[3] syn = new Exp2Syn(0.5) // HIPP syn based on Harney and Jones corrected for temp - syn.tau1 = 0.5 syn.tau2 = 6 syn.e = -70 - pre_list.append(syn) - - gcdend2[3] syn = new Exp2Syn(0.5) // HIPP syn based on Harney and Jones corrected for temp - syn.tau1 = 0.5 syn.tau2 = 6 syn.e = -70 - pre_list.append(syn) - - soma syn = new Exp2Syn(0.5) // BC syn based on Bartos - syn.tau1 = 0.26 syn.tau2 = 5.5 syn.e = -70 - pre_list.append(syn) - - gcdend1[1] syn = new Exp2Syn(0.5) // NOTE: SPROUTED SYNAPSE based on Molnar and Nadler - syn.tau1 = 1.5 syn.tau2 = 5.5 syn.e = 0 - pre_list.append(syn) - - gcdend2[1] syn = new Exp2Syn(0.5) // NOTE: SPROUTED SYNAPSE - syn.tau1 = 1.5 syn.tau2 = 5.5 syn.e = 0 - pre_list.append(syn) - - // Total of 7 synapses per GC 0,1 PP; 2,3 MC; 4,5 HIPP and 6 BC 7,8 Sprout -} - -func is_art() { return 0 } - -endtemplate GranuleCell diff --git a/GCinput.py b/GCinput.py deleted file mode 100644 index 6201cab..0000000 --- a/GCinput.py +++ /dev/null @@ -1,203 +0,0 @@ -### Analysis of DG network data ### -# This Python code extracts and plots the inputs to a selected GC. -# Enter the idname and cellID - -# ModelDB file along with publication: -# Yim MY, Hanuschkin A, Wolfart J (2015) Hippocampus 25:297-308. -# http://onlinelibrary.wiley.com/doi/10.1002/hipo.22373/abstract - -# modified and augmented by -# Man Yi Yim / 2015 -# Alexander Hanuschkin / 2011 - -import pylab -import numpy -import matplotlib as mpl - -# Font size in the figure -font_size = 12 # 20 -mpl.rcParams['axes.titlesize'] = font_size+10 -mpl.rcParams['xtick.labelsize'] = font_size+6 -mpl.rcParams['ytick.labelsize'] = font_size+6 -mpl.rcParams['axes.labelsize'] = font_size+8 -mpl.rcParams['legend.fontsize'] = font_size -mpl.rcParams['font.size'] = font_size+10 - -# Parameters -tstop = 200 # duration of the simulation to be analysed -is_all_Vt = 1 # 1 if the dgvtfile contains all GCs -cellID = 8 -npp = 100 # number of PP inputs -sid = 0 # PP_box_start_ -nid = 6 # number of active PPs -#idname = "-pp10-gaba1-kir1-st0" -print 'duration = ' + str(tstop) -print 'npp = ' + str(npp) -print 'idname = ' + idname - -# Import data -stiminfile = 'StimIn-'+str(sid)+'-'+str(nid)+idname+'.txt' -dgspfile = 'DGsp-'+str(sid)+'-'+str(nid)+idname+'.txt' -dgvtfile = 'DGVt-'+str(sid)+'-'+str(nid)+idname+'.txt' -con_file = 'DGNC-'+str(sid)+'-'+str(nid)+idname+'.txt' -con_file_id = 'id_DGNC.txt' - -f2 = open(stiminfile,'r') -if f2.readline() == '': - stimin = numpy.empty(0) -else: - stimin = numpy.loadtxt(stiminfile) -f3 = open(dgspfile,'r') -if f3.readline() == '': - dgsp = numpy.empty(0) -else: - dgsp = numpy.loadtxt(dgspfile) -dgvt = numpy.loadtxt(dgvtfile) -g1 = open('130408_GC_and_input.txt','w') -g2 = open('130408_GC_mp.txt','w') - -def extract_conn(): - """extract_conn() extracts connectivity and prints the global IDs for Santhakumar model by default. - Modify the number in this code if you change the network size or ordering. - ID 0-499: GC - ID 500-505: BC - ID 506-520: MC - ID 521-526: HC - ID 527-626: PP - ID >=527: BG """ - f_read = open(con_file) - f_write = open(con_file_id,'w') - for line in f_read: - if line[0] == 'G': - j = 13 - while line[j] != ']': - j += 1 - pre = int(line[12:j]) - elif line[0] == 'B': - j = 12 - while line[j] != ']': - j += 1 - pre = 500 + int(line[11:j]) - elif line[0] == 'M': - j = 11 - while line[j] != ']': - j += 1 - pre = 506 + int(line[10:j]) - elif line[0] == 'H': - j = 10 - while line[j] != ']': - j += 1 - pre = 521 + int(line[9:j]) - elif line[0] == 'N': - j = 9 #12 - while line[j] != ']': - j += 1 - if line[7] == 'B': # correlated PP - pre = 527 + int(line[11:j]) - elif line[7] == '1': # uncorrelated PP - pre = 527 + int(line[11:j]) - elif line[7] == '[': # BG - pre = 527 + npp + int(line[8:j]) - else: - print 'Please check the file: ' + con_file - j += 2 - if line[j] == 'G': - j += 13 - k = j - while line[j] != ']': - j += 1 - post = int(line[k-1:j]) - elif line[j] == 'B': - j += 12 - k = j - while line[j] != ']': - j += 1 - post = 500 + int(line[k-1:j]) - elif line[j] == 'M': - j += 11 - k = j - while line[j] != ']': - j += 1 - post = 506 + int(line[k-1:j]) - elif line[j] == 'H': - j += 10 - k = j - while line[j] != ']': - j += 1 - post = 521 + int(line[k-1:j]) - else: - print 'Please check the file: ' + con_file - f_write.write(str(pre)+'\t'+str(post)+'\n') - f_write.close() - -#extract_conn() -conn = numpy.loadtxt(con_file_id) -pylab.figure(figsize=(10,8)) -pylab.subplot(211) -ind = conn[conn[:,1]==cellID,0] -ind.sort() -j = 0 -ystring = ['G'+str(cellID)] -if dgsp.size == 2: - if dgsp[1] == cellID: - pylab.plot(dgsp[0]*numpy.ones(2),numpy.array([j-0.4,j+0.4]),'b',linewidth=2) -elif dgsp.size > 2: - if any(dgsp[:,1]==cellID): - for spk in dgsp[dgsp[:,1]==cellID,0]: - pylab.plot(spk*numpy.ones(2),numpy.array([j-0.4,j+0.4]),'b',linewidth=2) -spktime = dgsp[dgsp[:,1]==cellID,0] -for kk in range(spktime.size): - g1.write(str(spktime[kk])+'\t'+str(cellID)) - g1.write('\n') -for k in ind: - spktime = numpy.empty(0) - if k < 527: - if dgsp.size == 2 and dgsp[1] == k: - spktime = numpy.append(spktime,dgsp[0]) - elif dgsp.size > 2: - spktime = numpy.append(spktime,dgsp[dgsp[:,1]==k,0]) - if spktime.size > 0: - j += 1 - if k < 500: - color = 'b' - ystring.append('G'+str(int(k))+r'$\to$G'+str(cellID)) - elif k < 506: - color = 'g' - ystring.append('B'+str(int(k-500))+r'$\to$G'+str(cellID)) - elif k < 521: - color = 'r' - ystring.append('M'+str(int(k-506))+r'$\to$G'+str(cellID)) - elif k < 527: - color = 'c' - ystring.append('H'+str(int(k-521))+r'$\to$G'+str(cellID)) - for m in spktime: - pylab.plot(m*numpy.ones(2),numpy.array([j-0.4,j+0.4]),color,linewidth=2) - else: - if stimin.size == 2 and stimin[1] == k-527: - spktime = numpy.append(spktime,stimin[0]) - elif stimin.size > 2: - spktime = numpy.append(spktime,stimin[stimin[:,1]==k-527,0]) - if spktime.size > 0: - j += 1 - ystring.append('PP'+str(int(k-527))+r'$\to$G'+str(cellID)) - for m in spktime: - pylab.plot(m*numpy.ones(2),numpy.array([j-0.4,j+0.4]),'k',linewidth=2) - for kk in range(spktime.size): - g1.write(str(spktime[kk])+'\t'+str(int(k))) - g1.write('\n') -pylab.yticks(numpy.arange(0.,j+1,1),ystring) -pylab.axis([0,tstop,-0.5,j+0.5]) - -pylab.subplot(212) -pylab.plot(dgvt[:,0],dgvt[:,cellID+1],'b',linewidth=2) -for kk in range(dgvt[:,0].size): - g2.write(str(dgvt[kk,0])+'\t'+str(dgvt[kk,cellID+1])) - g2.write('\n') -g1.close() -g2.close() -pylab.ylabel('Voltage (mV)') -pylab.xlabel('Time (ms)') -pylab.xlim(0,tstop) -pylab.subplots_adjust(bottom=0.08,left=0.145) -pylab.savefig('DG_Fig1C.eps') -pylab.show() diff --git a/README.md b/README.md new file mode 100644 index 0000000..8e907e5 --- /dev/null +++ b/README.md @@ -0,0 +1,62 @@ +# Biophysical Model of Dentate Gyrus + +This is a reimplementation of the dentate gyrus model from which the following papers are based: + +1. Santhakumar et al. (2005) +2. Yim et al. (2015) + +Our model is modified primarily from the Yim et al. (2005) implementation. Changes are listed at the end of this file. + +## Contributing + +Please see the `CONTRIBUTING.md` file. + +## Usage + +Scripts were tested on Windows 10 and Arch Linux. + +__I recommend using Linux, and using the shell scripts, since this will create the necessary results/output directories.__ + +### Unix machines + +``` bash +sh run.sh + +``` +### Windows + +Note that here you will need to create the directories to store results. See `run.sh` for the directory names. + +```powershell +nrnivmodl mods +.\main.hoc +``` + +### Writing a Script + +```hoc +...[load templates here (see main.hoc) for example ]... + +strdef netlabel +netlabel = "MyDG" + +random_seed = 1 + +objref dg +dg = new DentateGyrus(netlabel, random_seed) +dg.run() + +``` + +## Changes + +Here we do not list all specific changes implemented, but rather a high-level overview of how the Yim et al. (2015) code was modified. + +1. Network is now encapsulated as a `template` to facilitate batch simulations +2. NEURON code was (largely) rewritten to achieve several goals + - Making code more modular, packing routines into smaller functions + - Making variable names more expressive (still work to do here) + - Attempting to maximize local scope for variables, rather than allowing for global definitions +3. Development of functions to "transform" network in order to test different effects +4. Analysis code written in Julia +5. Added capability for spontaneous action potential discharge of granule cells diff --git a/analysis.jl b/analysis.jl new file mode 100644 index 0000000..a873a22 --- /dev/null +++ b/analysis.jl @@ -0,0 +1,88 @@ + +#################################################################### +# SCRIPT TO RUN ANALYSIS +#################################################################### +include("utilities.jl"); +default(show=false) + +# HYPERPARAMS +n_runs = 10 +patterns = 0:12 +labels = ["HC", "LR", "NR"] +fig_ext = ".png" + +# CREATE NECESSARY DIRECTORIES +create_directories_if_not_exist() + +# IDENTIFY NEURON POPULATION RANGES +populations = Dict( + "DG" => [0, 500], + "BC" => [500,506], + "MC" => [506, 521], + "HIPP" => [521, 527] +) + +for run_ ∈ 1:n_runs + for i ∈ 1:length(labels) + spikes = load_spike_files(patterns, labels[i]*"-$run_", populations) + + # CREATE RASTER PLOTS + for p ∈ unique(spikes.Pattern) + stimin = spikes[(spikes.Population .== "PP") .& (spikes.Pattern .== p), :] + plots = [] + append!(plots, [raster_plot(stimin; ylab="PP")]) + for pop ∈ keys(populations) + lb, ub = populations[pop] + popspikes = spikes[(spikes.Population .== pop) .& (spikes.Pattern .== p),:] + append!(plots, [raster_plot(popspikes; xlab="", ylab=pop)]) + end + fig = plot(reverse(plots)..., layout=grid(5, 1, heights=[0.15, 0.15, 0.15, 0.4, 0.15]), size=(400, 500)) + savefig(fig, "figures/raster-plots/raster-"*string(p)*"-"*labels[i]*"-$run_"*fig_ext) + end + end +end + + +# PATTERN SEPARATION CURVES +colors=[:blue, :red, :green] +global psfig = plot([0;1], [0;1], ls=:dash, c=:black, + xlabel="Input Correlation "*L"(r_{in})", + ylabel="Output Correlation "*L"(r_{out})", + size=(300, 300), + label=nothing, legend=:topleft) + +psc = Dict("HC"=>[], "LR"=>[], "NR"=>[]) +for i ∈ 1:length(labels) + for run ∈ 1:n_runs + spikes = load_spike_files(patterns, labels[i]*"-$run", populations) + + out = pattern_separation_curve(spikes, 100, 500) + x, y = out[:,"Input Correlation"], out[:, "Output Correlation"] + + # Remove NaNs before fitting + idx_ = (.!isnan.(x) .& .!isnan.(y)) + x = x[idx_] + y = y[idx_] + + f = fit_power_law(x, y) + append!(psc[labels[i]], f(0.6)) + + if (run == n_runs) + psm = round(mean(psc[labels[i]]), digits=2) + psse = std(psc[labels[i]])/sqrt(n_runs) + pslci = round(psm - 1.96*psse, digits=2) + psuci = round(psm + 1.96*psse, digits=2) + if n_runs > 1 + psc_label = labels[i]*" (PS="*string(psm)*" ["*string(pslci)*", "*string(psuci)*"])" + else + psc_label = labels[i]*" (PS="*string(psm) + end + else + psc_label = nothing + end + global psfig = scatter!(x, y, c=colors[i], alpha=1/(2*n_runs), label=nothing) + global psfig = plot!(0:0.01:1, x -> f(x), c=colors[i], label=psc_label) + end +end +psfig +savefig(psfig, "figures/pattern-separation/pattern-separation-curve"*fig_ext) diff --git a/analyze-spike-properties.jl b/analyze-spike-properties.jl new file mode 100644 index 0000000..104225a --- /dev/null +++ b/analyze-spike-properties.jl @@ -0,0 +1,126 @@ +using CSV, DataFrames, Plots, LaTeXStrings, Peaks, Glob + + +V = DataFrame() +for i ∈ 0:2 + v0 = CSV.read("data/fi-curves/vm-$i.000000.txt", delim="\t", header=0, DataFrame) + v0[:,"ID"] .= i + V = [V; v0] +end + + +pks, vals = findmaxima(V[V.ID .== 0,20]) +pks = pks[vals .> 0.0] +ndrs, vals = findminima(V[V.ID .== 0, 20]) + +v0 = V[V.ID .== 0, :] +plot(v0[:,1], v0[:,20]) +scatter!(v0[pks,1], v0[pks,20]) +scatter!(v0[ndrs,1], v0[ndrs,20]) + + + +""" + find_spikes(v;thresh) + +Finds the spike peaks in a voltage tracing. Spikes are defined as those maxima that exceed some value `thresh` where default is that `thresh=0.0` +""" +function find_spikes(v::Vector{Float64};thresh::Float64=0.0) + +end + + +""" + mean_frequency(v) + +Mean frequency calculated as number of action potentials during +stimulation divided by time between stimulus onset and last spike in +Hz. +""" +function mean_frequency(v::Vector{Float64}) +end + +""" + isi_log_slope(v) + +Slope of loglog interspike intervals (ISI). +""" +function isi_log_slope(v::Vector{Float64}) +end + +""" + adaptation_index2(v) + +Normalized average difference of two consecutive ISI starting from +second ISI. +""" +function adaptation_index2(v::Vector{Float64}) +end + +""" + time_to_first_spike(v) + +Time from stimulus onset to peak time of first spike in ms. +""" +function time_to_first_spike(v::Vector{Float64}) +end + +""" + time_to_last_spike(v) + +Time from stimulus onset to peak time of last spike in ms. +""" +function time_to_last_spike(v::Vector{Float64}) +end + + +""" + AP_width(v) + +Mean of width at -20 mV of action potential (AP) in ms. Mean for all AP. +""" +function AP_width(v::Vector{Float64}) +end + + +""" + AP_height(v) + +Height at peak of action potential in mV. Mean for all AP. +""" +function AP_height(v::Vector{Float64}) +end + +""" + min_voltage_between_spikes(v) + +Minimum voltage between two action potentials in mV. Mean for all ISI. +""" +function min_voltage_between_spikes(v::Vector{Float64}) +end + +""" + steady_state_voltage_stimend(v) + +The average voltage during the last 90% of the stimulus duration in mV. +""" +function steady_state_voltage_stimend(v::Vector{Float64}) +end + +""" + voltage_base(v) + +The average voltage during the last 90% before stimulus onset in mV. +""" +function voltage_base(v::Vector{Float64}) +end + +""" + voltage_after_stim(v) + +The average voltage between 25% and 75% between end of stimulus and end of recording in mV. +""" +function voltage_after_stim(v::Vector{Float64}) +end + + \ No newline at end of file diff --git a/clamps.py b/clamps.py new file mode 100644 index 0000000..663e9fa --- /dev/null +++ b/clamps.py @@ -0,0 +1,184 @@ +# ======================================================================== +# DEFINES ELECTRODES FOR EVALUATING NEURONAL PROPERTIES +# +# ======================================================================== +import numpy as np +from netpyne import specs, sim + + +class VClamp(object): + def __init__(self, + cell, + delay=100, + duration=400, + T=600, + dt=0.025, + record_step=0.1, + verbose=False): + """ Defines a voltage clamp object for stimulating and recording at the soma + + Arguments: + cell: `dict`. Cellular properties specified in NetPyNE dictionary syntax + delay: `float`. Delay until voltage step is taken + duration: `float`. Duration of current injection [ms] + T: `float`. Total duration of simulation [ms] + dt: `float`. Integration timestep [ms] + record_step: `float`. Step size at which to save data [mS] + + TODO: + - [ ] Allow for changing the stimulation/recording location + """ + self.cell = cell + self.delay = delay + self.duration = duration + self.T = T + self.dt = dt + self.record_step = record_step + self.verbose = verbose + + self.netparams = specs.NetParams() + self._set_netparams_neuron() + self._set_netparams_stim() + self._set_simparams() + + def _set_netparams_neuron(self): + self.netparams.cellParams['neuron'] = self.cell + self.netparams.popParams['pop'] = {'cellType': 'neuron', 'numCells': 1} + + def _set_netparams_stim(self): + self.netparams.stimSourceParams['vclamp'] = { + 'type': 'VClamp', + 'dur': [ + self.delay, + self.duration, + self.T - (self.delay + self.duration), + ], + } + self.netparams.stimTargetParams['vclamp->neuron'] = { + 'source': 'vclamp', + 'sec': 'soma', + 'loc': 0.5, + 'conds': {'pop': 'pop', 'cellList': [0]}, + } + + def _set_simparams(self): + """ + TODO: + - [ ] Make it so that you can record many other types of currents + """ + self.simconfig = specs.SimConfig() + self.simconfig.duration = self.T + self.simconfig.dt = self.dt + self.simconfig.verbose = self.verbose + self.simconfig.recordCells = ["all"] + self.simconfig.recordTraces = { + 'i_na': {'sec': 'soma', 'loc': 0.5, 'var': 'ina'}, + 'i_k': {'sec': 'soma', 'loc': 0.5, 'var': 'ik'}, + } + self.simconfig.recordStep = self.record_step + + def __call__(self, amp): + """ + Arguments: + amp: `list` of `float`. Voltage at which membrane is to be maintained [nA] + + Returns: + `dict`. Simulation data with the following key-value pairs + - `t`: List representing time [ms] + - `V`: List representing membrane voltage [mV] + - `spkt`: List of spike times [ms] + - `avg_rate`: Average firing rate across whole recording [Hz] + - `rate`: Firing rate only during current injection [Hz] + """ + self.netparams.stimSourceParams['vclamp']['amp'] = amp + sim.createSimulateAnalyze(self.netparams, self.simconfig) + results = { + 't': sim.allSimData['t'], + 'i_na': np.array(sim.allSimData['i_na']['cell_0']), + 'i_k': np.array(sim.allSimData['i_k']['cell_0']), + } + return results + + +class IClamp(object): + def __init__(self, cell, delay=100, duration=200, T=400, dt=0.025, + record_step=0.1, verbose=False): + """ Runs a current-clamp experiment stimulating and recording at the soma + + Arguments: + cell: `dict`. Cellular properties specified in NetPyNE dictionary syntax + delay: `float`. Time after which current starts [ms] + duration: `float`. Duration of current injection [ms] + T: `float`. Total duration of simulation [ms] + dt: `float`. Integration timestep [ms] + record_step: `float`. Step size at which to save data [mS] + + TODO: + - [ ] Allow for changing the stimulation/recording location + """ + self.cell = cell + self.delay = delay + self.duration = duration + self.T = T + self.dt = dt + self.record_step = record_step + self.verbose = verbose + + self.netparams = specs.NetParams() + self._set_netparams_neuron() + self._set_netparams_stim() + self._set_simparams() + + def _set_netparams_neuron(self): + self.netparams.cellParams['neuron'] = self.cell + self.netparams.popParams['pop'] = {'cellType': 'neuron', 'numCells': 1} + + def _set_netparams_stim(self): + self.netparams.stimSourceParams['iclamp'] = {'type': 'IClamp', + 'del': self.delay, + 'dur': self.duration} + self.netparams.stimTargetParams['iclamp->neuron'] = { + 'source': 'iclamp', + 'sec': 'soma', + 'loc': 0.5, + 'conds': {'pop': 'pop', 'cellList': [0]}, + } + + def _set_simparams(self): + self.simconfig = specs.SimConfig() + self.simconfig.duration = self.T + self.simconfig.dt = self.dt + self.simconfig.verbose = self.verbose + self.simconfig.recordCells = ["all"] + self.simconfig.recordTraces = { + 'V_soma': {'sec': 'soma', 'loc': 0.5, 'var': 'v'}, + 'iclamp': {'sec': 'soma', 'loc': 0.5, 'stim': 'iclamp->neuron', 'var': 'i'} + } + self.simconfig.recordStep = self.record_step + + def __call__(self, amp): + """ + Arguments: + amp: `float`. Current to be injected [nA] + + Returns: + `dict`. Simulation data with the following key-value pairs + - `t`: List representing time [ms] + - `V`: List representing membrane voltage [mV] + - `spkt`: List of spike times [ms] + - `avg_rate`: Average firing rate across whole recording [Hz] + - `rate`: Firing rate only during current injection [Hz] + """ + self.netparams.stimSourceParams['iclamp']['amp'] = amp + sim.createSimulateAnalyze(self.netparams, self.simconfig) + + print(list(sim.allSimData.keys())) + results = { + 't': np.array(sim.allSimData['t']), + 'V': np.array(sim.allSimData['V_soma']['cell_0']), + 'spkt': np.array(sim.allSimData['spkt']), + 'avg_rate': sim.allSimData['avgRate'], + 'rate': len(sim.allSimData['spkt']) / (self.duration / 1000), + 'i': sim.allSimData['iclamp']['cell_0'] + } + return results diff --git a/cleanup.sh b/cleanup.sh new file mode 100644 index 0000000..f995c28 --- /dev/null +++ b/cleanup.sh @@ -0,0 +1,3 @@ +rm -R x86_64/ +rm -R data +rm -R figures diff --git a/ficurve.txt b/ficurve.txt new file mode 100644 index 0000000..84e12df --- /dev/null +++ b/ficurve.txt @@ -0,0 +1,366 @@ +0.000000 0.000000 1.000000 1.000000 1.000000 0.000400 0.000000 +0.020000 0.000000 1.000000 1.000000 1.000000 0.000400 0.000000 +0.040000 0.000000 1.000000 1.000000 1.000000 0.000400 0.000000 +0.060000 2.105263 1.000000 1.000000 1.000000 0.000400 0.000000 +0.080000 3.157895 1.000000 1.000000 1.000000 0.000400 0.000000 +0.100000 4.210526 1.000000 1.000000 1.000000 0.000400 0.000000 +0.120000 4.736842 1.000000 1.000000 1.000000 0.000400 0.000000 +0.140000 5.789474 1.000000 1.000000 1.000000 0.000400 0.000000 +0.160000 6.315789 1.000000 1.000000 1.000000 0.000400 0.000000 +0.180000 6.842105 1.000000 1.000000 1.000000 0.000400 0.000000 +0.200000 7.894737 1.000000 1.000000 1.000000 0.000400 0.000000 +0.220000 8.421053 1.000000 1.000000 1.000000 0.000400 0.000000 +0.240000 9.473684 1.000000 1.000000 1.000000 0.000400 0.000000 +0.260000 10.526316 1.000000 1.000000 1.000000 0.000400 0.000000 +0.280000 11.578947 1.000000 1.000000 1.000000 0.000400 0.000000 +0.300000 13.157895 1.000000 1.000000 1.000000 0.000400 0.000000 +0.320000 14.210526 1.000000 1.000000 1.000000 0.000400 0.000000 +0.340000 15.789474 1.000000 1.000000 1.000000 0.000400 0.000000 +0.360000 17.894737 1.000000 1.000000 1.000000 0.000400 0.000000 +0.380000 19.473684 1.000000 1.000000 1.000000 0.000400 0.000000 +0.400000 21.578947 1.000000 1.000000 1.000000 0.000400 0.000000 +0.000000 0.000000 1.030303 1.150000 1.200000 0.000800 0.000000 +0.020000 0.000000 1.030303 1.150000 1.200000 0.000800 0.000000 +0.040000 0.000000 1.030303 1.150000 1.200000 0.000800 0.000000 +0.060000 2.105263 1.030303 1.150000 1.200000 0.000800 0.000000 +0.080000 3.157895 1.030303 1.150000 1.200000 0.000800 0.000000 +0.100000 4.210526 1.030303 1.150000 1.200000 0.000800 0.000000 +0.120000 4.736842 1.030303 1.150000 1.200000 0.000800 0.000000 +0.140000 5.789474 1.030303 1.150000 1.200000 0.000800 0.000000 +0.160000 6.315789 1.030303 1.150000 1.200000 0.000800 0.000000 +0.180000 7.368421 1.030303 1.150000 1.200000 0.000800 0.000000 +0.200000 7.894737 1.030303 1.150000 1.200000 0.000800 0.000000 +0.220000 8.947368 1.030303 1.150000 1.200000 0.000800 0.000000 +0.240000 10.000000 1.030303 1.150000 1.200000 0.000800 0.000000 +0.260000 11.052632 1.030303 1.150000 1.200000 0.000800 0.000000 +0.280000 12.105263 1.030303 1.150000 1.200000 0.000800 0.000000 +0.300000 13.157895 1.030303 1.150000 1.200000 0.000800 0.000000 +0.320000 14.736842 1.030303 1.150000 1.200000 0.000800 0.000000 +0.340000 16.315789 1.030303 1.150000 1.200000 0.000800 0.000000 +0.360000 17.894737 1.030303 1.150000 1.200000 0.000800 0.000000 +0.380000 20.000000 1.030303 1.150000 1.200000 0.000800 0.000000 +0.400000 22.105263 1.030303 1.150000 1.200000 0.000800 0.000000 +0.000000 0.000000 0.787879 1.300000 1.300000 0.000400 0.000000 +0.020000 0.000000 0.787879 1.300000 1.300000 0.000400 0.000000 +0.040000 0.000000 0.787879 1.300000 1.300000 0.000400 0.000000 +0.060000 2.105263 0.787879 1.300000 1.300000 0.000400 0.000000 +0.080000 4.210526 0.787879 1.300000 1.300000 0.000400 0.000000 +0.100000 5.263158 0.787879 1.300000 1.300000 0.000400 0.000000 +0.120000 6.315789 0.787879 1.300000 1.300000 0.000400 0.000000 +0.140000 7.368421 0.787879 1.300000 1.300000 0.000400 0.000000 +0.160000 8.421053 0.787879 1.300000 1.300000 0.000400 0.000000 +0.180000 9.473684 0.787879 1.300000 1.300000 0.000400 0.000000 +0.200000 11.052632 0.787879 1.300000 1.300000 0.000400 0.000000 +0.220000 12.631579 0.787879 1.300000 1.300000 0.000400 0.000000 +0.240000 11.052632 0.787879 1.300000 1.300000 0.000400 0.000000 +0.260000 15.789474 0.787879 1.300000 1.300000 0.000400 0.000000 +0.280000 17.368421 0.787879 1.300000 1.300000 0.000400 0.000000 +0.300000 18.947368 0.787879 1.300000 1.300000 0.000400 0.000000 +0.320000 21.052632 0.787879 1.300000 1.300000 0.000400 0.000000 +0.340000 23.157895 0.787879 1.300000 1.300000 0.000400 0.000000 +0.360000 25.263158 0.787879 1.300000 1.300000 0.000400 0.000000 +0.380000 27.368421 0.787879 1.300000 1.300000 0.000400 0.000000 +0.400000 29.473684 0.787879 1.300000 1.300000 0.000400 0.000000 +0.000000 0.000000 1.000000 1.000000 1.000000 0.000400 0.000000 +0.004000 0.000000 1.000000 1.000000 1.000000 0.000400 0.000000 +0.008000 0.000000 1.000000 1.000000 1.000000 0.000400 0.000000 +0.012000 0.000000 1.000000 1.000000 1.000000 0.000400 0.000000 +0.016000 0.000000 1.000000 1.000000 1.000000 0.000400 0.000000 +0.020000 0.000000 1.000000 1.000000 1.000000 0.000400 0.000000 +0.024000 0.000000 1.000000 1.000000 1.000000 0.000400 0.000000 +0.028000 0.000000 1.000000 1.000000 1.000000 0.000400 0.000000 +0.032000 0.000000 1.000000 1.000000 1.000000 0.000400 0.000000 +0.036000 0.000000 1.000000 1.000000 1.000000 0.000400 0.000000 +0.040000 0.000000 1.000000 1.000000 1.000000 0.000400 0.000000 +0.044000 0.000000 1.000000 1.000000 1.000000 0.000400 0.000000 +0.048000 1.333333 1.000000 1.000000 1.000000 0.000400 0.000000 +0.052000 1.333333 1.000000 1.000000 1.000000 0.000400 0.000000 +0.056000 2.000000 1.000000 1.000000 1.000000 0.000400 0.000000 +0.060000 2.000000 1.000000 1.000000 1.000000 0.000400 0.000000 +0.064000 2.666667 1.000000 1.000000 1.000000 0.000400 0.000000 +0.068000 2.666667 1.000000 1.000000 1.000000 0.000400 0.000000 +0.072000 3.333333 1.000000 1.000000 1.000000 0.000400 0.000000 +0.076000 3.333333 1.000000 1.000000 1.000000 0.000400 0.000000 +0.080000 3.333333 1.000000 1.000000 1.000000 0.000400 0.000000 +0.084000 3.333333 1.000000 1.000000 1.000000 0.000400 0.000000 +0.088000 4.000000 1.000000 1.000000 1.000000 0.000400 0.000000 +0.092000 4.000000 1.000000 1.000000 1.000000 0.000400 0.000000 +0.096000 4.000000 1.000000 1.000000 1.000000 0.000400 0.000000 +0.100000 4.000000 1.000000 1.000000 1.000000 0.000400 0.000000 +0.104000 4.666667 1.000000 1.000000 1.000000 0.000400 0.000000 +0.108000 4.666667 1.000000 1.000000 1.000000 0.000400 0.000000 +0.112000 4.666667 1.000000 1.000000 1.000000 0.000400 0.000000 +0.116000 4.666667 1.000000 1.000000 1.000000 0.000400 0.000000 +0.120000 4.666667 1.000000 1.000000 1.000000 0.000400 0.000000 +0.124000 5.333333 1.000000 1.000000 1.000000 0.000400 0.000000 +0.128000 5.333333 1.000000 1.000000 1.000000 0.000400 0.000000 +0.132000 5.333333 1.000000 1.000000 1.000000 0.000400 0.000000 +0.136000 5.333333 1.000000 1.000000 1.000000 0.000400 0.000000 +0.140000 5.333333 1.000000 1.000000 1.000000 0.000400 0.000000 +0.144000 6.000000 1.000000 1.000000 1.000000 0.000400 0.000000 +0.148000 6.000000 1.000000 1.000000 1.000000 0.000400 0.000000 +0.152000 6.000000 1.000000 1.000000 1.000000 0.000400 0.000000 +0.156000 6.000000 1.000000 1.000000 1.000000 0.000400 0.000000 +0.160000 6.000000 1.000000 1.000000 1.000000 0.000400 0.000000 +0.164000 6.666667 1.000000 1.000000 1.000000 0.000400 0.000000 +0.168000 6.666667 1.000000 1.000000 1.000000 0.000400 0.000000 +0.172000 6.666667 1.000000 1.000000 1.000000 0.000400 0.000000 +0.176000 6.666667 1.000000 1.000000 1.000000 0.000400 0.000000 +0.180000 7.333333 1.000000 1.000000 1.000000 0.000400 0.000000 +0.184000 7.333333 1.000000 1.000000 1.000000 0.000400 0.000000 +0.188000 7.333333 1.000000 1.000000 1.000000 0.000400 0.000000 +0.192000 7.333333 1.000000 1.000000 1.000000 0.000400 0.000000 +0.196000 7.333333 1.000000 1.000000 1.000000 0.000400 0.000000 +0.200000 8.000000 1.000000 1.000000 1.000000 0.000400 0.000000 +0.204000 8.000000 1.000000 1.000000 1.000000 0.000400 0.000000 +0.208000 8.000000 1.000000 1.000000 1.000000 0.000400 0.000000 +0.212000 8.000000 1.000000 1.000000 1.000000 0.000400 0.000000 +0.216000 8.666667 1.000000 1.000000 1.000000 0.000400 0.000000 +0.220000 8.666667 1.000000 1.000000 1.000000 0.000400 0.000000 +0.224000 8.666667 1.000000 1.000000 1.000000 0.000400 0.000000 +0.228000 8.666667 1.000000 1.000000 1.000000 0.000400 0.000000 +0.232000 8.666667 1.000000 1.000000 1.000000 0.000400 0.000000 +0.236000 9.333333 1.000000 1.000000 1.000000 0.000400 0.000000 +0.240000 9.333333 1.000000 1.000000 1.000000 0.000400 0.000000 +0.244000 9.333333 1.000000 1.000000 1.000000 0.000400 0.000000 +0.248000 10.000000 1.000000 1.000000 1.000000 0.000400 0.000000 +0.252000 10.000000 1.000000 1.000000 1.000000 0.000400 0.000000 +0.256000 10.000000 1.000000 1.000000 1.000000 0.000400 0.000000 +0.260000 10.000000 1.000000 1.000000 1.000000 0.000400 0.000000 +0.264000 10.666667 1.000000 1.000000 1.000000 0.000400 0.000000 +0.268000 10.666667 1.000000 1.000000 1.000000 0.000400 0.000000 +0.272000 10.666667 1.000000 1.000000 1.000000 0.000400 0.000000 +0.276000 11.333333 1.000000 1.000000 1.000000 0.000400 0.000000 +0.280000 11.333333 1.000000 1.000000 1.000000 0.000400 0.000000 +0.284000 11.333333 1.000000 1.000000 1.000000 0.000400 0.000000 +0.288000 12.000000 1.000000 1.000000 1.000000 0.000400 0.000000 +0.292000 12.000000 1.000000 1.000000 1.000000 0.000400 0.000000 +0.296000 12.000000 1.000000 1.000000 1.000000 0.000400 0.000000 +0.300000 12.666667 1.000000 1.000000 1.000000 0.000400 0.000000 +0.304000 12.666667 1.000000 1.000000 1.000000 0.000400 0.000000 +0.308000 12.666667 1.000000 1.000000 1.000000 0.000400 0.000000 +0.312000 13.333333 1.000000 1.000000 1.000000 0.000400 0.000000 +0.316000 13.333333 1.000000 1.000000 1.000000 0.000400 0.000000 +0.320000 14.000000 1.000000 1.000000 1.000000 0.000400 0.000000 +0.324000 14.000000 1.000000 1.000000 1.000000 0.000400 0.000000 +0.328000 14.666667 1.000000 1.000000 1.000000 0.000400 0.000000 +0.332000 14.666667 1.000000 1.000000 1.000000 0.000400 0.000000 +0.336000 14.666667 1.000000 1.000000 1.000000 0.000400 0.000000 +0.340000 15.333333 1.000000 1.000000 1.000000 0.000400 0.000000 +0.344000 15.333333 1.000000 1.000000 1.000000 0.000400 0.000000 +0.348000 16.000000 1.000000 1.000000 1.000000 0.000400 0.000000 +0.352000 16.000000 1.000000 1.000000 1.000000 0.000400 0.000000 +0.356000 16.666667 1.000000 1.000000 1.000000 0.000400 0.000000 +0.360000 16.666667 1.000000 1.000000 1.000000 0.000400 0.000000 +0.364000 17.333333 1.000000 1.000000 1.000000 0.000400 0.000000 +0.368000 17.333333 1.000000 1.000000 1.000000 0.000400 0.000000 +0.372000 18.000000 1.000000 1.000000 1.000000 0.000400 0.000000 +0.376000 18.666667 1.000000 1.000000 1.000000 0.000400 0.000000 +0.380000 18.666667 1.000000 1.000000 1.000000 0.000400 0.000000 +0.384000 19.333333 1.000000 1.000000 1.000000 0.000400 0.000000 +0.388000 19.333333 1.000000 1.000000 1.000000 0.000400 0.000000 +0.392000 20.000000 1.000000 1.000000 1.000000 0.000400 0.000000 +0.396000 20.000000 1.000000 1.000000 1.000000 0.000400 0.000000 +0.400000 20.666667 1.000000 1.000000 1.000000 0.000400 0.000000 +0.000000 0.000000 1.030303 1.150000 1.200000 0.000800 0.000000 +0.004000 0.000000 1.030303 1.150000 1.200000 0.000800 0.000000 +0.008000 0.000000 1.030303 1.150000 1.200000 0.000800 0.000000 +0.012000 0.000000 1.030303 1.150000 1.200000 0.000800 0.000000 +0.016000 0.000000 1.030303 1.150000 1.200000 0.000800 0.000000 +0.020000 0.000000 1.030303 1.150000 1.200000 0.000800 0.000000 +0.024000 0.000000 1.030303 1.150000 1.200000 0.000800 0.000000 +0.028000 0.000000 1.030303 1.150000 1.200000 0.000800 0.000000 +0.032000 0.000000 1.030303 1.150000 1.200000 0.000800 0.000000 +0.036000 0.000000 1.030303 1.150000 1.200000 0.000800 0.000000 +0.040000 0.000000 1.030303 1.150000 1.200000 0.000800 0.000000 +0.044000 0.000000 1.030303 1.150000 1.200000 0.000800 0.000000 +0.048000 0.000000 1.030303 1.150000 1.200000 0.000800 0.000000 +0.052000 0.666667 1.030303 1.150000 1.200000 0.000800 0.000000 +0.056000 1.333333 1.030303 1.150000 1.200000 0.000800 0.000000 +0.060000 2.000000 1.030303 1.150000 1.200000 0.000800 0.000000 +0.064000 2.000000 1.030303 1.150000 1.200000 0.000800 0.000000 +0.068000 2.666667 1.030303 1.150000 1.200000 0.000800 0.000000 +0.072000 2.666667 1.030303 1.150000 1.200000 0.000800 0.000000 +0.076000 3.333333 1.030303 1.150000 1.200000 0.000800 0.000000 +0.080000 3.333333 1.030303 1.150000 1.200000 0.000800 0.000000 +0.084000 3.333333 1.030303 1.150000 1.200000 0.000800 0.000000 +0.088000 4.000000 1.030303 1.150000 1.200000 0.000800 0.000000 +0.092000 4.000000 1.030303 1.150000 1.200000 0.000800 0.000000 +0.096000 4.000000 1.030303 1.150000 1.200000 0.000800 0.000000 +0.100000 4.000000 1.030303 1.150000 1.200000 0.000800 0.000000 +0.104000 4.666667 1.030303 1.150000 1.200000 0.000800 0.000000 +0.108000 4.666667 1.030303 1.150000 1.200000 0.000800 0.000000 +0.112000 4.666667 1.030303 1.150000 1.200000 0.000800 0.000000 +0.116000 4.666667 1.030303 1.150000 1.200000 0.000800 0.000000 +0.120000 4.666667 1.030303 1.150000 1.200000 0.000800 0.000000 +0.124000 5.333333 1.030303 1.150000 1.200000 0.000800 0.000000 +0.128000 5.333333 1.030303 1.150000 1.200000 0.000800 0.000000 +0.132000 5.333333 1.030303 1.150000 1.200000 0.000800 0.000000 +0.136000 5.333333 1.030303 1.150000 1.200000 0.000800 0.000000 +0.140000 6.000000 1.030303 1.150000 1.200000 0.000800 0.000000 +0.144000 6.000000 1.030303 1.150000 1.200000 0.000800 0.000000 +0.148000 6.000000 1.030303 1.150000 1.200000 0.000800 0.000000 +0.152000 6.000000 1.030303 1.150000 1.200000 0.000800 0.000000 +0.156000 6.000000 1.030303 1.150000 1.200000 0.000800 0.000000 +0.160000 6.666667 1.030303 1.150000 1.200000 0.000800 0.000000 +0.164000 6.666667 1.030303 1.150000 1.200000 0.000800 0.000000 +0.168000 6.666667 1.030303 1.150000 1.200000 0.000800 0.000000 +0.172000 6.666667 1.030303 1.150000 1.200000 0.000800 0.000000 +0.176000 6.666667 1.030303 1.150000 1.200000 0.000800 0.000000 +0.180000 7.333333 1.030303 1.150000 1.200000 0.000800 0.000000 +0.184000 7.333333 1.030303 1.150000 1.200000 0.000800 0.000000 +0.188000 7.333333 1.030303 1.150000 1.200000 0.000800 0.000000 +0.192000 7.333333 1.030303 1.150000 1.200000 0.000800 0.000000 +0.196000 8.000000 1.030303 1.150000 1.200000 0.000800 0.000000 +0.200000 8.000000 1.030303 1.150000 1.200000 0.000800 0.000000 +0.204000 8.000000 1.030303 1.150000 1.200000 0.000800 0.000000 +0.208000 8.000000 1.030303 1.150000 1.200000 0.000800 0.000000 +0.212000 8.666667 1.030303 1.150000 1.200000 0.000800 0.000000 +0.216000 8.666667 1.030303 1.150000 1.200000 0.000800 0.000000 +0.220000 8.666667 1.030303 1.150000 1.200000 0.000800 0.000000 +0.224000 8.666667 1.030303 1.150000 1.200000 0.000800 0.000000 +0.228000 9.333333 1.030303 1.150000 1.200000 0.000800 0.000000 +0.232000 9.333333 1.030303 1.150000 1.200000 0.000800 0.000000 +0.236000 9.333333 1.030303 1.150000 1.200000 0.000800 0.000000 +0.240000 9.333333 1.030303 1.150000 1.200000 0.000800 0.000000 +0.244000 10.000000 1.030303 1.150000 1.200000 0.000800 0.000000 +0.248000 10.000000 1.030303 1.150000 1.200000 0.000800 0.000000 +0.252000 10.000000 1.030303 1.150000 1.200000 0.000800 0.000000 +0.256000 10.666667 1.030303 1.150000 1.200000 0.000800 0.000000 +0.260000 10.666667 1.030303 1.150000 1.200000 0.000800 0.000000 +0.264000 10.666667 1.030303 1.150000 1.200000 0.000800 0.000000 +0.268000 11.333333 1.030303 1.150000 1.200000 0.000800 0.000000 +0.272000 11.333333 1.030303 1.150000 1.200000 0.000800 0.000000 +0.276000 11.333333 1.030303 1.150000 1.200000 0.000800 0.000000 +0.280000 12.000000 1.030303 1.150000 1.200000 0.000800 0.000000 +0.284000 12.000000 1.030303 1.150000 1.200000 0.000800 0.000000 +0.288000 12.000000 1.030303 1.150000 1.200000 0.000800 0.000000 +0.292000 12.666667 1.030303 1.150000 1.200000 0.000800 0.000000 +0.296000 12.666667 1.030303 1.150000 1.200000 0.000800 0.000000 +0.300000 12.666667 1.030303 1.150000 1.200000 0.000800 0.000000 +0.304000 13.333333 1.030303 1.150000 1.200000 0.000800 0.000000 +0.308000 13.333333 1.030303 1.150000 1.200000 0.000800 0.000000 +0.312000 13.333333 1.030303 1.150000 1.200000 0.000800 0.000000 +0.316000 14.000000 1.030303 1.150000 1.200000 0.000800 0.000000 +0.320000 14.000000 1.030303 1.150000 1.200000 0.000800 0.000000 +0.324000 14.666667 1.030303 1.150000 1.200000 0.000800 0.000000 +0.328000 14.666667 1.030303 1.150000 1.200000 0.000800 0.000000 +0.332000 15.333333 1.030303 1.150000 1.200000 0.000800 0.000000 +0.336000 15.333333 1.030303 1.150000 1.200000 0.000800 0.000000 +0.340000 16.000000 1.030303 1.150000 1.200000 0.000800 0.000000 +0.344000 16.000000 1.030303 1.150000 1.200000 0.000800 0.000000 +0.348000 16.666667 1.030303 1.150000 1.200000 0.000800 0.000000 +0.352000 16.666667 1.030303 1.150000 1.200000 0.000800 0.000000 +0.356000 17.333333 1.030303 1.150000 1.200000 0.000800 0.000000 +0.360000 17.333333 1.030303 1.150000 1.200000 0.000800 0.000000 +0.364000 18.000000 1.030303 1.150000 1.200000 0.000800 0.000000 +0.368000 18.000000 1.030303 1.150000 1.200000 0.000800 0.000000 +0.372000 18.666667 1.030303 1.150000 1.200000 0.000800 0.000000 +0.376000 18.666667 1.030303 1.150000 1.200000 0.000800 0.000000 +0.380000 19.333333 1.030303 1.150000 1.200000 0.000800 0.000000 +0.384000 19.333333 1.030303 1.150000 1.200000 0.000800 0.000000 +0.388000 20.000000 1.030303 1.150000 1.200000 0.000800 0.000000 +0.392000 20.000000 1.030303 1.150000 1.200000 0.000800 0.000000 +0.396000 20.666667 1.030303 1.150000 1.200000 0.000800 0.000000 +0.400000 21.333333 1.030303 1.150000 1.200000 0.000800 0.000000 +0.000000 0.000000 0.787879 1.300000 1.300000 0.000400 0.000000 +0.004000 0.000000 0.787879 1.300000 1.300000 0.000400 0.000000 +0.008000 0.000000 0.787879 1.300000 1.300000 0.000400 0.000000 +0.012000 0.000000 0.787879 1.300000 1.300000 0.000400 0.000000 +0.016000 0.000000 0.787879 1.300000 1.300000 0.000400 0.000000 +0.020000 0.000000 0.787879 1.300000 1.300000 0.000400 0.000000 +0.024000 0.000000 0.787879 1.300000 1.300000 0.000400 0.000000 +0.028000 0.000000 0.787879 1.300000 1.300000 0.000400 0.000000 +0.032000 0.000000 0.787879 1.300000 1.300000 0.000400 0.000000 +0.036000 0.000000 0.787879 1.300000 1.300000 0.000400 0.000000 +0.040000 0.000000 0.787879 1.300000 1.300000 0.000400 0.000000 +0.044000 0.000000 0.787879 1.300000 1.300000 0.000400 0.000000 +0.048000 0.000000 0.787879 1.300000 1.300000 0.000400 0.000000 +0.052000 0.666667 0.787879 1.300000 1.300000 0.000400 0.000000 +0.056000 2.000000 0.787879 1.300000 1.300000 0.000400 0.000000 +0.060000 2.666667 0.787879 1.300000 1.300000 0.000400 0.000000 +0.064000 2.666667 0.787879 1.300000 1.300000 0.000400 0.000000 +0.068000 3.333333 0.787879 1.300000 1.300000 0.000400 0.000000 +0.072000 3.333333 0.787879 1.300000 1.300000 0.000400 0.000000 +0.076000 4.000000 0.787879 1.300000 1.300000 0.000400 0.000000 +0.080000 4.000000 0.787879 1.300000 1.300000 0.000400 0.000000 +0.084000 4.000000 0.787879 1.300000 1.300000 0.000400 0.000000 +0.088000 4.666667 0.787879 1.300000 1.300000 0.000400 0.000000 +0.092000 4.666667 0.787879 1.300000 1.300000 0.000400 0.000000 +0.096000 4.666667 0.787879 1.300000 1.300000 0.000400 0.000000 +0.100000 5.333333 0.787879 1.300000 1.300000 0.000400 0.000000 +0.104000 5.333333 0.787879 1.300000 1.300000 0.000400 0.000000 +0.108000 5.333333 0.787879 1.300000 1.300000 0.000400 0.000000 +0.112000 6.000000 0.787879 1.300000 1.300000 0.000400 0.000000 +0.116000 6.000000 0.787879 1.300000 1.300000 0.000400 0.000000 +0.120000 6.000000 0.787879 1.300000 1.300000 0.000400 0.000000 +0.124000 6.000000 0.787879 1.300000 1.300000 0.000400 0.000000 +0.128000 6.666667 0.787879 1.300000 1.300000 0.000400 0.000000 +0.132000 6.666667 0.787879 1.300000 1.300000 0.000400 0.000000 +0.136000 6.666667 0.787879 1.300000 1.300000 0.000400 0.000000 +0.140000 7.333333 0.787879 1.300000 1.300000 0.000400 0.000000 +0.144000 7.333333 0.787879 1.300000 1.300000 0.000400 0.000000 +0.148000 7.333333 0.787879 1.300000 1.300000 0.000400 0.000000 +0.152000 8.000000 0.787879 1.300000 1.300000 0.000400 0.000000 +0.156000 8.000000 0.787879 1.300000 1.300000 0.000400 0.000000 +0.160000 8.000000 0.787879 1.300000 1.300000 0.000400 0.000000 +0.164000 8.666667 0.787879 1.300000 1.300000 0.000400 0.000000 +0.168000 8.666667 0.787879 1.300000 1.300000 0.000400 0.000000 +0.172000 8.666667 0.787879 1.300000 1.300000 0.000400 0.000000 +0.176000 9.333333 0.787879 1.300000 1.300000 0.000400 0.000000 +0.180000 9.333333 0.787879 1.300000 1.300000 0.000400 0.000000 +0.184000 9.333333 0.787879 1.300000 1.300000 0.000400 0.000000 +0.188000 10.000000 0.787879 1.300000 1.300000 0.000400 0.000000 +0.192000 10.000000 0.787879 1.300000 1.300000 0.000400 0.000000 +0.196000 10.000000 0.787879 1.300000 1.300000 0.000400 0.000000 +0.200000 10.666667 0.787879 1.300000 1.300000 0.000400 0.000000 +0.204000 10.666667 0.787879 1.300000 1.300000 0.000400 0.000000 +0.208000 11.333333 0.787879 1.300000 1.300000 0.000400 0.000000 +0.212000 11.333333 0.787879 1.300000 1.300000 0.000400 0.000000 +0.216000 11.333333 0.787879 1.300000 1.300000 0.000400 0.000000 +0.220000 12.000000 0.787879 1.300000 1.300000 0.000400 0.000000 +0.224000 12.000000 0.787879 1.300000 1.300000 0.000400 0.000000 +0.228000 12.666667 0.787879 1.300000 1.300000 0.000400 0.000000 +0.232000 9.333333 0.787879 1.300000 1.300000 0.000400 0.000000 +0.236000 9.333333 0.787879 1.300000 1.300000 0.000400 0.000000 +0.240000 10.000000 0.787879 1.300000 1.300000 0.000400 0.000000 +0.244000 10.666667 0.787879 1.300000 1.300000 0.000400 0.000000 +0.248000 10.666667 0.787879 1.300000 1.300000 0.000400 0.000000 +0.252000 11.333333 0.787879 1.300000 1.300000 0.000400 0.000000 +0.256000 11.333333 0.787879 1.300000 1.300000 0.000400 0.000000 +0.260000 12.000000 0.787879 1.300000 1.300000 0.000400 0.000000 +0.264000 12.000000 0.787879 1.300000 1.300000 0.000400 0.000000 +0.268000 12.666667 0.787879 1.300000 1.300000 0.000400 0.000000 +0.272000 13.333333 0.787879 1.300000 1.300000 0.000400 0.000000 +0.276000 14.000000 0.787879 1.300000 1.300000 0.000400 0.000000 +0.280000 16.666667 0.787879 1.300000 1.300000 0.000400 0.000000 +0.284000 17.333333 0.787879 1.300000 1.300000 0.000400 0.000000 +0.288000 17.333333 0.787879 1.300000 1.300000 0.000400 0.000000 +0.292000 18.000000 0.787879 1.300000 1.300000 0.000400 0.000000 +0.296000 18.000000 0.787879 1.300000 1.300000 0.000400 0.000000 +0.300000 18.666667 0.787879 1.300000 1.300000 0.000400 0.000000 +0.304000 18.666667 0.787879 1.300000 1.300000 0.000400 0.000000 +0.308000 19.333333 0.787879 1.300000 1.300000 0.000400 0.000000 +0.312000 19.333333 0.787879 1.300000 1.300000 0.000400 0.000000 +0.316000 20.000000 0.787879 1.300000 1.300000 0.000400 0.000000 +0.320000 20.000000 0.787879 1.300000 1.300000 0.000400 0.000000 +0.324000 20.666667 0.787879 1.300000 1.300000 0.000400 0.000000 +0.328000 21.333333 0.787879 1.300000 1.300000 0.000400 0.000000 +0.332000 21.333333 0.787879 1.300000 1.300000 0.000400 0.000000 +0.336000 22.000000 0.787879 1.300000 1.300000 0.000400 0.000000 +0.340000 22.000000 0.787879 1.300000 1.300000 0.000400 0.000000 +0.344000 22.666667 0.787879 1.300000 1.300000 0.000400 0.000000 +0.348000 23.333333 0.787879 1.300000 1.300000 0.000400 0.000000 +0.352000 23.333333 0.787879 1.300000 1.300000 0.000400 0.000000 +0.356000 24.000000 0.787879 1.300000 1.300000 0.000400 0.000000 +0.360000 24.000000 0.787879 1.300000 1.300000 0.000400 0.000000 +0.364000 24.666667 0.787879 1.300000 1.300000 0.000400 0.000000 +0.368000 25.333333 0.787879 1.300000 1.300000 0.000400 0.000000 +0.372000 25.333333 0.787879 1.300000 1.300000 0.000400 0.000000 +0.376000 26.000000 0.787879 1.300000 1.300000 0.000400 0.000000 +0.380000 26.000000 0.787879 1.300000 1.300000 0.000400 0.000000 +0.384000 26.666667 0.787879 1.300000 1.300000 0.000400 0.000000 +0.388000 27.333333 0.787879 1.300000 1.300000 0.000400 0.000000 +0.392000 27.333333 0.787879 1.300000 1.300000 0.000400 0.000000 +0.396000 28.000000 0.787879 1.300000 1.300000 0.000400 0.000000 +0.400000 28.666667 0.787879 1.300000 1.300000 0.000400 0.000000 diff --git a/ficurves.hoc b/ficurves.hoc new file mode 100644 index 0000000..49072db --- /dev/null +++ b/ficurves.hoc @@ -0,0 +1,99 @@ +// Load cell templates +err_ = load_file("objects/GC.hoc") + +print_voltage_tracings = 1 + +tstop = 3000 +delay = 1000 +duration = tstop - 2*delay + +n_steps = 100 +v_init = -77.71 +i_min = 0 +i_max = 0.4 +i_step = (i_max - i_min)/n_steps + +objref tdata, vdata[n_steps+1], vdatalist, spikes[n_steps+1] +objref current_steps +objref cells[n_steps+1], stim[n_steps+1] +objref vmfile, fifile +strdef vm_filename, ficurve_filename +vm_filename = "ficurve-vm.txt" +ficurve_filename = "ficurve.txt" + + +objref na_scales, kdr_scales, ka_scales, ht_gbars, lt_gbars, size_scales +na_scales = new Vector() +na_scales.append(1, 0.034/0.033, 0.026/0.033) + +kdr_scales = new Vector() +kdr_scales.append(1, 1.15, 1.3) + +ka_scales = new Vector() +ka_scales.append(1, 0.0096/0.008, 0.0104/0.008) + +ht_gbars = new Vector() +ht_gbars.append(0.0004, 0.0008, 0.0004) +//ht_gbars.append(0.000, 0.000, 0.000) + +lt_gbars = new Vector() +//lt_gbars.append(0.0002, 0.0002, 0.0002) +lt_gbars.append(0.000, 0.000, 0.000) + + +for s=0, 2 { + current_steps = new Vector() + for i = 0, n_steps { + cells[i] = new GranuleCell(i, na_scales.x[s], kdr_scales.x[s], ka_scales.x[s], ht_gbars.x[s]) + cells[i].soma stim[i] = new IClamp(0.5) + stim[i].del = delay + stim[i].dur = duration + stim[i].amp = i_min + i*i_step + current_steps.append(i_min + i*i_step) + } + + tdata = new Vector() + tdata.record(&t) + for i = 0, n_steps { + vdata[i] = new Vector() + vdata[i].record(&cells[i].soma.v(0.5)) + spikes[i] = new Vector() + } + + finitialize(v_init) + while (t < tstop) { + fadvance() + } + + + // FUNCTIONS TO PRINT FILES + // Save FI CURVE + fifile = new File() + fifile.aopen(ficurve_filename) + for i = 0, n_steps { + spikes[i].spikebin(vdata[i], 0) + fifile.printf("%f\t%f\t%f\t%f\t%f\t%f\t%f\n", current_steps.x[i], spikes[i].sum()/(duration/1000), na_scales.x[s], kdr_scales.x[s], ka_scales.x[s], ht_gbars.x[s], lt_gbars.x[s]) + } + fifile.close() + + if (print_voltage_tracings == 1) { + // Save voltage traces + vmfile = new File() + sprint(vm_filename, "%svm-%f.txt", "data/fi-curves/", s) + vmfile.aopen(vm_filename) + for j = 0, tdata.size - 1 { + for i = 0, n_steps { + if (i == 0) { + vmfile.printf("%f\t%f\t", tdata.x[j], vdata[i].x[j]) + } + if (i == n_steps) { + vmfile.printf("%f\n", vdata[i].x[j]) + } + if (i > 0 && i < n_steps ) { + vmfile.printf("%f\t", vdata[i].x[j]) + } + } + } + vmfile.close() + } +} diff --git a/fig1.jpg b/fig1.jpg deleted file mode 100644 index 1d6fe49..0000000 Binary files a/fig1.jpg and /dev/null differ diff --git a/fig1.py b/fig1.py deleted file mode 100644 index fd1c9d8..0000000 --- a/fig1.py +++ /dev/null @@ -1,18 +0,0 @@ -### Analysis of DG network data ### -# This Python code creates a scatter plot of output vs input sim scores. -# Enter the idname - -# ModelDB file along with publication: -# Yim MY, Hanuschkin A, Wolfart J (2015) Hippocampus 25:297-308. -# http://onlinelibrary.wiley.com/doi/10.1002/hipo.22373/abstract - -# modified and augmented by -# Man Yi Yim / 2015 -# Alexander Hanuschkin / 2011 - -idname = "-pp10-gaba1-kir1-st0" -execfile('plot_DG_all.py') -execfile('GCinput.py') -execfile('inout_pattern.py') -execfile('sim_score.py') -print 'Set the idname as ', idname, ' in FitSimScore_ForallFigs.m and run the file in Matlab for data fitting' \ No newline at end of file diff --git a/fig2.jpg b/fig2.jpg deleted file mode 100644 index 4ad2246..0000000 Binary files a/fig2.jpg and /dev/null differ diff --git a/fig2.py b/fig2.py deleted file mode 100644 index 2e4ac3e..0000000 --- a/fig2.py +++ /dev/null @@ -1,16 +0,0 @@ -### Analysis of DG network data ### -# This Python code creates a scatter plot of output vs input sim scores. -# Enter the idname - -# ModelDB file along with publication: -# Yim MY, Hanuschkin A, Wolfart J (2015) Hippocampus 25:297-308. -# http://onlinelibrary.wiley.com/doi/10.1002/hipo.22373/abstract - -# modified and augmented by -# Man Yi Yim / 2015 -# Alexander Hanuschkin / 2011 - -idname = "-pp10-gaba1-kir1-st30" -execfile('plot_DG_all.py') -execfile('GCinput.py') -execfile('inout_pattern.py') \ No newline at end of file diff --git a/fig3.jpg b/fig3.jpg deleted file mode 100644 index 6f76a63..0000000 Binary files a/fig3.jpg and /dev/null differ diff --git a/fig3.py b/fig3.py deleted file mode 100644 index bd10ead..0000000 --- a/fig3.py +++ /dev/null @@ -1,16 +0,0 @@ -### Analysis of DG network data ### -# This Python code creates a scatter plot of output vs input sim scores. -# Enter the idname - -# ModelDB file along with publication: -# Yim MY, Hanuschkin A, Wolfart J (2015) Hippocampus 25:297-308. -# http://onlinelibrary.wiley.com/doi/10.1002/hipo.22373/abstract - -# modified and augmented by -# Man Yi Yim / 2015 -# Alexander Hanuschkin / 2011 - -idname = "-pp10-gaba4-kir4-st30" -execfile('plot_DG_all.py') -execfile('inout_pattern.py') -print "Change step = 1 in inout_pattern.py to produce the same result as in Yim et al. (2015)" \ No newline at end of file diff --git a/fig4.py b/fig4.py deleted file mode 100644 index 65d9aa6..0000000 --- a/fig4.py +++ /dev/null @@ -1,18 +0,0 @@ -### Analysis of DG network data ### -# This Python code creates a scatter plot of output vs input sim scores. -# Enter the idname - -# ModelDB file along with publication: -# Yim MY, Hanuschkin A, Wolfart J (2015) Hippocampus 25:297-308. -# http://onlinelibrary.wiley.com/doi/10.1002/hipo.22373/abstract - -# modified and augmented by -# Man Yi Yim / 2015 -# Alexander Hanuschkin / 2011 - -idname = "-pp16-gaba1-kir1-st10" -execfile('plot_DG_all.py') -execfile('GCinput.py') -execfile('inout_pattern.py') -execfile('sim_score.py') -print 'Set the idname as ', idname, ' in FitSimScore_ForallFigs.m and run the file in Matlab for data fitting' \ No newline at end of file diff --git a/fig5.py b/fig5.py deleted file mode 100644 index 8add789..0000000 --- a/fig5.py +++ /dev/null @@ -1,18 +0,0 @@ -### Analysis of DG network data ### -# This Python code creates a scatter plot of output vs input sim scores. -# Enter the idname - -# ModelDB file along with publication: -# Yim MY, Hanuschkin A, Wolfart J (2015) Hippocampus 25:297-308. -# http://onlinelibrary.wiley.com/doi/10.1002/hipo.22373/abstract - -# modified and augmented by -# Man Yi Yim / 2015 -# Alexander Hanuschkin / 2011 - -idname = "-pp16-gaba4-kir4-st10" -execfile('plot_DG_all.py') -execfile('GCinput.py') -execfile('inout_pattern.py') -execfile('sim_score.py') -print 'Set the idname as ', idname, ' in FitSimScore_ForallFigs.m and run the file in Matlab for data fitting' \ No newline at end of file diff --git a/find-rheobase.py b/find-rheobase.py new file mode 100644 index 0000000..7004f0e --- /dev/null +++ b/find-rheobase.py @@ -0,0 +1,114 @@ +import numpy as np +import pandas as pd +import matplotlib.pyplot as plt +from netpyne import sim, specs +from scipy.optimize import golden +from clamps import VClamp, IClamp + + +netparams = specs.NetParams() +gc = netparams.importCellParams( + label="GC", + conds={"cellType": "GranuleCell", "cellModel": "GranuleCell"}, + fileName="objects/GC.hoc", + cellName="GranuleCell", + cellArgs=[1], + importSynMechs=False +) + +class ElectrophysiologicalPhenotype(object): + """ This object computes a wide variety of electrophysiological properties for + a given neuron object. + """ + def __init__(self, cell): + self.cell_dict = {"secs": cell["secs"]} + self.fi_data = {} + self.fi_curve = None + self.rheo_current_brack = None + + def step_current(self, current, delay=250, duration=500): + """ Injects a level of current and returns the number of spikes emitted + + Arguments: + current: `float`. Amount of current injected [nA] + delay: `float`. Time after recording starts where current is injected [ms] + duration: `float`. Total duration of the current injection [ms] + + Returns: + `dict`. Results of the step current injection simulation + + """ + iclamp = IClamp(self.cell_dict, delay=delay, duration=duration, T=duration + delay*2) + res = iclamp(current) + return res + + def compute_fi_curve(self, ilow=0, ihigh=0.5, n_steps=100, delay=250, duration=500): + """ Computes the fi curve, and stores the raw data + + Arguments: + ilow: `float`. Starting current injection + ihigh: `float`. Top current injection + n_steps: `int`. Number of current injection steps + delay: `float`. Time after recording starts where current is injected [ms] + duration: `float`. Total duration of the current injection [ms] + + Returns: + `pandas.DataFrame`. + """ + self.fi_curve = pd.DataFrame(np.zeros((n_steps, 2)), columns=["I", "F"]) + current_steps = np.linspace(ilow, ihigh, n_steps) + self.fi_curve["I"] = current_steps + + for j, current in enumerate(current_steps): + self.fi_data[j] = self.step_current(current, delay=delay, duration=duration) + self.fi_data[j]["current"] = current + self.fi_curve.iloc[j, 0] = current + self.fi_curve.iloc[j, 1] = self.fi_data[j]["rate"] + + self._get_rheobase_bracket() + + return self.fi_curve + + def _get_rheobase_bracket(self): + """ Finds the initial bracket for the rheobase calculation based on the fI curve """ + ilow = np.max(self.fi_curve.loc[self.fi_curve.F == 0, "I"].values) + ihigh = np.min(self.fi_curve.loc[self.fi_curve.F > 0, "I"].values) + self.rheo_current_brack = (ilow, ihigh) + + def find_rheobase(self, current_brack=(0, 1), delay=250, duration=500, tol=1e-3): + """ Finds the rheobase of the clamped neuron using Golden Section Search, + minimizing the loss function ((n_spikes-1) + I)^2, where I is the + injected current level + + Arguments: + current_brack: `tuple` or `None`. (low, mid, high) for golden section search + delay: `float`. Time after recording starts where current is injected [ms] + duration: `float`. Total duration of the current injection [ms] + tol: `float`. The tolerance level + + Returns: + `float`. Rheobase + """ + # Use current bracket defined from F-I curve preferentially + if self.rheo_current_brack is None: + self.rheo_current_brack = current_brack + + # Define the loss function + rheo_loss = lambda i: (( len(self.step_current(i, delay, duration)["spkt"]) - 1) + i)**2 + + # Perform golden section search + self.rheobase = golden(rheo_loss, brack=self.rheo_current_brack, tol=tol) + + return self.rheobase + + +# gcpheno = ElectrophysiologicalPhenotype(gc) +# fi = gcpheno.compute_fi_curve(n_steps=20, delay=500, duration=1000) +# gcpheno.find_rheobase() + +# i = 2 +# fig, ax = plt.subplots(nrows=2) +# ax[0].plot(gcpheno.fi_data[i]["t"], gcpheno.fi_data[i]["V"]) +# ax[1].plot(gcpheno.fi_data[i]["t"], gcpheno.fi_data[i]["i"]) +# plt.show() +# plt.close() \ No newline at end of file diff --git a/inout_pattern.py b/inout_pattern.py deleted file mode 100644 index 86dc4d5..0000000 --- a/inout_pattern.py +++ /dev/null @@ -1,81 +0,0 @@ -### Analysis of DG network data ### -# This Python code creates a scatter plot of output vs input sim scores. -# Enter the idname - -# ModelDB file along with publication: -# Yim MY, Hanuschkin A, Wolfart J (2015) Hippocampus 25:297-308. -# http://onlinelibrary.wiley.com/doi/10.1002/hipo.22373/abstract - -# modified and augmented by -# Man Yi Yim / 2015 -# Alexander Hanuschkin / 2011 - -import pylab -import numpy -import matplotlib as mpl -font_size = 12 # 20 -mpl.rcParams['axes.titlesize'] = font_size+10 -mpl.rcParams['xtick.labelsize'] = font_size+6 -mpl.rcParams['ytick.labelsize'] = font_size+6 -mpl.rcParams['axes.labelsize'] = font_size+8 -mpl.rcParams['legend.fontsize'] = font_size -mpl.rcParams['font.size'] = font_size+10 - -sid = 0 -sidm = 12 -nid = 6 -step = 2 -#delta = numpy.array([5,10,20,30])/2. # half the width of the kernel base -delta = numpy.array([20])/2. -dur = 200. # duration of the simulation -bin = 1. -ninput = 100 # number of input sources -ncell = 500 # number of neurons (GCs) -#idname = "-pp10-gaba1-kir1-st0" - -pylab.figure(figsize=(15,16)) -for k in range(4): #number to show - pylab.subplot(4,1,k+1) - stimin = numpy.loadtxt('StimIn-'+str(step*k)+'-'+str(nid)+idname+'.txt') - for j in range(stimin.size/2): - if stimin.size == 2: - pylab.plot(stimin[0]*numpy.ones(2),numpy.array([stimin[1]-0.4,stimin[1]+0.4]),'k',linewidth=2) - else: - pylab.plot(stimin[j,0]*numpy.ones(2),numpy.array([stimin[j,1]-0.4,stimin[j,1]+0.4]),'k',linewidth=2) - pylab.xticks([]) - pylab.yticks([]) - pylab.xlim(0.,dur) - pylab.ylabel(str(step*k+1)) - pylab.axis([0,dur,-0.5,19-0.5]) - if k == 0: - pylab.title('Pattern'+', '+idname) - if k == 12: - pylab.xticks(range(0,201,100)) - pylab.xlabel('Time (ms)') -pylab.xticks(range(0,201,50)) -pylab.savefig('PP_'+idname+'.eps') - -############### -pylab.figure(figsize=(15,16)) -for k in range(4): - pylab.subplot(4,1,k+1) - dgsp = numpy.loadtxt('DGsp-'+str(step*k)+'-'+str(nid)+idname+'.txt') - for j in range(dgsp.size/2): - if dgsp.size == 2: - pylab.plot(dgsp[0]*numpy.ones(2),numpy.array([dgsp[1]-0.4,dgsp[1]+0.4]),color,linewidth=2) - else: - pylab.plot(dgsp[j,0]*numpy.ones(2),numpy.array([dgsp[j,1]-0.4,dgsp[j,1]+0.4]),'b',linewidth=2) # color for GC - pylab.xticks([]) - pylab.yticks([]) - pylab.xlim(0.,dur) - pylab.ylabel(str(step*k+1)) - pylab.axis([0,dur,-0.5,500-0.5]) - if k == 0: - pylab.title('GC'+', '+idname) - if k == 12: - pylab.xticks(range(0,201,100)) - pylab.xlabel('Time (ms)') -pylab.xticks(range(0,201,50)) -pylab.savefig('GCoutput_'+idname+'.eps') -############### -pylab.show() diff --git a/main.hoc b/main.hoc index 254c90c..18945f5 100644 --- a/main.hoc +++ b/main.hoc @@ -1,1099 +1,33 @@ -// This is a fully wired network that functions with 500 Granule Cells, 15 Mossy cells, 6 Basket cells and 6 HIPP Cells - -// modified version of -// Dentate gyrus network model -// Santhakumar V, Aradi I, Soltesz I (2005) J Neurophysiol 93:437-53 -// https://senselab.med.yale.edu/ModelDB/showModel.cshtml?model=51781&file=\dentategyrusnet2005\DG500_M7.hoc - -// ModelDB file along with publication: -// Yim MY, Hanuschkin A, Wolfart J (2015) Hippocampus 25:297-308. -// http://onlinelibrary.wiley.com/doi/10.1002/hipo.22373/abstract - -// modified by -// Man Yi Yim / 2015 -// Alexander Hanuschkin / 2011 - -// Which figure (1-5) do you want to simulate? -// Please use the original parameters to produce the figures. -// If you want to run the simulation with your own parameters, set fig = 0 -fig = 1 - -// SIMULATION: parameters -dt = 0.1 // simulation time [ms] -tstop = 200 // total simulation time [ms] -trial = 1 // trialnumber = seed of simulation -trial_noise = 1 // seed for noise to PP - -// NOTE: // all NetStims share the same random number generator -// used solution given in http://www.neuron.yale.edu/neuron/node/61 -err_ = load_file("ranstream.hoc") - -objref cvode -cvode = new CVode() -using_cvode_ = 1 -debug_ = 2 // print output - // everything: debug==1 - // some: debug==0 - -// FLAGS: Input -p_sprouted_ = 0 // 10% sprouting (i.e. each GC connects to 50 other GCs randomly...) // original p = 0.02 -PP_nstimulus_ = 1 // number of stimulated GC (100 in original file), set to 1 to have one input per GC (e.g. uncorrelated Poisson or correlated input (Myers and Scharfman, 2009) -PP_input2MC_ = 0 // stimulate MC with PP -PP_rate_ = 10. // rate of PP Poisson input -PP_rate_noise_ = 0. // rate of PP Poisson noise -PP_box_nr_ = 6 // number of active PPs -PP_box_stop_ = 35. // time of box [ms] -PP_box_start_ = 5. // shift of box [ms] -PP_box_nspk_ = 3 // number of spike per active PP - - -// FLAGS: Output -print_Vtrace = 1 // print voltage trace to file -print_Raster = 1 // print spike raster to file -print_template = 1 // write out spike sequence into template file -print_GC_sprouting_input_ = 0 // calculates and print out GC sprouting input number distribution -print_stim_ = 1 // print stimulus to GCs for debug purpose.. -print_stim_noise = 1 // print noise to GCs for debug purpose.. (same output file as Poisson input to GCs!) - - -// FLAGS: Scaling -scale_gpas_dg_ = 100 // scaling [%] of gpas of th GC model (100% in the original file) -scale_sk_dg_ = 100 // scaling of Ca-dependent K (SK) of th GC model (100% in the original file) -scale_kir_ = 100 // scaling of KIR conductance in GC model -scale_gabaa_ = 100 // scaling of GABAA in GC model -scale_PP_strength = 10 // scaling of synaptic weight of PP -scale_PP2BC_strength = 100 // scaling of synaptic weight PP->BC (100% for original value which was 50% of PP->GC strength..) -scale_PP2HC_strength = 0 // scaling of synaptic weight PP->HIPP (set to 0% for Santhakumar because there no PP->HIPP connections...) -scale_HC2GC_strength = 100 // scaling of synaptic weight HC->GC (beta_HIPP in Myers and Scharfman, 2009) -scale_MC2GC_strength = 100 // scaling of synaptic weight MC->GC (beta_HIPP in Myers and Scharfman, 2009) -scale_GC2BC_strength = 300 // scaling of synaptic weight GC->BC -scale_BC2GC_strength = 300 // scaling of synaptic weight BC->GC - - -init_from_file_ = 0 // read in initial potential of all cells/compartments from file -init_to_file_ = 0 // save potential of all cells/compartments to file @ end of run (to be used for init the network later again) - -strdef suffix, idname // define file names for simulation results output -suffix = "txt" -idname = "-pp10-gaba1-kir1-st0" -strdef init_state_fname_ -init_state_fname_ = "NetworkState_init_10Hz.dat" - -// Reset parameters to generate the figures in Yim et al (2015) -if (fig == 1) { - scale_PP_strength = 10 - scale_kir_ = 100 - scale_gabaa_ = 100 - p_sprouted_ = 0 - idname = "-pp10-gaba1-kir1-st0" -} - -if (fig == 2) { - scale_PP_strength = 10 - scale_kir_ = 100 - scale_gabaa_ = 100 - p_sprouted_ = 30 - idname = "-pp10-gaba1-kir1-st30" -} - -if (fig == 3) { - scale_PP_strength = 10 - scale_kir_ = 400 - scale_gabaa_ = 400 - p_sprouted_ = 30 - idname = "-pp10-gaba4-kir4-st30" -} - -if (fig == 4) { - scale_PP_strength = 16 - scale_kir_ = 100 - scale_gabaa_ = 100 - p_sprouted_ = 10 - idname = "-pp16-gaba1-kir1-st10" -} - -if (fig == 5) { - scale_PP_strength = 16 - scale_kir_ = 400 - scale_gabaa_ = 400 - p_sprouted_ = 10 - idname = "-pp16-gaba4-kir4-st10" -} -// End of parameter reset - -// NETWORK: parameters -ngcell = 500 // number of GCs -nbcell = 6 // number of BCs (6 for original Samathakumar 2005) -nmcell = 15 // number of MCs -nhcell = 6 // number of HCs (6 for original Samathakumar 2005) -npp = 100 // 100 enorhinal layer II Neurons (Myers and Scharfman, 2009) - -ntotal = ngcell + nbcell + nmcell + nhcell + npp + npp // two times npp because of PP input and PP noise... - -if (PP_input2MC_ != 0) { - print "PP stimulation not yet implemented if stimulation of also MCs...." - quit() -} - -objref PPSt[npp], PPSt_noise[npp] -create aacell // artificial cell if windows (boxed) activation -objref pp_start -objref netcon_d[npp] - -// for debug purpose only... -objref netcon, nil -objref vec_debug_ -objref PP_init_rnd // random number generator to select an PP stimulus site.. -objref nslist, rslist, rs_noiselist -objref ns, rs - -//Generate spike matrix and save to file -objref Spike[ntotal-1], Spike_times[ntotal-1] -for i=0, (ngcell+nbcell +nmcell +nhcell -1) { - Spike[i] = new Vector() // vector of spikes for each cell - Spike_times[i] = new Vector() // vector of spikes for each cell -} -strdef Spkstr -objref dfile -dfile = new File() - -//********************** GRANULE CELL **************************************** -err_ = load_file("GC.hoc") // moved the granule cell definition to an external file => can be substituted by own reconstructed cell/ different model etc.. - -// ********* BASKET CELL ************************************** -err_ = load_file("BC.hoc") - -//**************** MOSSY CELL *********************************************************** -err_ = load_file("MC.hoc") - -//*************** HIPP CELL **************************** -err_ = load_file("HIPP.hoc") - -//********** Perforant Path Stimulus *********************************************** -err_ = load_file("PP.hoc") - -//############################################################################################################### -//############## CONNECTING THE CELLS ############################# - -// NETWORK SPECIFICATION INTERFACE -for i=0, ngcell-1 {Gcell[i] = new GranuleCell(i)} -for i=0, nbcell-1 {Bcell[i] = new BasketCell(i)} -for i=0, nmcell-1 {Mcell[i] = new MossyCell(i)} -for i=0, nhcell-1 {Hcell[i] = new HIPPCell(i)} -for i=0, npp-1 {PPSt[i] = new PPstim(i)} -for i=0, npp-1 {PPSt_noise[i] = new PPstim(i)} - -objref nclist, netcon, cells, net_c, net_d, net_gr, net_bc, net_mc, net_hc, vbc2gc, vmc2gc, vhc2gc -{ cells = new List() - nclist = new List() -} - -// Include all cells in cells -func cell_append() { - cells.append($o1) - return cells.count -1 +// Hyperparameters +N_RUNS = 10 // Number of runs per group + +// Load cell templates and network +err_ = load_file("objects/GC.hoc") +err_ = load_file("objects/BC.hoc") +err_ = load_file("objects/MC.hoc") +err_ = load_file("objects/HIPP.hoc") +err_ = load_file("objects/PP.hoc") +err_ = load_file("objects/DentateGyrus.hoc") + +// Instantiate cell objects and labels +objref hc, lr, nr +strdef hclab, lrlab, nrlab + +// Do N_RUNS Runs for each network setting +for random_state = 1, N_RUNS { + sprint(hclab, "%s-%d", "HC", random_state) + sprint(lrlab, "%s-%d", "LR", random_state) + sprint(nrlab, "%s-%d", "NR", random_state) + hc = new DentateGyrus(hclab, random_state, "HC") + lr = new DentateGyrus(lrlab, random_state, "LR") + nr = new DentateGyrus(nrlab, random_state, "NR") + + print "RUNNING DENTATE GYRUS ", hclab + hc.run() + + print "RUNNING DENTATE GYRUS ", lrlab + lr.run() + + print "RUNNING DENTATE GYRUS ", nrlab + nr.run() } - -func nc_append() { // neuron connect $1 with $2.pre_list.object($3), weight $4, delay $5, threshold $6 - // connects: - // cells.object($1) with - // $o1 = cells.object($2).pre_list.object($3) and - // returns: - // netcon = $o2 - - if ($3 >= 0 ) { - cells.object($1).connect_pre(cells.object($2).pre_list.object($3),netcon) // connect_pre is function in the respective cell definition - netcon.weight = $4 netcon.delay = $5 netcon.threshold = $6 - } - nclist.append(netcon) - return nclist.count-1 -} - -func nc_append_rec() { // neuron connect $1 with $2.pre_list.object($3), weight $4, delay $5, threshold $6 - // connects: - // cells.object($1) with - // $o1 = cells.object($2).pre_list.object($3) and - // returns: - // netcon = $o2 - // record events to $o7 - - if ($3 >= 0 ) { - cells.object($1).connect_pre(cells.object($2).pre_list.object($3),netcon) // connect_pre is function in the respective cell definition - netcon.weight = $4 netcon.delay = $5 netcon.threshold = $6 - netcon.record($o7) // the difference is here! - } - nclist.append(netcon) - return nclist.count-1 -} - - -// To check for preexisting connections between randomly selected cells -// to avoid multiple contacts between same 2 cells -func is_connected() {local i, c - c=0 - for i=0, nclist.count-1 { - net_c= nclist.object(i) - if (($o1 == net_c.postcell()) && ($o2 == net_c.precell())) {c=1} - } - return c -} - - -objref vbc2gc, vmc2gc, vhc2gc, vgc2bc, vbc2bc, vmc2bc, vhc2bc, vgc2mc, vbc2mc, vmc2mc, vhc2mc, vgc2hc, vmc2hc,vgc2gc -//To delete certain randomly selected cells from net - in this case 8 of 15 MCs -objref killMC - { - vgc2bc = new Vector(nbcell, 0) // defines vectors for each "pre2post" pair, vector length is same as number of post cells fills with 0 - vbc2bc = new Vector(nbcell, 0) // increment x[i] for each connection made! - vmc2bc = new Vector(nbcell, 0) - vhc2bc = new Vector(nbcell, 0) - - killMC = new Vector(8, 0) //To delete certain randomly selected cells from net - in this case 8 of 15 MCs - vgc2mc = new Vector(nmcell, 0) - vbc2mc = new Vector(nmcell, 0) - vmc2mc = new Vector(nmcell, 0) - vhc2mc = new Vector(nmcell, 0) - - vgc2hc = new Vector(nhcell, 0) - vmc2hc = new Vector(nhcell, 0) - - vbc2gc = new Vector(ngcell, 0) - vmc2gc = new Vector(ngcell, 0) - vhc2gc = new Vector(ngcell, 0) - vgc2gc = new Vector(ngcell, 0) - } - - - - -//initiating random number generator for each pre2post pair -//also randomly select MC to kill "deadMC" - -objref rdsynb, rdsyna, rdgc2hc, rdgc2bc, rdgc2mc, rdbc2gc, rdbc2bc, rdbc2mc, deadMC -objref rdmc2gc1, rdmc2gc2, rdmc2bc, rdmc2mc, rdmc2mc1, rdmc2hc, rdhc2gc, rdhc2bc, rdhc2mc, rdgc2gc -objref pp2gc // connectivity vector PP -> GC -objref rnd_pp2gc // new random number generator for drawing connections -objref pp2bc // connectivity vector PP -> BC -objref rnd_pp2bc // new random number generator for drawing connections -objref pp2hc // connectivity vector PP -> HC -objref rnd_pp2hc // new random number generator for drawing connections - - -err_ = ropen("/proc/uptime") // get a seed that is changing based on the processing time - { - //rseed = fscan() // so simulation will not start with the same seed (original) // Fscan() reads the next number from the file opened with the ropen() command - rseed = trial // set the seed to trialnumber, to get reproduceable connectivity and responses - ropen() - } - - -//************************************ PP *********************************************** -rnd_pp2gc = new Random(rseed) -proc new_rnd_pp2gc() {rnd_pp2gc.discunif(0,npp-1)} -new_rnd_pp2gc() - -rnd_pp2bc = new Random(rseed) -proc new_rnd_pp2bc() {rnd_pp2bc.discunif(0,npp-1)} -new_rnd_pp2bc() - -rnd_pp2hc = new Random(rseed) -proc new_rnd_pp2hc() {rnd_pp2hc.discunif(0,npp-1)} -new_rnd_pp2hc() - - -//************************************ GC *********************************************** -rdgc2bc = new Random(rseed) // use for syn.connections -proc new_rdgc2bc() {rdgc2bc.discunif(-1,1)} // range is based on spread of connections on either side of RING- precell loc =0 -new_rdgc2bc() -rdgc2mc = new Random(rseed) // use for syn.connections -proc new_rdgc2mc() {rdgc2mc.discunif(0,2)} -new_rdgc2mc() -rdgc2hc = new Random(rseed) // use for syn.connections -proc new_rdgc2hc() {rdgc2hc.discunif(-2,2)} -new_rdgc2hc() -rdgc2gc = new Random(rseed) // use for syn.connections - -//proc new_rdgc2gc() {rdgc2gc.discunif(-50, 50)} // search around a GC ... However this has to be changed for large sprounting values -// n_spout_ = int (p_sprouted_ ) -1 // NOTE: 100% Sprouting = 100 new synapses!!! (compare p. 443 in Santhakumar paper) -proc new_rdgc2gc() {rdgc2gc.discunif(-50, 50)} - - -new_rdgc2gc() - -//************************************ BC *********************************************** -rdbc2gc = new Random(rseed) // use for syn.connections -proc new_rdbc2gc() {rdbc2gc.discunif(-70, 70)} // range is based on spread of connections on either side of RING- precell loc =0 -new_rdbc2gc() -rdbc2bc = new Random(rseed) // use for syn.connections -proc new_rdbc2bc() {rdbc2bc.discunif(-1, 1)} -new_rdbc2bc() -rdbc2mc = new Random(rseed) // use for syn.connections -proc new_rdbc2mc() {rdbc2mc.discunif(-3, 3)} -new_rdbc2mc() - -//************************************* MC ******************************************** -deadMC = new Random(rseed) // randomly select MC to kill "deadMC" -proc new_deadMC() { - deadMC.discunif(ngcell+nbcell, ngcell+nbcell+nmcell-1) -} -new_deadMC() - -rdmc2gc1 = new Random(rseed) // use for syn.connections -proc new_rdmc2gc1() {rdmc2gc1.discunif(25, 175)} // range is based on spread of connections on either side of RING- precell loc =0 -new_rdmc2gc1() -rdmc2gc2 = new Random(rseed) // use for syn.connections -proc new_rdmc2gc2() {rdmc2gc2.discunif(-175, -25)} -new_rdmc2gc2() -rdmc2bc = new Random(rseed) // use for syn.connections -proc new_rdmc2bc() {rdmc2bc.discunif(-3,3)} -new_rdmc2bc() -rdmc2mc = new Random(rseed) // use for syn.connections -proc new_rdmc2mc() {rdmc2mc.discunif(ngcell+nbcell, ngcell+nbcell+nmcell-1)} -new_rdmc2mc() -rdmc2mc1 = new Random(rseed) // use for syn.connections -proc new_rdmc2mc1() {rdmc2mc1.discunif(-3, 3)} -new_rdmc2mc1() -rdmc2hc = new Random(rseed) // use for syn.connections -proc new_rdmc2hc() {rdmc2hc.discunif(-2, 2)} -new_rdmc2hc() -//************************************* HC ******************************************** - -rdhc2gc = new Random(rseed) // use for syn.connections -proc new_rdhc2gc() {rdhc2gc.discunif(-130, 130)} -new_rdhc2gc() -rdhc2bc = new Random(rseed) // use for syn.connections -proc new_rdhc2bc() {rdhc2bc.discunif(-2, 2)} -new_rdhc2bc() -rdhc2mc = new Random(rseed) // use for syn.connections -proc new_rdhc2mc() {rdhc2mc.discunif(-2, 2)} -new_rdhc2mc() - -//**************** randomizer for the dendritic location of synapse ************************************** - -rdsyna = new Random(rseed) // initialize random distr. -proc new_rdsyna() {rdsyna.discunif(0, 1)} // randomize among 2 dendrites -new_rdsyna() - -rdsynb = new Random(rseed) // initialize random distr. -proc new_rdsynb() {rdsynb.discunif(0, 3)} // randomize among 4 dendrites -new_rdsynb() - - -// NETWORK INITIATION - for i = 0, ngcell-1 {cell_append(Gcell[i])} // cells 0-499 GCs - for i = 0, nbcell-1 {cell_append(Bcell[i])} // cells 500-505 BC - for i = 0, nmcell-1 {cell_append(Mcell[i])} // cells 506-520 MC - for i = 0, nhcell-1 {cell_append(Hcell[i])} // cells 521-526 HC - for i = 0, npp-1 {cell_append(PPSt[i])} // 527 - xxx PP artificial cell - for i = 0, npp-1 {cell_append(PPSt_noise[i])} // 527 - xxx PP artificial cell for noise - - -// STIM -// record input stimulus for debug purpose -objref vec_stim[npp] -for i = 0, npp-1 { vec_stim[i] = new Vector() } -// record input noise for debug purpose -objref vec_stim_noise[npp] -for i = 0, npp-1 { vec_stim_noise[i] = new Vector() } -objref pprec // vector if recorded from PP -// record Poisson in for debug purpose - - -objref distri_gc_input_ -distri_gc_input_ = new Vector() - - - -proc initNet_GC() { local i,j - -//**************Preforant Path synaptic connections ****************************** -if (debug_ == 2 ) { - print "Preforant Path synaptic connections" - print "Correlated input from 100 EC Layer 2 cells" -} - - pprec = new Vector(npp, 0) - - // ################ PP -> GC - for i=0, ngcell-1 { - pp2gc = new Vector(npp, 0) // init connectivity vector to prevent multiple stimulation connections from the same PP cell - for nr_conn = 0, int(npp/5.)-1 { // select randomly 20% of the PP for input - j = rnd_pp2gc.repick() // id of PP cell - while (pp2gc.x[j] == 1) { // random drawing until unconnected PP was found - j = rnd_pp2gc.repick() // id of PP cell - } - pp2gc.x[j] = 1 - if ((print_stim_==1) && (i NOTE both synapses are equal, ie. weight delay (except position) not important - } - } - if (debug_ == 3) { - print "\nGC: ",i - for ii=0,pp2gc.size()-1{ printf ("%d, ",pp2gc.x[ii])} - } - } - - - // ################ PP -> BC (input to two BC cells as in Santhakumar 2005) - - if (scale_PP2BC_strength != 0) { // if PP2BC connections - for i=ngcell, ngcell+5 { - pp2bc = new Vector(npp, 0) // init connectivity vector to prevent multiple stimulation connections from the same PP cell - for nr_conn = 0, int(npp/5.)-1 { // select randomly 20% of the PP for input - j = rnd_pp2bc.repick() // id of PP cell - // print j - while (pp2bc.x[j] == 1) { // random drawing until unconnected PP was found - j = rnd_pp2bc.repick() // id of PP cell - } - pp2bc.x[j] = 1 - - nc_append(j+ngcell+nbcell+nmcell+nhcell, i, 0, 1e-2*scale_PP_strength/100.*scale_PP2BC_strength/100., 3, 10) // connect PP to GC[j],syn[0],wt,del,threshold NOTE both synapses are equal -// nc_append(j+ngcell+nbcell+nmcell+nhcell+npp, i, 0, 1e-2*scale_PP_strength/100.*scale_PP2BC_strength/100., 3, 10) // add noise to the GCs - } - if (debug_ == 3) { - print "\nBC: ",i - for ii=0,pp2bc.size()-1{ printf ("%d, ",pp2gc.x[ii])} - } - } - } - - // ################ PP -> HIPP (input to 20% HIPP cells as in Myers and Scharfman 2009) - - if (scale_PP2HC_strength != 0) { // if PP2HC connections - if (debug_ == 2) {print "PP -> HIPP (input to two HIPP cells as in Myers and Scharfman 2009)"} - for i=0, nhcell-1 { - pp2hc = new Vector(npp, 0) // init connectivity vector to prevent multiple stimulation connections from the same PP cell - for nr_conn = 0, int(npp/5.)-1 { // select randomly 20% of the PP for input - j = rnd_pp2hc.repick() // id of PP cell - // print j - while (pp2hc.x[j] == 1) { // random drawing until unconnected PP was found - j = rnd_pp2hc.repick() // id of PP cell - } - pp2hc.x[j] = 1 - - nc_append(j+ngcell+nbcell+nmcell+nhcell, i+ngcell+nbcell+nmcell, 0, 0.5e-3*scale_PP_strength/100.*scale_PP2HC_strength/100., 3, 10) // used here synaptic strength of GC->HIPP / connections PP-> HIPP not in Santhakumar 2005 -// nc_append(j+ngcell+nbcell+nmcell+nhcell+npp,i+ngcell+nbcell+nmcell, 0, 0.5e-3*scale_PP_strength/100.*scale_PP2HC_strength/100., 3, 10) // add noise to the PP -> HIPP - - } - if (debug_ == 3) { - print "\nHIPP: ",i - for ii=0,pp2bc.size()-1{ printf ("%d, ",pp2hc.x[ii])} - } - } - if (debug_ == 3) { print "PP -> HIPP -DONE-" } - } - - -//****************************************************************************************** - - - -//**************Granule Cell post synaptic connections ****************************** - - - - -if (debug_ == 2 ) { print "Granule Cell post synaptic connections"} -if (debug_ == 0) {print "n_sprout_ = ",n_sprout_} -if (p_sprouted_>0 && debug_ == 2) {print "NOTE: we have sprouting connections..."} - -for i=0, ngcell-1 { - - for j=0, 0 { - // Based on the lamellar distribution of the GCs to BCs - 500 GCs were divided into 6 groups proximal to each of the 6(!) BCs - if (i < ngcell/6) { a=0} - if ((i > ngcell/6-1) && (i < ngcell*2/6)) { a=1} - if ((i > ngcell*2/6-1) && (i < ngcell*3/6)) { a=2} - if ((i > ngcell*3/6-1) && (i < ngcell*4/6)) { a=3} - if ((i > ngcell*4/6-1) && (i < ngcell*5/6)) { a=4} - if ((i > ngcell*5/6-1) && (i < ngcell)) { a=5} - - Gauz3 = rdgc2bc.repick() // randomly pick location of post synaptic Bcell from distribution [-1:1] - if (a+Gauz3 > nbcell-1) {npost = a+Gauz3-nbcell } // determine appropriate post syn BC - if (a+Gauz3 < 0) {npost = a+Gauz3+nbcell} - if ((a+Gauz3 > -1) && (a+Gauz3 < nbcell)) {npost = a+Gauz3} - dbr = rdsynb.repick() // randomly pick the dendrite to connect to from [0:3] (i.e. // randomize among 4 dendrites) - if (debug_ == 1 ) { print "GC \t",i,"\t to BC\t\t",npost, a} - if (vgc2bc.x[npost] < 90) { // check to make sure that post syn BC does not get more than 90 GC synapses - nc_append(i, ngcell+npost, dbr+2, 4.7e-3*scale_GC2BC_strength/100., .8, 10) // connect GC[i] to BC[j],syn[2]+dendritic_var,wt,del,threshold - if (debug_ == 0 ) { print "GC \t",i,"\t to BC\t\t",npost,"\t",dbr+2} - vgc2bc.x[npost] +=1 // increment the no of synapses to the post cell - } else { - j -= 1 - if (debug_ == 1 ) {print "nogc2bc"} - } // for connection that is not made reconnect axon to another cell - } - - for j=0, 0 { - // Based on the lamellar distribution of the GCs to MCs - 500 GCs were divided into 5 groups, 3 MCs were distributed in each lamella - // print "Based on the lamellar distribution of the GCs to MCs..." - if (i < ngcell/5) { a=0} - if ((i > ngcell/5-1) && (i < ngcell*2/5)) { a=1} - if ((i > ngcell*2/5-1) && (i < ngcell*3/5)) { a=2} - if ((i > ngcell*3/5-1) && (i < ngcell*4/5)) { a=3} - if ((i > ngcell*4/5-1) && (i < ngcell)) { a=4} - b=a*3 // from [0:12] - npost = rdgc2mc.repick() // from [0:2] - dbr = rdsynb.repick() // from [0:2] - //print npost, b - - - if (vgc2mc.x[npost+b] < 38){ - nc_append(i, ngcell+nbcell+npost+b, dbr+4, 0.2e-3, 1.5, 10) // Gcell[3] to Bcell[1] - if (debug_ == 1 ) {print "GC \t",i,"\t to MC\t\t", npost+b, "\t", dbr+4} - vgc2mc.x[npost+b] +=1 - } else { - j -= 1 - if (debug_ == 1 ) {print "nogc2mc"} - } - } - - - for j=0, 2 { - // comment added: - // Based on the lamellar distribution of the GCs to HIPPs - 500 GCs were divided into 6 groups - // print "Based on the lamellar distribution of the GCs to HIPPs..." - if (i < ngcell/6) { a=0} - if ((i > ngcell/6-1) && (i < ngcell*2/6)) { a=1} - if ((i > ngcell*2/6-1) && (i < ngcell*3/6)) { a=2} - if ((i > ngcell*3/6-1) && (i < ngcell*4/6)) { a=3} - if ((i > ngcell*4/6-1) && (i < ngcell*5/6)) { a=4} - if ((i > ngcell*5/6-1) && (i < ngcell)) { a=5} - - Gauz3 = rdgc2hc.repick() - if (a+Gauz3 > nhcell-1) {npost = a+Gauz3-nhcell } // change here to allow more then 6 HIPP cells [-2:2] - if (a+Gauz3 < 0) {npost = a+Gauz3+nhcell} - if ((a+Gauz3 > -1) && (a+Gauz3 < nhcell)) {npost = a+Gauz3} - dbr = rdsynb.repick() - if ((is_connected(HIPPCell[npost], GranuleCell[i]) == 0) && (vgc2hc.x[npost] < 275)) { - nc_append(i, ngcell+nbcell+nmcell+npost, dbr, 0.5e-3, 1.5, 10) // Bcell -> Hcell - if (debug_ == 0 ) {print "GC \t",i,"\t to HC\t\t",npost, "\t", dbr} - vgc2hc.x[npost] +=1 - } else { - j -= 1 - } - } - - // NOTE: THIS IS FOR SPROUTED SYNAPSES - n_sprout_ = p_sprouted_ -1 // NOTE: 100% Sprouting = 100 new synapses! (compare p. 443 in Santhakumar paper) - - for j=0,n_sprout_ { // 9 in original file // each GC is recurrent connected to 10 GC (i.e. 10/500 => p = 0.02) but sprouting is diff different -> see above - Gauz3 = rdgc2gc.repick() - //print Gauz3 - if (i+Gauz3 > 499) {npost = i+Gauz3-500 } - if (i+Gauz3 < 0) {npost = i+Gauz3+500} - if ((i+Gauz3 > -1) && (i+Gauz3 < 500)) {npost = i+Gauz3} - //print npost - dbr = rdsyna.repick() // [0:1] - if ((is_connected(GranuleCell[npost], GranuleCell[i]) == 0) && (vgc2gc.x[npost] < (n_sprout_*1.5+2) )) { // if is connected AND not more than 14 incoming connections... - // (original file < 15) (assume 1.5 times average connections for upper limit...) - nc_append(i, npost, dbr+7, 2e-3, .8, 10) // Gcell[3] to Bcell[1] - // print npost, dbr+8 - if (debug_ == 0 ) {print "GC \t",i,"\t to GC\t\t",npost, "\t", dbr+7} - vgc2gc.x[npost] +=1 - } else { - j -= 1 - if (debug_ == 0) {print "gc2gc"} - } - } -} - -if (print_GC_sprouting_input_ == 1) { - // objref distri_gc_input_ - // distri_gc_input_ = new Vector() - max_gc_input_ = 0 - if (debug_ ==2) { print "Calculate GC-GC Input Distribution"} - - for zz = 0, int(n_sprout_*1.5+2) {distri_gc_input_.append(0)} - for npost=0,ngcell-1 { - distri_gc_input_.x[vgc2gc.x[npost]]+=1 - if (vgc2gc.x[npost]>max_gc_input_) { max_gc_input_ = vgc2gc.x[npost]} // find max input number - } - for zz = 0, int(n_sprout_*1.5+2) {print zz,"\t",distri_gc_input_.x[zz]} - print "maximum input number is:\t",max_gc_input_ -} - - - -//****************************************************************************************** -} - -proc initNet_BC() { local i,j - -//**************Basket Cell post synaptic connections ****************************** - -if (debug_ ==2) {print "Basket Cell post synaptic connections ... "} - -for i=0, nbcell-1 { - for j=0, 99 { - Gauz3 = rdbc2gc.repick() // [-70:70] - if (debug_ == 1 ) {print Gauz3} - if (i*83+41+Gauz3 > 499) {npost = i*83+41+Gauz3-500 } - if (i*83+41+Gauz3 < 0) {npost = i*83+41+Gauz3+500} - if ((i*83+41+Gauz3 > -1) && (i*83+41+Gauz3 < 500)) {npost = i*83+41+Gauz3} - if (debug_ == 1 ) {print i, npost} - if (nbcell != 6) {max_conn = 4} else {max_conn = 2} // if not original setup use more spread... - if ((is_connected(GranuleCell[npost], BasketCell[i]) == 0) && (vbc2gc.x[npost] < max_conn)) { // change < 2 to < 4 - nc_append(i+ngcell, npost, 6, 1.6e-3*scale_BC2GC_strength/100, .85, -10) - vbc2gc.x[npost] +=1 - if (debug_ == 1 ) {print i, npost, 6} - } else { - j -= 1 - if (debug_ == 1 ) {print "BC2GC"} - } - } - - for j=0, 1 { - Gauz3 = rdbc2bc.repick() // [-1,0,1] (postsyn spread around single id...) - //print Gauz3 - if (i+Gauz3 > nbcell-1) {npost = i+Gauz3-nbcell } // periodic boundary conditions - if (i+Gauz3 < 0) {npost = i+Gauz3+nbcell} - if ((i+Gauz3 >-1) && (i+Gauz3 < nbcell)) {npost = i+Gauz3} - dbr = rdsyna.repick() // [0:1] - if (nbcell != 6) {max_conn = 4} else {max_conn = 3} - if ((is_connected(BasketCell[npost], BasketCell[i]) == 0) && (vbc2bc.x[npost] < max_conn)) { // change < 3 to < 4 - nc_append(i+ngcell, npost+ngcell, dbr+8, 7.6e-3, .8, -10) - if (debug_ == 1 ) {print npost, dbr+8} - vbc2bc.x[npost] +=1 - } else { - j -= 1 - if (debug_ == 1 ) {print "bc2bc"} - } - } - - for j=0, 2 { - Gauz3 = rdbc2mc.repick() // [-3:3] - if (i*2+2+Gauz3 > 14) {npost = i*2+2+Gauz3-15 } - if (i*2+2+Gauz3 < 0) {npost = i*2+2+Gauz3+15} - if ((i*2+2+Gauz3 >-1) && (i*2+2+Gauz3 < 15)) {npost = i*2+2+Gauz3} -// if ((is_connected(MossyCell[npost], BasketCell[i]) == 0) && (vbc2mc.x[npost] < 3) && (killMC.contains(ngcell+nbcell+npost) == 0)) { // use if killing MC - if (nbcell != 6) {max_conn = 4} else {max_conn = 3} - if ((is_connected(MossyCell[npost], BasketCell[i]) == 0) && (vbc2mc.x[npost] < max_conn)) { // change < 3 to < 4 - nc_append(i+ngcell, npost+ngcell+nbcell, 12, 1.5e-3, 1.5, -10) // Gcell[3] to Bcell[1] - if (debug_ == 1 ) {print npost, 12} - vbc2mc.x[npost] +=1 - } else { - j -= 1 - if (debug_ == 1 ) {print "bc2mc"} - } -// if (killMC.contains(ngcell+nbcell+npost) == 1) {j +=1 if (debug_ == 1 ) {print "dead MC"}} // use if killing MC - } - - -} -//****************************************************************************************** -} - - -proc initNet_rest() { local i,j - - - -//**************Mossy Cell post synaptic connections ****************************** - -if (debug_ ==2 ) { print "Mossy Cell post synaptic connections"} -for i=0, nmcell-1 { - // MC -> GC - //if (killMC.contains(ngcell+nbcell+i) == 0) // use if killing MC - - if (i < 3) { y=0} - if ((i > 2) && (i < 6)) { y=1} - if ((i > 5) && (i < 9)) { y=2} - if ((i > 8) && (i < 12)) { y=3} - if ((i > 11) && (i < 15)) { y=4} - - for j=0, 99 { - Gauz1 = rdmc2gc1.repick() // [25:175] - //if (debug_ == 1 ) {print Gauz1} - - if (i*33+17+Gauz1 > 499) { - npost1 = i*33+17+Gauz1-500 - } else {npost1 =i*33+17+Gauz1} - //if (debug_ == 1 ) {print npost1} - dbr = rdsyna.repick() // [0:1] - if ((is_connected(GranuleCell[npost1], MossyCell[i]) == 0) && (vmc2gc.x[npost1] < 7)) { - nc_append(i+ngcell+nbcell, npost1, dbr+2, 0.3e-3*scale_MC2GC_strength/100., 3, 10) // Gcell[3] to Bcell[1] - vmc2gc.x[npost1] +=1 - } else { j -= 1 } // if (debug_ == 1 ) {print "MC2GC1"} - } - - - for j=0, 99 { - Gauz2 = rdmc2gc2.repick() // [-175:25] - if (i*33+17+Gauz2 < 0) { - npost2 =i*33+17+Gauz2+500 - } else {npost2 =i*33+17+Gauz2} - dbr = rdsyna.repick() // [0:1] - if ((is_connected(GranuleCell[npost2], MossyCell[i]) == 0) && (vmc2gc.x[npost2] < 7)) { - nc_append(i+ngcell+nbcell, npost2, dbr+2, 0.3e-3*scale_MC2GC_strength/100., 3, 10) // Gcell[3] to Bcell[1] - vmc2gc.x[npost2] +=1 - } else { j -= 1 } - // if (debug_ == 1 ) {print "MC2GC2"} - } - - - // MC -> BC - for j=0, 0 { - Gauz3 = rdmc2bc.repick() // Gauz3 = [-3:3] - if (y+Gauz3 > nbcell-1) {npost = y+Gauz3-nbcell} // y = [0:4] => y+Gaus3 = [-3:7] - if (y+Gauz3 < 0) {npost = y+Gauz3+nbcell} - if ((y+Gauz3 > -1) && (y+Gauz3 < nbcell)) {npost = y+Gauz3} - dbr = rdsyna.repick() - if ((vmc2bc.x[npost] < 4) && (Gauz3 !=0)) { - nc_append(i+ngcell+nbcell, ngcell+npost, dbr+6, 0.3e-3, 3, 10) // Gcell[3] to Bcell[1] - vmc2bc.x[npost] += 1 - } else { - j -= 1 - // if (debug_ == 1 ) {print "mc2bc"} - } - } - - - // MC -> MC - for j=0, 2 { - Gauz3 = rdmc2mc1.repick() //[-3:3] - if (i+Gauz3 > 14) {npost = i+Gauz3-15 } - if (i+Gauz3 < 0) {npost = i+Gauz3+15} - if ((i+Gauz3 >-1) && (i+Gauz3 < 15)) {npost = i+Gauz3} - dbr = rdsynb.repick() -// if ((is_connected(MossyCell[npost], MossyCell[i]) == 0) && (vmc2mc.x[npost] < 4) && (Gauz3 != 0) && (killMC.contains(ngcell+nbcell+npost) == 0)) // use if killing MC - if ((is_connected(MossyCell[npost], MossyCell[i]) == 0) && (vmc2mc.x[npost] < 4) && (Gauz3 != 0)) { - nc_append(i+ngcell+nbcell, npost+ngcell+nbcell, dbr+8, 0.5e-3, 2, 10) // Gcell[3] to Bcell[1] -// print npost, dbr+8 - vmc2mc.x[npost] +=1 - } else {j -= 1 }// if (debug_ == 1 ) {print "mc2mc"} -// if (killMC.contains(ngcell+nbcell+npost) == 1){ j += 1 print "dead MC"} // use if killing MC - } - - - // MC -> HC - for j=0, 1 { - Gauz3 = rdmc2hc.repick() // [-2:2] - if (y+Gauz3 > nhcell-1) {npost = y+Gauz3-nhcell} // changed code here to allow > 6 HCells - if (y+Gauz3 < 0) {npost = y+Gauz3+nhcell} - if ((y+Gauz3 > -1) && (y+Gauz3 < nhcell)) {npost = y+Gauz3} - dbr = rdsynb.repick() - if ((is_connected(HIPPCell[npost], MossyCell[i]) == 0) && (vmc2hc.x[npost] < 7) && (Gauz3 != 0)) { - nc_append(i+ngcell+nbcell, ngcell+nbcell+nmcell+npost, dbr+4, 0.2e-3, 3, 10) // Gcell[3] to Bcell[1] - // print npost, dbr+4 - vmc2hc.x[npost] +=1 - } else { - j -= 1 - // if (debug_ == 1 ) {print y, Gauz3, "mc2hc"} - } - } -} - - - - -//****************************************************************************************** -//**************HIPP Cell post synaptic connections ****************************** - - -if (debug_ ==2 ) { print "HIPP Cell post synaptic connections"} - for i=0, nhcell-1 { - for j=0, 159 { // NOTE number of connections explicitly coded here! - Gauz3 = rdhc2gc.repick() // [-130:130] - //print Gauz3 - if (i*83+41+Gauz3 > 499) {npost = i*83+41+Gauz3-500 } // NOTE: 500 is explicitly coded here! - if (i*83+41+Gauz3 < 0) {npost = i*83+41+Gauz3+500} - if ((i*83+41+Gauz3 > -1) && (i*83+41+Gauz3 < 500)) {npost = i*83+41+Gauz3} - //print npost - dbr = rdsyna.repick() - if (nhcell != 6) {max_conn = 5} else {max_conn = 3} - if ((is_connected(GranuleCell[npost], HIPPCell[i]) == 0) && (vhc2gc.x[npost] < max_conn)) { // NOTE < 3 coded explicitly here! -> change to 5 - nc_append(i+ngcell+nbcell+nmcell, npost, dbr+4, 0.5e-3*scale_HC2GC_strength/100., 1.6, 10) // Hcell to Gcell - vhc2gc.x[npost] +=1 - if (debug_ == 1 ) {print i, npost, dbr+4} - } else { - j -= 1 - if (debug_ == 1 ) {print "HC2GC"} - } - } - - for j=0, 3 { - Gauz3 = rdhc2bc.repick() // [-2:2] - if (i+Gauz3 > nbcell-1) {npost = i+Gauz3-nbcell} - if (i+Gauz3 < 0) {npost = i+Gauz3+nbcell} - if ((i+Gauz3 > -1) && (i+Gauz3 < nbcell)) {npost = i+Gauz3} - dbr = rdsyna.repick() - if ((is_connected(BasketCell[npost], HIPPCell[i]) == 0) && (vhc2bc.x[npost] < 5)) { // NOTE < 5 coded explicitly here! - nc_append(i+ngcell+nbcell+nmcell, npost+ngcell, dbr+10, 0.5e-3, 1.6, 10) // Hcell to Bcell - if (debug_ == 1 ) {print npost, dbr+10} - vhc2bc.x[npost] += 1 - } else { - j -= 1 - if (debug_ == 1 ) {print "HC2BC"} - } - } - - - for j=0, 3 { - Gauz3 = rdhc2mc.repick() // [-2:2] - //print Gauz3 - if (i*2+2+Gauz3 > 14) {npost = i*2+2+Gauz3-15 } - if (i*2+2+Gauz3 < 0) {npost = i*2+2+Gauz3+15} - if ((i*2+2+Gauz3 >-1) && (i*2+2+Gauz3 < 15)) {npost = i*2+2+Gauz3} - //print npost - dbr = rdsynb.repick() - // if ((is_connected(MossyCell[npost], HIPPCell[i]) == 0) && (vhc2mc.x[npost] < 2) && (killMC.contains(ngcell+nbcell+npost) == 0)) //use if killing MC - if (nhcell != 6) {max_conn = 4} else {max_conn = 2} - if ((is_connected(MossyCell[npost], HIPPCell[i]) == 0) && (vhc2mc.x[npost] < max_conn)) { // NOTE < 2 coded explicitly here!! => changed to 4 - nc_append(i+ngcell+nbcell+nmcell, npost+ngcell+nbcell, dbr+13, 1.5e-3, 1, 10) // Hcell to Mcell - //if (debug_ == 1 ) {print npost, dbr+13} - vhc2mc.x[npost] += 1 - } else { - j -= 1 - if (debug_ == 1 ) {print "HC2MC"} - } - // if (killMC.contains(ngcell+nbcell+npost) == 1){ j += 1 print "dead MC"} //use if killing MC - } - } - - if (debug_ == 2 ) {print "Connections drawn -> Start Simulation"} -} - -objref VmT // vector of time points -objref VmMat[cells.count] // for each cell vector of Membrane potential - - -VmT = new Vector() -for i=0, cells.count-1 { - VmMat[i] = new Vector() -} - -proc VecMx() { local i // procedure called in every time step to add voltage of recorded cell at time t - VmT.append(t) - for i=0, (ngcell+nbcell +nmcell +nhcell -1) { - VmMat[i].append( cells.object[i].soma.v(0.5)) - } -} - - -proc ConvVT2Sp() { local i, j, max_id - max_id = npp - if ( print_template ==1) { max_id = npp } - if ( print_Raster ==1) { max_id = ngcell+nbcell +nmcell +nhcell-1 } - - for i=0, (max_id) { - Spike[i].spikebin(VmMat[i], 0) // Used to make a binary version of a spike train. - // is a vector of membrane potential. - // is the voltage threshold for spike detection. - // is set to all zeros except at the onset of spikes (the first dt which the spike crosses threshold) - } -} -// ********* Print files ************************** -err_ = load_file("printfile.hoc") - -//################################################################################################ - -// Define Random Generator for each PPStim object... -print "Define Random Generator for each PPStim object..." - -// nslist = new List() -rslist = new List() -rs_noiselist = new List() - -// ####### Input Window (box) - random_stream_offset_ = PP_rate_ * 1000 - for i=0, npp-1 { - rs = new RandomStream(i) // (stream nr) - rslist.append(rs) - PPSt[i].pp.noiseFromRandom(rs.r) - rs.r.uniform(0,1) // use uniform distribution - rs.start() - } - for i=0, npp-1 { - rs = new RandomStream(i+npp) // (stream nr) - rs_noiselist.append(rs) - PPSt_noise[i].pp.noiseFromRandom(rs.r) - rs.r.uniform(0,1) // use uniform distribution - rs.start() - } - - aacell pp_start = new NetStim125(.5) // artificial puls to PPSt[i].pp in order to become active... - pp_start.interval = 1e-5 - pp_start.number = 1 - pp_start.start = 0 - pp_start.forcestop = 1. - pp_start.noise = 1 - - rs = new RandomStream(npp) // assign Random Generator to init.. - return_val_ = rslist.append(rs) // save returned value in return_val_ to suppress screen output - return_val_ = pp_start.noiseFromRandom(rs.r) - return_val_ = rs.r.negexp(1) - return_val_ = rs.start() - - for i=0, npp-1 { - netcon_d[i] = new NetCon( pp_start,PPSt[i].pp) // each PPSt[i].pp needs to receive a netevent to become active... - netcon_d[i].weight = 10. - netcon_d[i].delay = 0.001 - netcon_d[i].threshold = 10. - } -print "Finished Random Generator for each PPStim object" - - - -objref dataset, data_ -proc read_custominit() { local i - data_ = new Vector() - dataset = new File() - err_ = dataset.ropen(init_state_fname_) - if (err_ ==0) { - print "IO code:",err_ - print "Initial State File not found!" - quit() - } - err_ = data_.scanf(dataset) - err_ = dataset.close() - cnt = 0 - for i = 0, ngcell -1 { - forsec Gcell[i].all { // iterate over all sections of this cell - for (x,0) { // iterate over internal nodes - v = data_.x[cnt] // f() returns the desired initial potential - if (v>-30.) { v = -30. } // cut off initial potential -> no spike artifacts ... - cnt += 1 - } - } - } - for i = 0, nbcell -1 { - forsec Bcell[i].all { // iterate over all sections of this cell - for (x,0) { // iterate over internal nodes - v = data_.x[cnt] - cnt += 1 - } - } - } - for i = 0, nmcell -1 { - forsec Mcell[i].all { // iterate over all sections of this cell - for (x,0) { // iterate over internal nodes - v = data_.x[cnt] - cnt += 1 - } - } - } - for i = 0, nhcell -1 { - forsec Hcell[i].all { // iterate over all sections of this cell - for (x,0) { // iterate over internal nodes - v = data_.x[cnt] - cnt += 1 - } - } - } - finitialize() -} - - -proc init() { local dtsav, temp, secsav //localobj ns, rs // init with a int value =! 100 reset not all random number generators equally... - v_init = -77.71 //-70.45 - finitialize(v_init) - t = -1000 // negative t step to initialize to steady state - dtsav = dt // original simulation step size - dt= 10 // larger step size - // if cvode is on, turn it off to do large fixed step - temp= cvode.active() - if (temp!=0) {cvode.active(0)} - while(t<-100) { - fadvance() - } - //restore cvode if reqd - if (temp!=0) {cvode.active(1)} - dt = dtsav - t = 0 - if (cvode.active()){ - cvode.re_init() - }else{ - fcurrent() - } - - // restart number generators for PPStim - t = 0 - finitialize(v_init) - trial_old = trial - - if (debug_ ==2 ) { print "number of Poisson inputs:\t",rslist.count()} - if ($1 == 100) { - for i = 0, rslist.count()-1 rslist.o(i).start() - // for j=0, 5 print j, rslist.o(0).r.repick - if (debug_ ==2 ) {print "reset all Poisson inputs with the correct seed"} - } else { - trial = trial + 1 - for i = 0, int((rslist.count()-1)*(1-$1/100.)) rslist.o(i).start() - if (debug_ ==2 ) {print "reset ",int((rslist.count()-1)*(1-$1/100.))," Poisson inputs with a different seed"} - trial = trial_old - for i = int((rslist.count()-1)*(1-$1/100.)), rslist.count()-1 rslist.o(i).start() - if (debug_ ==2 ) {print "reset ", rslist.count()-int((rslist.count()-1)*(1-$1/100.))," Poisson inputs with the correct seed"} - } - - // init noise Generators (with seed trial_noise ...) - trial = trial_noise - for i = 0, rs_noiselist.count()-1 rs_noiselist.o(i).start() - trial = trial_old - if (debug_ ==2 ) {print "reset all Poisson noise inputs with the trial_noise seed"} - - - VmT = new Vector() - for i=0, cells.count-1 { - VmMat[i] = new Vector() - } - - for i=0, (ngcell+nbcell +nmcell +nhcell -1) { - Spike[i] = new Vector() // vector of spikes for each cell - Spike_times[i] = new Vector() // vector of spikes for each cell - } - - t = 0 - finitialize(v_init) // finalize at the end - - if (init_from_file_ ==1 ) { - if (debug_ == 2) {print "Read in Init State of Network"} - read_custominit() // read in initial states from file - } - // custominit() - - frecord_init() - - - - - // ### select input pattern.... - for i=PP_box_startid_,PP_box_startid_+PP_box_nr_-1 { - PPSt[i].pp.status = 1 - PPSt[i].pp.start = PP_box_start_ - PPSt[i].pp.forcestop = PP_box_stop_ - PPSt[i].pp.nspk = PP_box_nspk_ - } -} - -proc run(){ - initNet_GC() - initNet_BC() - initNet_rest() - for (PP_box_startid_=0; PP_box_startid_<=12;PP_box_startid_+=1) { - saveNet() - init(100) - if (print_Vtrace == 1) { sMatrix_init() } // file header voltage output file - while (t - ek (mV) + ek (mV) gl (mho/cm2) : leak (gbar and reversal poti) el (mV) @@ -115,11 +115,11 @@ ASSIGNED { dt (ms) : 2) - gna (mho/cm2) : Na - ina (mA/cm2) + gna (mho/cm2) : Na + ina (mA/cm2) - gkf (mho/cm2) : K - gks (mho/cm2) + gkf (mho/cm2) : K + gks (mho/cm2) ik (mA/cm2) il (mA/cm2) : leak @@ -135,12 +135,12 @@ ASSIGNED { : This block is evaluated every time step. BREAKPOINT { SOLVE states : here the state variables are updated - gna = gnatbar*m*m*m*h : calculated g at timepoint t - gkf = gkfbar*nf*nf*nf*nf - gks = gksbar*ns*ns*ns*ns + gna = gnatbar*m*m*m*h : calculated g at timepoint t + gkf = gkfbar*nf*nf*nf*nf + gks = gksbar*ns*ns*ns*ns - ina = gna*(v - ena) : calculated currents flowing - ik = gkf*(v-ek) + gks*(v-ek) + ina = gna*(v - ena) : calculated currents flowing + ik = gkf*(v-ek) + gks*(v-ek) il = gl*(v-el) igabaa = ggabaa*(v-egabaa) } @@ -153,12 +153,8 @@ INITIAL { m = minf h = hinf - nf = nfinf - ns = nsinf - - VERBATIM - return 0; - ENDVERBATIM + nf = nfinf + ns = nsinf } : discreticed versions of the differential equations, hence a PROCEDURE and not DERIVATIVE block @@ -168,9 +164,6 @@ PROCEDURE states() { : Computes state variables m, h, and n h = h + hexp*(hinf-h) nf = nf + nfexp*(nfinf-nf) ns = ns + nsexp*(nsinf-ns) - VERBATIM - return 0; - ENDVERBATIM } : moved this to assign block @@ -178,46 +171,51 @@ PROCEDURE states() { : Computes state variables m, h, and n :Computes rate and other constants at current v. PROCEDURE rates(v) { - LOCAL alpha, beta, sum - q10 = 3^((celsius - 6.3)/10) - :"m" sodium activation system - act and inact cross at -40 : shifted by 68mV compared to in Aradi 1999/2002 + LOCAL alpha, beta, sum + q10 = 3^((celsius - 6.3)/10) + + :"m" sodium activation system - act and inact cross at -40 : shifted by 68mV compared to in Aradi 1999/2002 alpha = -0.3*vtrap((v+60-17),-5) : in Aradi 1999: alpha = -0.3*vtrap((v-25),-5); in Aradi 2002: alpha = 0.3*vtrap((v-25),-5) beta = 0.3*vtrap((v+60-45),5) : in Aradi 1999: beta = 0.3*vtrap((v-53),5); in Aradi 2002: beta = -0.3*vtrap((v-53),5) sum = alpha+beta - mtau = 1/sum minf = alpha/sum - :"h" sodium inactivation system : shifted by 68mV compared to in Aradi 1999/2002 + mtau = 1/sum + minf = alpha/sum + + :"h" sodium inactivation system : shifted by 68mV compared to in Aradi 1999/2002 alpha = 0.23/exp((v+60+5)/20) : in Aradi 1999/2002: alpha = 0.23/exp((v-3)/20) - beta = 3.33/(1+exp((v+60-47.5)/-10)) : in Aradi 1999/2002: beta = 3.33/(1+exp((v-55.5)/-10)) + beta = 3.33/(1+exp((v+60-47.5)/-10)) : in Aradi 1999/2002: beta = 3.33/(1+exp((v-55.5)/-10)) sum = alpha+beta htau = 1/sum - hinf = alpha/sum - + hinf = alpha/sum - :"ns" sKDR activation system : shifted by 65mV compared to Aradi 1999 - alpha = -0.028*vtrap((v+65-35),-6) : in Aradi 1999: alpha = -0.028*vtrap((v-35),-6) + :"ns" sKDR activation system : shifted by 65mV compared to Aradi 1999 + alpha = -0.028*vtrap((v+65-35),-6) : in Aradi 1999: alpha = -0.028*vtrap((v-35),-6) beta = 0.1056/exp((v+65-10)/40) : in Aradi 1999: beta = 0.1056/exp((v-10)/40) sum = alpha+beta - nstau = 1/sum nsinf = alpha/sum - :"nf" fKDR activation system : shifted by 65mV compared to Aradi 1999/2002 - alpha = -0.07*vtrap((v+65-47),-6) : in Aradi 1999: alpha = -0.07*vtrap((v-47),-6); in Aradi 2002: alpha = 0.07*vtrap((v-47),-6) + nstau = 1/sum + nsinf = alpha/sum + + :"nf" fKDR activation system : shifted by 65mV compared to Aradi 1999/2002 + alpha = -0.07*vtrap((v+65-47),-6) : in Aradi 1999: alpha = -0.07*vtrap((v-47),-6); in Aradi 2002: alpha = 0.07*vtrap((v-47),-6) beta = 0.264/exp((v+65-22)/40) : in Aradi 1999/2002: beta = 0.264/exp((v-22)/40) // probably typo in Aradi & Soltez 2002 there: beta = 0.264/exp((v-22)/4) sum = alpha+beta - nftau = 1/sum nfinf = alpha/sum + nftau = 1/sum + nfinf = alpha/sum } : Computes rate and other constants at current v. PROCEDURE trates(v) { LOCAL tinc - : TABLE minf, mexp, hinf, hexp, nfinf, nfexp, nsinf, nsexp, mtau, htau, nftau, nstau : + : TABLE minf, mexp, hinf, hexp, nfinf, nfexp, nsinf, nsexp, mtau, htau, nftau, nstau : : DEPEND dt, celsius FROM -100 TO 100 WITH 200 : rates(v) : not consistently executed from here if usetable_hh == 1 : so don't expect the tau values to be tracking along with : the inf values in hoc - tinc = -dt * q10 - mexp = 1 - exp(tinc/mtau) - hexp = 1 - exp(tinc/htau) + tinc = -dt * q10 + mexp = 1 - exp(tinc/mtau) + hexp = 1 - exp(tinc/htau) nfexp = 1 - exp(tinc/nftau) nsexp = 1 - exp(tinc/nstau) } diff --git a/mods/km.mod b/mods/km.mod new file mode 100644 index 0000000..e9a6415 --- /dev/null +++ b/mods/km.mod @@ -0,0 +1,92 @@ +TITLE CA1 KM channel from Mala Shah +: M. Migliore June 2006 + +UNITS { + (mA) = (milliamp) + (mV) = (millivolt) + +} + +PARAMETER { + v (mV) + ek + celsius (degC) + gbar=.0001 (mho/cm2) + vhalfl=-40 (mV) + kl=-10 + vhalft=-42 (mV) + a0t=0.003 (/ms) + zetat=7 (1) + gmt=.4 (1) + q10=5 + b0=60 + st=1 + sh =24 +} + + +NEURON { + SUFFIX km + USEION k READ ek WRITE ik + RANGE gbar,ik, sh + GLOBAL inf, tau +} + +STATE { + m +} + +ASSIGNED { + ik (mA/cm2) + inf + tau + taua + taub +} + +INITIAL { + rate(v) + m=inf +} + + +BREAKPOINT { + SOLVE state METHOD cnexp + ik = gbar*m^st*(v-ek) +} + + +FUNCTION alpt(v(mV)) { + alpt = exp(0.0378*zetat*(v-vhalft-sh)) +} + +FUNCTION bett(v(mV)) { + bett = exp(0.0378*zetat*gmt*(v-vhalft-sh)) +} + +DERIVATIVE state { + rate(v) + m' = (inf - m)/tau +} + +PROCEDURE rate(v (mV)) { :callable from hoc + LOCAL a,qt + qt=q10^((celsius-35)/10) + inf = (1/(1 + exp((v-vhalfl-sh)/kl))) + a = alpt(v) + tau = b0 + bett(v)/(a0t*(1+a)) +} + + + + + + + + + + + + + + diff --git a/netstim125.mod b/mods/netstim125.mod similarity index 100% rename from netstim125.mod rename to mods/netstim125.mod diff --git a/netstimbox.mod b/mods/netstimbox.mod similarity index 97% rename from netstimbox.mod rename to mods/netstimbox.mod index fbfb27a..2b927e0 100644 --- a/netstimbox.mod +++ b/mods/netstimbox.mod @@ -86,7 +86,7 @@ NET_RECEIVE (w) { if (event < 0) { event = 0 } - printf ("NetStimBox: Send spike at: %f\n",event) + : printf ("NetStimBox: Send spike at: %f\n",event) net_event(event) } status = 0 : switch it off diff --git a/BC.hoc b/objects/BC.hoc similarity index 99% rename from BC.hoc rename to objects/BC.hoc index 0beaed8..14fdfd7 100644 --- a/BC.hoc +++ b/objects/BC.hoc @@ -10,11 +10,10 @@ // http://onlinelibrary.wiley.com/doi/10.1002/hipo.22373/abstract // modified by +// Abraham Nunes / 2022 // Man Yi Yim / 2015 // Alexander Hanuschkin / 2011 -objref Bcell[nbcell] - begintemplate BasketCell ndend1=4 diff --git a/objects/DentateGyrus.hoc b/objects/DentateGyrus.hoc new file mode 100644 index 0000000..27cfa74 --- /dev/null +++ b/objects/DentateGyrus.hoc @@ -0,0 +1,1479 @@ +/* DENTATE GYRUS NETWORK TEMPLATE + TODO: + - Consider allowing seed change via argument to init() + - Change the scaling to not be from 0-100 but from 0-1 + - Account for different networks in file output names + - Implement PP stimulation of MCs + - Might rename `global_seed` + + - Where to put the following? + // for debug purpose only... + objref netcon, nil + objref vec_debug_ + objref PP_init_rnd // random number generator to select an PP stimulus site.. + objref nslist, rslist, rs_noiselist + objref ns, rs + + - Update docstrings on nc_append_rec() and nc_append() + - Clean the following procs/funcs, which were not changed much from Yim et al.: + - onnect_pp_to_gc() + - connect_pp_to_bc() + - connect_pp_to_hipp() + + - Address the `KillMC` aspects that are not integrated in Yim, but from Santhakumar + + - Functions to be renamed for consistent style: + - nc_append() + - nc_append_rec() + + - Variables to be renamed for consistent style: + - PP_box_startid_ + - PP_box_nr_ + - DGVt_name_ + - print_Vtrace + - print_Raster +*/ + + +// Load ranstream +err_ = load_file("ranstream.hoc") + +begintemplate DentateGyrus + // Public methods + public init, run + + // Simulation related declarations + // Random number generators + objref rdgc2gc, rdgc2bc, rdgc2mc, rdgc2hc + objref rdbc2gc, rdbc2bc, rdbc2mc + objref rdmc2gc1, rdmc2gc2, rdmc2bc, rdmc2mc, rdmc2mc1, rdmc2hc + objref rdhc2gc, rdhc2bc, rdhc2mc + objref rnd_pp2gc, rnd_pp2bc, rnd_pp2hc + objref rdsynb, rdsyna + objref rdsa + + // Declare cvode + objref cvode + + // Data management related declarations + public VmT, VmMat, idname + objref Spike[726] // Binary representations of spikes + objref Spike_times[726] // Spike times for each cell + objref VmT // Vector storing time index + objref VmMat[726] // One vector per cell storing membrane potential + objref efile // Membrane potential storage file + + // Declarations related to debugging + objref vec_stim[100], vec_stim_noise[500] + objref distri_gc_input_ + + // Network related declarations + public ngcell, nbcell, nmcell, nhcell, npp, nmcell_per_lamella, ntotal + public p_sprouted_, scale_gpas_dg_, scale_sk_dg_, scale_kir_, scale_gabaa_ + public scale_PP_strength, scale_PP2BC_strength, scale_PP2HC_strength + public scale_HC2GC_strength, scale_MC2GC_strength, scale_GC2BC_strength + public scale_BC2GC_strength + public cells, nclist + objref cells , Gcell[500], Bcell[6], Mcell[15], Hcell[6] + objref PPSt[100], PPSt_noise[500] + objref rs, rs_noise, rslist, rs_noiselist, pp_start, pp_noise_start + objref nclist, netcon, netcon_d[100], netcon_d_noise[500] + objref vgc2gc, vgc2bc, vgc2mc, vgc2hc + objref vbc2gc, vbc2bc, vbc2mc + objref vmc2gc, vmc2bc, vmc2mc, vmc2hc + objref vhc2gc, vhc2bc, vhc2mc + objref pp2gc, pp2bc, pp2hc + create aacell // artificial cell if windows (boxed) activation [ TODO ] - Study this more + + proc init() { + /* Initialization of the DentateGyrus object + + Arguments: + $s1 : str : A label for this simulation + $2 : int : Random state (global RNG seed) + $s3 : str : Name of the group being simulated + */ + strdef idname, group + idname = $s1 + random_state = $2 + group = $s3 + + // Specify directory to save data + // [ TODO ] - Allow this to be changed optionally + // without having to specify it + // every time. I.e. need to look at + // how one can specify default + // arguments to functions in NEURON + strdef datadir + datadir = "data/dgnetwork/" + + // Simulation and data storage parameters + set_simulation_parameters(random_state) + set_data_management_params() + + // Set parameters related to network + set_neuron_population_params() + set_perforant_path_input_params() + set_connectivity_params() + + // Initialize pseudorandom number generators + set_pseudorandom_number_generators() + + // Parameter modifications + bd_param_settings(group) + //yim_param_modifications(fig) + + } + + // ######################################################################## + // FUNCTIONS FOR SETTING PARAMETERS, HOUSEKEEPING, ETC + // ######################################################################## + proc set_simulation_parameters() { + /* Sets parameters for numerical aspects of the simulation + + Arguments: + $1 : int : Global RNG seed + */ + n_patterns = 12 // number of patterns to simulate (should be 12 for the full sim) + dt = 0.1 // simulation time step [ms] + tstop = 200 // total simulation time [ms] + trial = $1 // trialnumber = seed of simulation [ GENERALIZE ] + trial_noise = 1 // seed for noise to PP + rseed = $1 // pseudorandom number generator seed + debug_ = 2 + + spontaneous_activity_ = 1 // whether the GC neuron is spontaneously active + test_sa = 0 // 1 if we are testing spontaneous GC activity (cuts all inputs to GCs) + + } + + proc set_pseudorandom_number_generators(){ + /* Create the pseudorandom number generators for the network + + TODO: + - Pass argument to ranstream rather than requiring external + */ + + // Create RNGs for network connectivity establishment + // GC -> {GC, BC, MC, HC} + rdgc2gc = new Random(rseed) + rdgc2bc = new Random(rseed) + rdgc2mc = new Random(rseed) + rdgc2hc = new Random(rseed) + rdgc2gc.discunif(-50, 50) + rdgc2bc.discunif(-1,1) + rdgc2mc.discunif(0,2) + rdgc2hc.discunif(-2,2) + + // BC -> {GC, BC, MC} + rdbc2gc = new Random(rseed) + rdbc2bc = new Random(rseed) + rdbc2mc = new Random(rseed) + rdbc2gc.discunif(-70, 70) + rdbc2bc.discunif(-1, 1) + rdbc2mc.discunif(-3, 3) + + // MC -> {GC, BC, MC, HC} + rdmc2gc1 = new Random(rseed) + rdmc2gc2 = new Random(rseed) + rdmc2bc = new Random(rseed) + rdmc2mc = new Random(rseed) + rdmc2mc1 = new Random(rseed) + rdmc2hc = new Random(rseed) + rdmc2gc1.discunif(25, 175) + rdmc2gc2.discunif(-175, -25) + rdmc2bc.discunif(-3,3) + rdmc2mc.discunif(ngcell+nbcell, ngcell+nbcell+nmcell-1) + rdmc2mc1.discunif(-3, 3) + rdmc2hc.discunif(-2, 2) + + // HC -> {GC, BC, MC} + rdhc2gc = new Random(rseed) + rdhc2bc = new Random(rseed) + rdhc2mc = new Random(rseed) + rdhc2gc.discunif(-130, 130) + rdhc2bc.discunif(-2, 2) + rdhc2mc.discunif(-2, 2) + + // PP -> {GC, BC, HC} + rnd_pp2gc = new Random(rseed) + rnd_pp2bc = new Random(rseed) + rnd_pp2hc = new Random(rseed) + rnd_pp2gc.discunif(0,npp-1) + rnd_pp2bc.discunif(0,npp-1) + rnd_pp2hc.discunif(0,npp-1) + + // RNGs for synapses + rdsyna = new Random(rseed) + rdsynb = new Random(rseed) + rdsyna.discunif(0, 1) + rdsynb.discunif(0, 3) + + // RNG for spontaneous activity + rdsa = new Random(rseed) + rdsa.discunif(0, ngcell-1) + } + + proc set_data_management_params(){ + /* Parameters for how data are managed during the simulation. + */ + + // Define file names for simulation results output + strdef suffix + suffix = "txt" + + // Flags about what to print + print_Vtrace = 0 // print voltage trace to file + print_Raster = 1 // print spike raster to file + print_template = 1 // write out spike sequence into template file + print_GC_sprouting_input_ = 0 // calculates and print out GC sprouting input number distribution + print_stim_ = 1 // print stimulus to GCs for debug purpose.. + print_stim_noise = 1 // print noise to GCs for debug purpose.. + } + + proc set_neuron_population_params(){ + /* Initialize the number of cells of each type in the networks + + TODO: + - Generalize to allow these to vary + */ + ngcell = 500 // number of GCs + nbcell = 6 // number of BCs (Samathakumar 2005) + nmcell = 15 // number of MCs + nhcell = 6 // number of HCs (Samathakumar 2005) + npp = 100 // ECII Neurons (Myers and Scharfman, 2009) + nmcell_per_lamella = 3 // Number of mossy cells per lamella + + // Add 2*npp because of PP input and PP noise... + ntotal = ngcell + nbcell + nmcell + nhcell + npp + ngcell + + // Neuronal params + v_init = -77.71 + + } + + proc set_perforant_path_input_params() { + /* Set parameters for the perforant path stimulation + */ + PP_nstimulus_ = 1 // one input per GC + PP_input2MC_ = 0 // stimulate MC with PP + PP_rate_ = 10. // rate of PP Poisson input + PP_rate_noise_ = 0. // rate of PP Poisson noise + PP_box_nr_ = 6 // number of active PPs + PP_box_stop_ = 35. // time of box [ms] + PP_box_start_ = 5. // shift of box [ms] + PP_box_nspk_ = 3 // number of spike per active PP + + if (PP_input2MC_ != 0) { + print "PP stimulation not yet implemented if stimulation of also MCs...." + quit() + } + } + + proc set_connectivity_params() { + /* Parameters related to network connectivity + + Arguments: + + */ + p_sprouted_ = 0 // proportion of other GCs to which a GC will connect + scale_gpas_dg_ = 1 // scaling [%/100] of gpas of th GC model + scale_sk_dg_ = 1 // scaling of Ca-dependent K (SK) of th GC model + scale_kir_ = 1 // scaling of KIR conductance in GC model + scale_gabaa_ = 1 // scaling of GABAA in GC model + scale_PP_strength = 0.1*(1-test_sa) // scaling of synaptic weight of PP + scale_PP2BC_strength = 1*(1-test_sa) // scaling of synaptic weight PP->BC (100% for original value which was 50% of PP->GC strength..) + scale_PP2HC_strength = 0 // scaling of synaptic weight PP->HIPP (set to 0% because there are no PP->HIPP connections...) + scale_HC2GC_strength = 1*(1-test_sa) // scaling of synaptic weight HC->GC (beta_HIPP in Myers and Scharfman, 2009) + scale_MC2GC_strength = 1*(1-test_sa) // scaling of synaptic weight MC->GC (beta_HIPP in Myers and Scharfman, 2009) + scale_GC2BC_strength = 3*(1-test_sa) // scaling of synaptic weight GC->BC + scale_BC2GC_strength = 3*(1-test_sa) // scaling of synaptic weight BC->GC + spontaneous_activity_strength = 100 // scaling of the strength of synaptic weight for the spontaneous AP generators + + } + + proc yim_param_modifications() { + /* Reset parameters to generate the figures in Yim et al (2015) + + Arguments: + $1 : int : Which figure you are trying to generate data for + + */ + fig = $1 + if (fig == 1) { + scale_PP_strength = 0.10 + scale_kir_ = 1 + scale_gabaa_ = 1 + p_sprouted_ = 0 + } + + if (fig == 2) { + scale_PP_strength = 0.1 + scale_kir_ = 1 + scale_gabaa_ = 1 + p_sprouted_ = 0.3 + } + + if (fig == 3) { + scale_PP_strength = 0.1 + scale_kir_ = 4. + scale_gabaa_ = 4. + p_sprouted_ = 0.3 + } + + if (fig == 4) { + scale_PP_strength = 0.16 + scale_kir_ = 1. + scale_gabaa_ = 1. + p_sprouted_ = 0.1 + } + + if (fig == 5) { + scale_PP_strength = 0.16 + scale_kir_ = 4. + scale_gabaa_ = 4. + p_sprouted_ = 0.1 + } + } + + proc bd_param_settings() { + /* Sets different parameters for bipolar disorder related + experiments + + Arguments: + $s1 : str : Which group we are simulating + + Notes: + - We can change HT and LT later. + */ + + // Healthy control neurons + if (strcmp($s1, "HC") == 0) { + scale_na_conductances = 1 + scale_kdr_conductances = 1 + scale_ka_conductances = 1 + gbar_ht_ = 0.0004 + gbar_lt_ = 0 + scale_size_ = 1 + scale_gpas_dg_ = 1 + scale_sk_dg_ = 1 + scale_gabaa_ = 1 + scale_kir_ = 0 + spontaneous_activity_rate = 0.25 // In Hz + } + + // Lithium Responders + if (strcmp($s1, "LR") == 0) { + scale_na_conductances = 1.03 + scale_kdr_conductances = 1.15 + scale_ka_conductances = 0.0096/0.008 + gbar_ht_ = 0.0008 + gbar_lt_ = 0 + scale_size_ = 1 + scale_gpas_dg_ = 1 + scale_sk_dg_ = 1 + scale_gabaa_ = 1 + scale_kir_ = 0 + spontaneous_activity_rate = 1 // In Hz + } + + // Lithium Nonresponders + if (strcmp($s1, "NR") == 0) { + scale_na_conductances = 0.026/0.033 + scale_kdr_conductances = 1.3 + scale_ka_conductances = 0.0104/0.008 + gbar_ht_ = 0.0004 + gbar_lt_ = 0 + scale_size_ = 1 + scale_gpas_dg_ = 1 + scale_sk_dg_ = 1 + scale_gabaa_ = 1 + scale_kir_ = 0 + spontaneous_activity_rate = 0.25 // In Hz + } + } + + // ######################################################################## + // FUNCTIONS FOR SETTING UP CELLS/INPUTS + // ######################################################################## + proc make_cells() {local i + // Need to re-declare cell arrays with proper sizes + objref Gcell[ngcell], Bcell[nbcell], Mcell[nmcell], Hcell[nhcell] + objref PPSt[npp], PPSt_noise[ngcell] + cells = new List() + + // Create populations + for i=0, ngcell-1 {Gcell[i] = new GranuleCell(i,scale_na_conductances, scale_kdr_conductances, scale_ka_conductances, gbar_ht_, gbar_lt_, scale_size_, scale_gpas_dg_, scale_sk_dg_, scale_gabaa_, scale_kir_)} + for i=0, nbcell-1 {Bcell[i] = new BasketCell(i)} + for i=0, nmcell-1 {Mcell[i] = new MossyCell(i)} + for i=0, nhcell-1 {Hcell[i] = new HIPPCell(i)} + for i=0, npp-1 {PPSt[i] = new PPstim(i, PP_rate_, tstop, PP_box_start_, PP_box_stop_)} + + // [TESTING] Adding individual stims to each cell in order to simulate spontaneous activity + //for i=0, npp-1 {PPSt_noise[i] = new PPstim(i, PP_rate_, tstop, PP_box_start_, PP_box_stop_)} + for i=0, ngcell-1 { PPSt_noise[i] = new PPstim(i, PP_rate_, tstop, 1., tstop-1) } + + // Append to cells list + for i = 0, ngcell-1 {cells.append(Gcell[i])} // cells 0-499 GCs + for i = 0, nbcell-1 {cells.append(Bcell[i])} // cells 500-505 BC + for i = 0, nmcell-1 {cells.append(Mcell[i])} // cells 506-520 MC + for i = 0, nhcell-1 {cells.append(Hcell[i])} // cells 521-526 HC + for i = 0, npp-1 {cells.append(PPSt[i])} // 527 - xxx PP artificial cell + + // [TESTING] Append GC spontaneous activity generators to list + //for i = 0, npp-1 {cells.append(PPSt_noise[i])} + if (spontaneous_activity_ == 1) { + for i = 0, ngcell-1 { cells.append(PPSt_noise[i]) } + } + } + + proc make_input_stimulation() { local i + /* Creates the perforant path input stimulation objects + */ + // Print update to console + print "\tDefining Random Generator for each PPStim object." + objref pp_start, pp_noise_start + objref netcon_d[npp] + objref netcon_d_noise[ngcell] + + + rslist = new List() + rs_noiselist = new List() + + // Input window (box) + random_stream_offset_ = PP_rate_ * 1000 + for i=0, npp-1 { + rs = new RandomStream(i, random_state) + rslist.append(rs) + PPSt[i].pp.noiseFromRandom(rs.r) + rs.r.uniform(0,1) + rs.start(random_state) + } + + // artificial pulse to PPSt[i].pp in order to become active... + aacell pp_start = new NetStim125(.5) + pp_start.interval = 1e-5 + pp_start.number = 1 + pp_start.start = 0 + pp_start.forcestop = 1. + pp_start.noise = 1 + + // assign Random Generator to init.. + // [TODO] - Study more what is happening here + rs = new RandomStream(npp, random_state) + return_val_ = rslist.append(rs) // save returned value in return_val_ to suppress screen output + return_val_ = pp_start.noiseFromRandom(rs.r) + return_val_ = rs.r.negexp(1) + return_val_ = rs.start(random_state) + + // each PPSt[i].pp needs to receive a netevent to become active... + // [ TODO ] - allow for these parameters to be changed more generally + for i=0, npp-1 { + netcon_d[i] = new NetCon(pp_start,PPSt[i].pp) + netcon_d[i].weight = 10. + netcon_d[i].delay = 0.001 + netcon_d[i].threshold = 10. + } + + + if (spontaneous_activity_ == 1) { + // Input window (box) + for i=0, ngcell-1 { + rs_noise = new RandomStream(i, random_state) + rs_noiselist.append(rs_noise) + PPSt_noise[i].pp.noiseFromRandom(rs_noise.r) + rs_noise.r.uniform(0,1) + rs_noise.start(random_state) + } + + + // artificial pulse to PPSt_noise[i].pp in order to become active... + aacell pp_noise_start = new NetStim125(.5) + pp_noise_start.interval = 1e-5 + pp_noise_start.number = 1 + pp_noise_start.start = 0 + pp_noise_start.forcestop = 1. + pp_noise_start.noise = 1 + + // assign random number generator + rs_noise = new RandomStream(ngcell, random_state) + return_val_ = rslist.append(rs_noise) // save returned value in return_val_ to suppress screen output + return_val_ = pp_noise_start.noiseFromRandom(rs_noise.r) + return_val_ = rs_noise.r.negexp(1) + return_val_ = rs_noise.start(random_state) + + // each PPSt_noise[i].pp needs to receive a netevent to become active... + for i=0, ngcell-1 { + netcon_d_noise[i] = new NetCon(pp_noise_start,PPSt_noise[i].pp) + netcon_d_noise[i].weight = 10. + netcon_d_noise[i].delay = 0.001 + netcon_d_noise[i].threshold = 10. + } + } + + + } + + // ######################################################################## + // FUNCTIONS FOR MAKING NETWORK CONNECTIVITY + // ######################################################################## + proc make_connections(){ + /* Instantiates connectivity in the network + */ + // Instantiate the connectivity list + nclist = new List() + + // Instantiate vectors for connections + // GC -> {GC, BC, MC, HC} + vgc2gc = new Vector(ngcell, 0) + vgc2bc = new Vector(nbcell, 0) + vgc2mc = new Vector(nmcell, 0) + vgc2hc = new Vector(nhcell, 0) + + // BC -> {GC, BC, MC} + vbc2gc = new Vector(ngcell, 0) + vbc2bc = new Vector(nbcell, 0) + vbc2mc = new Vector(nmcell, 0) + + // MC -> {GC, BC, MC, HC} + vmc2gc = new Vector(ngcell, 0) + vmc2bc = new Vector(nbcell, 0) + vmc2mc = new Vector(nmcell, 0) + vmc2hc = new Vector(nhcell, 0) + + // HC -> {GC, BC, MC} + vhc2gc = new Vector(ngcell, 0) + vhc2bc = new Vector(nbcell, 0) + vhc2mc = new Vector(nmcell, 0) + + + // Connect areas + connect_pp_to_gc() + connect_pp_to_bc() + connect_pp_to_hipp() + connect_granule_cells() + connect_basket_cells() + connect_mossy_cells() + connect_hipp_cells() + + if (spontaneous_activity_ == 1) { + add_spontaneous_gc_activity() + } + } + + func nc_append() { + // neuron connect $1 with $2.pre_list.object($3), weight $4, delay $5, threshold $6 + // connects: + // cells.object($1) with + // $o1 = cells.object($2).pre_list.object($3) and + // returns: + // netcon = $o2 + + if ($3 >= 0 ) { + // connect_pre is function in the respective cell definition + cells.object($1).connect_pre(cells.object($2).pre_list.object($3),netcon) + netcon.weight = $4 + netcon.delay = $5 + netcon.threshold = $6 + } + nclist.append(netcon) + return nclist.count-1 + } + + func nc_append_rec() { + /* neuron connect $1 with $2.pre_list.object($3), weight $4, delay $5, threshold $6 + // connects: + // cells.object($1) with + // $o1 = cells.object($2).pre_list.object($3) and + // returns: + // netcon = $o2 + // record events to $o7 + */ + + if ($3 >= 0 ) { + // connect_pre is function in the respective cell definition + cells.object($1).connect_pre(cells.object($2).pre_list.object($3),netcon) + netcon.weight = $4 + netcon.delay = $5 + netcon.threshold = $6 + netcon.record($o7) + } + nclist.append(netcon) + return nclist.count-1 + } + + func is_connected() {local i, c localobj net_c + /* Checks for preexisting connections between randomly selected cells + to avoid multiple contacts between same 2 cells + */ + c=0 + for i=0, nclist.count-1 { + net_c = nclist.object(i) + if (($o1 == net_c.postcell()) && ($o2 == net_c.precell())) {c=1} + } + return c + } + + proc connect_pp_to_gc() { local i, j localobj pprec + /* Function that connects perforant path inputs to granule cells. + */ + // Print message to indicate procedure + print "Connecting PP to post-synaptic targets." + print "\tConnecting PP -> GC." + + // Create vector that marks PP neurons that already project to a GC + pprec = new Vector(npp, 0) + + // Objects to record inputs/noise for debugging [TODO] verify necessity... + // Must be re-declared with proper value of npp + objref vec_stim[npp] + for i = 0, npp-1 { vec_stim[i] = new Vector() } + + // Make connections for each individual GC + for i=0, ngcell-1 { + // Re-initialize connectivity vector to prevent multiple stimulation + // connections from the same PP cell + pp2gc = new Vector(npp, 0) + + // Each GC receives input from random set of 20% of PP neurons + // [ TODO ] - ALLOW TO HAVE 20% CHANGED TO DIFFERENT VALUES + for nr_conn = 0, int(npp/5.) - 1 { + // Pick a random PP cell to connect to the GC + // Repeat until j is a PP cell that is not yet connected + // Then mark j as connected so it is not chosen again + j = rnd_pp2gc.repick() + while (pp2gc.x[j] == 1) { + j = rnd_pp2gc.repick() + } + pp2gc.x[j] = 1 + + if ((print_stim_==1) && (i NOTE both synapses are equal, ie. weight delay (except position) not important + nc_append(j+ngcell+nbcell+nmcell+nhcell, i, 0, 2e-2*scale_PP_strength, 3, 10) + } + } + + // Print debugging message + if (debug_ == 3) { + print "\nGC: ",i + for ii=0,pp2gc.size()-1 { printf ("%d, ",pp2gc.x[ii]) } + } + } + } + + proc add_spontaneous_gc_activity() { local i,j + /* Function that adds spontaneous hyperactivity to the granule cells + */ + // Print message to indicate procedure + print "Adding spontaneous GC activity." + + // Objects to record inputs/noise for debugging [TODO] verify necessity... + // Must be re-declared with proper value of npp + objref vec_stim_noise[ngcell] + for i = 0, ngcell-1 { vec_stim_noise[i] = new Vector() } + + + // Add spontaneous activity generator for each individual GC + for i=0, ngcell-1 { + if ((print_stim_ == 1)) { + nc_append_rec(i+ngcell+nbcell+nmcell+nhcell+npp, i, 0, 2e-2*spontaneous_activity_strength, 0, 0, vec_stim_noise[i]) + } else { + nc_append(i+ngcell+nbcell+nmcell+nhcell+npp, i, 0, 2e-2*spontaneous_activity_strength, 0, 0) + } + } + + } + + proc connect_pp_to_bc() { local i,j + /* Connect perforant path inputs to basket cells + */ + // Print message to indicate procedure + print "\tConnecting PP -> BC." + + // Only proceed if there are PP to BC connections + if (scale_PP2BC_strength != 0) { + for i=ngcell, ngcell+nbcell-1 { + pp2bc = new Vector(npp, 0) + + // Each BC receives input from random set of 20% of PP neurons + // [ TODO ] - ALLOW TO HAVE 20% CHANGED TO DIFFERENT VALUES + for nr_conn = 0, int(npp/5.)-1 { + j = rnd_pp2bc.repick() + while (pp2bc.x[j] == 1) { j = rnd_pp2bc.repick() } + pp2bc.x[j] = 1 + + // connect PP to BC[j],syn[0],wt,del,threshold NOTE both synapses are equal + nc_append(j+ngcell+nbcell+nmcell+nhcell, i, 0, 1e-2*scale_PP_strength*scale_PP2BC_strength, 3, 10) + } + + // Print debugging message + if (debug_ == 3) { + print "\nBC: ",i + for ii=0,pp2bc.size()-1{ printf ("%d, ",pp2gc.x[ii])} + } + } + } + } + + proc connect_pp_to_hipp() { local i,j + /* Connecting perorant path to HIPP cells + */ + if (scale_PP2HC_strength != 0) { + print "\tConnecting PP -> HIPP." + + for i=0, nhcell-1 { + pp2hc = new Vector(npp, 0) + + // Each HIPP receives input from random set of 20% of PP neurons + // [ TODO ] - ALLOW TO HAVE 20% CHANGED TO DIFFERENT VALUES + for nr_conn = 0, int(npp/5.)-1 { + j = rnd_pp2hc.repick() + while (pp2hc.x[j] == 1) { j = rnd_pp2hc.repick() } + pp2hc.x[j] = 1 + + // used here synaptic strength of GC->HIPP / connections PP-> HIPP not in Santhakumar 2005 + nc_append(j+ngcell+nbcell+nmcell+nhcell, i+ngcell+nbcell+nmcell, 0, 0.5e-3*scale_PP_strength*scale_PP2HC_strength, 3, 10) + + } + if (debug_ == 3) { + print "\nHIPP: ",i + for ii=0,pp2bc.size()-1{ printf ("%d, ",pp2hc.x[ii])} + } + } + if (debug_ == 3) { print "PP -> HIPP -DONE-" } + } + } + + proc connect_granule_cells() { local i,j localobj c1, c2 + /* Connect granule cells to post-synaptic targets + */ + // Print to indicate step + print "Connecting GCs to post-synaptic targets." + + if (debug_ == 0) { print "n_sprout_ = ",n_sprout_ } + if (p_sprouted_>0 && debug_ == 2) {print "\tNOTE: Implementing mossy fibre sprouting."} + + for i=0, ngcell-1 { + // CONNECT GC'S TO BC'S + for j=0, 0 { + // Which lamella is GC[i] in wrt basket cells + // Note: one basket cell per lamella + // [ TODO ] - Allow for more than one BC per lamella + if (i < ngcell/nbcell) { a=0 } + if ((i > ngcell/nbcell-1) && (i < ngcell*2/nbcell)) { a=1 } + if ((i > ngcell*2/nbcell-1) && (i < ngcell*3/nbcell)) { a=2 } + if ((i > ngcell*3/nbcell-1) && (i < ngcell*4/nbcell)) { a=3 } + if ((i > ngcell*4/nbcell-1) && (i < ngcell*5/nbcell)) { a=4 } + if ((i > ngcell*5/nbcell-1) && (i < ngcell)) { a=5} + + // Randomly pick location of post synaptic Bcell from distribution [-1:1] + Gauz3 = rdgc2bc.repick() + + // Determine appropriate post syn BC + if (a+Gauz3 > nbcell-1) { npost = a+Gauz3-nbcell } + if (a+Gauz3 < 0) { npost = a+Gauz3+nbcell } + if ((a+Gauz3 > -1) && (a+Gauz3 < nbcell)) {npost = a+Gauz3} + + // randomly pick the dendrite to connect to from [0:3] ( i.e. // randomize among 4 dendrites ) + dbr = rdsynb.repick() + if (debug_ == 1 ) { + print "GC \t",i,"\t to BC\t\t",npost, a + } + + // check to make sure that post syn BC does not get more than 90 GC synapses + // [ TODO ] - ALLOW FOR DIFFERENT NUMBER OF MAX SYNAPSES + if (vgc2bc.x[npost] < 90) { + // connect GC[i] to BC[j],syn[2]+dendritic_var,wt,del,threshold + nc_append(i, ngcell+npost, dbr+2, 4.7e-3*scale_GC2BC_strength, .8, 10) + + if (debug_ == 0 ) { + print "GC \t",i,"\t to BC\t\t",npost,"\t",dbr+2 + } + + // increment the no of synapses to the post cell + vgc2bc.x[npost] +=1 + + } else { + j -= 1 + if (debug_ == 1 ) {print "nogc2bc"} + } // for connection that is not made reconnect axon to another cell + } + + // CONNECT GC'S TO MC'S + for j=0, 0 { + // Based on the lamellar distribution of the GCs to MCs - 500 GCs were divided into 5 groups, 3 MCs were distributed in each lamella + // print "Based on the lamellar distribution of the GCs to MCs..." + // [ TODO ] - Allow for different number of MC's per lamella + if (i < ngcell/5) { a=0} + if ((i > ngcell/5-1) && (i < ngcell*2/5)) { a=1} + if ((i > ngcell*2/5-1) && (i < ngcell*3/5)) { a=2} + if ((i > ngcell*3/5-1) && (i < ngcell*4/5)) { a=3} + if ((i > ngcell*4/5-1) && (i < ngcell)) { a=4} + b=a*3 // from [0:12] + npost = rdgc2mc.repick() // from [0:2] + dbr = rdsynb.repick() // from [0:2] + + + // [ TODO ] - ALLOW FOR DIFFERENT NUMBER OF MAX SYNAPSES + if (vgc2mc.x[npost+b] < 38){ + nc_append(i, ngcell+nbcell+npost+b, dbr+4, 0.2e-3, 1.5, 10) + + if (debug_ == 1 ) { + print "GC \t",i,"\t to MC\t\t", npost+b, "\t", dbr+4 + } + + vgc2mc.x[npost+b] +=1 + } else { + j -= 1 + if (debug_ == 1 ) { + print "nogc2mc" + } + } + } + + // CONNECT GC's to HIPP + for j=0, 2 { + // comment added: + // Based on the lamellar distribution of the GCs to HIPPs - 500 GCs were divided into 6 groups + if (i < ngcell/6) { a=0} + if ((i > ngcell/6-1) && (i < ngcell*2/6)) { a=1} + if ((i > ngcell*2/6-1) && (i < ngcell*3/6)) { a=2} + if ((i > ngcell*3/6-1) && (i < ngcell*4/6)) { a=3} + if ((i > ngcell*4/6-1) && (i < ngcell*5/6)) { a=4} + if ((i > ngcell*5/6-1) && (i < ngcell)) { a=5} + + Gauz3 = rdgc2hc.repick() + if (a+Gauz3 > nhcell-1) { npost = a+Gauz3-nhcell } + if (a+Gauz3 < 0) { npost = a+Gauz3+nhcell } + if ((a+Gauz3 > -1) && (a+Gauz3 < nhcell)) { npost = a+Gauz3 } + + dbr = rdsynb.repick() + c1 = cells.object(ngcell+nbcell+nmcell+npost) + c2 = cells.object(i) + if ((is_connected(c1, c2) == 0) && (vgc2hc.x[npost] < 275)) { + nc_append(i, ngcell+nbcell+nmcell+npost, dbr, 0.5e-3, 1.5, 10) + if (debug_ == 1 ) { + print "GC \t",i,"\t to HC\t\t",npost, "\t", dbr + } + vgc2hc.x[npost] +=1 + } else { + j -= 1 + } + } + + // GCs -> GCs + // NOTE: THIS IS FOR SPROUTED SYNAPSES + // NOTE: 100% Sprouting = 100 new synapses! (compare p. 443 in Santhakumar paper) + n_sprout_ = p_sprouted_ - 1 + + // 9 in original file // each GC is recurrent connected to 10 GC + // (i.e. 10/500 => p = 0.02) but sprouting is diff different -> see above + for j=0, n_sprout_ { + Gauz3 = rdgc2gc.repick() + if (i+Gauz3 > 499) { npost = i+Gauz3-500 } + if (i+Gauz3 < 0) { npost = i+Gauz3+500 } + if ((i+Gauz3 > -1) && (i+Gauz3 < 500)) { npost = i+Gauz3 } + + dbr = rdsyna.repick() // [0:1] + c1 = cells.object(npost) + c2 = cells.object(i) + if ((is_connected(c1, c2) == 0) && (vgc2gc.x[npost] < (n_sprout_*1.5+2) )) { // if is connected AND not more than 14 incoming connections... + // (original file < 15) (assume 1.5 times average connections for upper limit...) + nc_append(i, npost, dbr+7, 2e-3, .8, 10) // Gcell[3] to Bcell[1] + // print npost, dbr+8 + if (debug_ == 0 ) {print "GC \t",i,"\t to GC\t\t",npost, "\t", dbr+7} + vgc2gc.x[npost] +=1 + } else { + j -= 1 + if (debug_ == 0) {print "gc2gc"} + } + } + + } + + if (print_GC_sprouting_input_ == 1) { + distri_gc_input_ = new Vector() + max_gc_input_ = 0 + if (debug_ ==2) { print "Calculate GC-GC Input Distribution"} + + for zz = 0, int(n_sprout_*1.5+2) {distri_gc_input_.append(0)} + for npost=0,ngcell-1 { + distri_gc_input_.x[vgc2gc.x[npost]]+=1 + if (vgc2gc.x[npost]>max_gc_input_) { max_gc_input_ = vgc2gc.x[npost]} // find max input number + } + for zz = 0, int(n_sprout_*1.5+2) {print zz,"\t",distri_gc_input_.x[zz]} + print "maximum input number is:\t",max_gc_input_ + } + } + + proc connect_basket_cells() { local i,j localobj c1, c2 + print "Connecting BCs to post-synaptic targets. " + + for i=0, nbcell-1 { + // BC -> GC + // [ TODO ] - Allow for different lamella sizes + for j=0, 99 { + Gauz3 = rdbc2gc.repick() // [-70:70] + if (debug_ == 1 ) {print Gauz3} + if (i*83+41+Gauz3 > 499) {npost = i*83+41+Gauz3-500 } + if (i*83+41+Gauz3 < 0) {npost = i*83+41+Gauz3+500} + if ((i*83+41+Gauz3 > -1) && (i*83+41+Gauz3 < 500)) {npost = i*83+41+Gauz3} + if (debug_ == 1 ) {print i, npost} + if (nbcell != 6) {max_conn = 4} else {max_conn = 2} // if not original setup use more spread... + c1 = cells.object(npost) + c2 = cells.object(i+ngcell) + if ((is_connected(c1, c2) == 0) && (vbc2gc.x[npost] < max_conn)) { // change < 2 to < 4 + nc_append(i+ngcell, npost, 6, 1.6e-3*scale_BC2GC_strength, .85, -10) + vbc2gc.x[npost] +=1 + if (debug_ == 1 ) { print i, npost, 6 } + } else { + j -= 1 + if (debug_ == 1 ) {print "BC2GC"} + } + } + + // BC -> BC + for j=0, 1 { + Gauz3 = rdbc2bc.repick() // [-1,0,1] (postsyn spread around single id...) + //print Gauz3 + if (i+Gauz3 > nbcell-1) {npost = i+Gauz3-nbcell } // periodic boundary conditions + if (i+Gauz3 < 0) {npost = i+Gauz3+nbcell} + if ((i+Gauz3 >-1) && (i+Gauz3 < nbcell)) {npost = i+Gauz3} + dbr = rdsyna.repick() // [0:1] + if (nbcell != 6) {max_conn = 4} else {max_conn = 3} + c1 = cells.object(npost+ngcell) + c2 = cells.object(i+ngcell) + if ((is_connected(c1, c2) == 0) && (vbc2bc.x[npost] < max_conn)) { // change < 3 to < 4 + nc_append(i+ngcell, npost+ngcell, dbr+8, 7.6e-3, .8, -10) + if (debug_ == 1 ) {print npost, dbr+8} + vbc2bc.x[npost] +=1 + } else { + j -= 1 + if (debug_ == 1 ) {print "bc2bc"} + } + } + + // BC -> MC + for j=0, 2 { + Gauz3 = rdbc2mc.repick() // [-3:3] + if (i*2+2+Gauz3 > 14) {npost = i*2+2+Gauz3-15 } + if (i*2+2+Gauz3 < 0) {npost = i*2+2+Gauz3+15} + if ((i*2+2+Gauz3 >-1) && (i*2+2+Gauz3 < 15)) {npost = i*2+2+Gauz3} + // if ((is_connected(MossyCell[npost], BasketCell[i]) == 0) && (vbc2mc.x[npost] < 3) && (killMC.contains(ngcell+nbcell+npost) == 0)) { // use if killing MC + if (nbcell != 6) {max_conn = 4} else {max_conn = 3} + + c1 = cells.object(npost+ngcell+nbcell) + c2 = cells.object(i+ngcell) + if ((is_connected(c1, c2) == 0) && (vbc2mc.x[npost] < max_conn)) { // change < 3 to < 4 + nc_append(i+ngcell, npost+ngcell+nbcell, 12, 1.5e-3, 1.5, -10) // Gcell[3] to Bcell[1] + if (debug_ == 1 ) {print npost, 12} + vbc2mc.x[npost] +=1 + } else { + j -= 1 + if (debug_ == 1 ) {print "bc2mc"} + } + // if (killMC.contains(ngcell+nbcell+npost) == 1) {j +=1 if (debug_ == 1 ) {print "dead MC"}} // use if killing MC + } + } + } + + proc connect_mossy_cells() {local i,j localobj c1, c2 + /* Connect mossy cells to post-synaptic targets + */ + // Print to indicate step + print "Connecting MCs to post-synaptic targets." + + for i=0, nmcell-1 { + //if (killMC.contains(ngcell+nbcell+i) == 0) // use if killing MC + if (i < 3) { y=0 } + if ((i > 2) && (i < 6)) { y=1 } + if ((i > 5) && (i < 9)) { y=2 } + if ((i > 8) && (i < 12)) { y=3 } + if ((i > 11) && (i < 15)) { y=4 } + + // MC -> GC1 + for j=0, 99 { + Gauz1 = rdmc2gc1.repick() // [25:175] + if (i*33+17+Gauz1 > 499) { + npost1 = i*33+17+Gauz1-500 + } else { + npost1 =i*33+17+Gauz1 + } + + dbr = rdsyna.repick() // [0:1] + c1 = cells.object(npost1) + c2 = cells.object(i+ngcell+nbcell) + if ((is_connected(c1, c2) == 0) && (vmc2gc.x[npost1] < 7)) { + nc_append(i+ngcell+nbcell, npost1, dbr+2, 0.3e-3*scale_MC2GC_strength, 3, 10) + vmc2gc.x[npost1] +=1 + } else { + j -= 1 + } + } + + // MC -> GC2 + for j=0, 99 { + Gauz2 = rdmc2gc2.repick() // [-175:25] + if (i*33+17+Gauz2 < 0) { + npost2 =i*33+17+Gauz2+500 + } else {npost2 =i*33+17+Gauz2} + dbr = rdsyna.repick() // [0:1] + c1 = cells.object(npost2) + c2 = cells.object(i+ngcell+nbcell) + if ((is_connected(c1, c2) == 0) && (vmc2gc.x[npost2] < 7)) { + nc_append(i+ngcell+nbcell, npost2, dbr+2, 0.3e-3*scale_MC2GC_strength, 3, 10) // Gcell[3] to Bcell[1] + vmc2gc.x[npost2] +=1 + } else { + j -= 1 + } + } + + + // MC -> BC + for j=0, 0 { + Gauz3 = rdmc2bc.repick() // Gauz3 = [-3:3] + if (y+Gauz3 > nbcell-1) {npost = y+Gauz3-nbcell} // y = [0:4] => y+Gaus3 = [-3:7] + if (y+Gauz3 < 0) {npost = y+Gauz3+nbcell} + if ((y+Gauz3 > -1) && (y+Gauz3 < nbcell)) {npost = y+Gauz3} + dbr = rdsyna.repick() + if ((vmc2bc.x[npost] < 4) && (Gauz3 !=0)) { + nc_append(i+ngcell+nbcell, ngcell+npost, dbr+6, 0.3e-3, 3, 10) // Gcell[3] to Bcell[1] + vmc2bc.x[npost] += 1 + } else { + j -= 1 + } + } + + + // MC -> MC + for j=0, 2 { + Gauz3 = rdmc2mc1.repick() //[-3:3] + if (i+Gauz3 > 14) {npost = i+Gauz3-15 } + if (i+Gauz3 < 0) {npost = i+Gauz3+15} + if ((i+Gauz3 >-1) && (i+Gauz3 < 15)) {npost = i+Gauz3} + dbr = rdsynb.repick() + // if ((is_connected(MossyCell[npost], MossyCell[i]) == 0) && (vmc2mc.x[npost] < 4) && (Gauz3 != 0) && (killMC.contains(ngcell+nbcell+npost) == 0)) // use if killing MC + + c1 = cells.object(npost+ngcell+nbcell) + c2 = cells.object(i+ngcell+nbcell) + if ((is_connected(c1, c2) == 0) && (vmc2mc.x[npost] < 4) && (Gauz3 != 0)) { + nc_append(i+ngcell+nbcell, npost+ngcell+nbcell, dbr+8, 0.5e-3, 2, 10) + // print npost, dbr+8 + vmc2mc.x[npost] +=1 + } else { + j -= 1 + } + // if (killMC.contains(ngcell+nbcell+npost) == 1){ j += 1 print "dead MC"} // use if killing MC + } + + + // MC -> HC + for j=0, 1 { + Gauz3 = rdmc2hc.repick() // [-2:2] + if (y+Gauz3 > nhcell-1) {npost = y+Gauz3-nhcell} // changed code here to allow > 6 HCells + if (y+Gauz3 < 0) {npost = y+Gauz3+nhcell} + if ((y+Gauz3 > -1) && (y+Gauz3 < nhcell)) {npost = y+Gauz3} + dbr = rdsynb.repick() + + c1 = cells.object(ngcell+nbcell+nmcell+npost) + c2 = cells.object(i+ngcell+nbcell) + if ((is_connected(c1, c2) == 0) && (vmc2hc.x[npost] < 7) && (Gauz3 != 0)) { + nc_append(i+ngcell+nbcell, ngcell+nbcell+nmcell+npost, dbr+4, 0.2e-3, 3, 10) // Gcell[3] to Bcell[1] + // print npost, dbr+4 + vmc2hc.x[npost] +=1 + } else { + j -= 1 + } + } + } + } + + proc connect_hipp_cells() {local i,j localobj c1, c2 + /* Connect HIPP cells to post-synaptic targets + */ + // Print to indicate step + print "Connecting HIPP cells to post-synaptic targets." + + for i=0, nhcell-1 { + // HC -> GC + for j=0, 159 { // NOTE number of connections explicitly coded here! + Gauz3 = rdhc2gc.repick() // [-130:130] + //print Gauz3 + if (i*83+41+Gauz3 > 499) {npost = i*83+41+Gauz3-500 } // NOTE: 500 is explicitly coded here! + if (i*83+41+Gauz3 < 0) {npost = i*83+41+Gauz3+500} + if ((i*83+41+Gauz3 > -1) && (i*83+41+Gauz3 < 500)) {npost = i*83+41+Gauz3} + //print npost + dbr = rdsyna.repick() + if (nhcell != 6) {max_conn = 5} else {max_conn = 3} + + c1 = cells.object(npost) + c2 = cells.object(i+ngcell+nbcell+nmcell) + if ((is_connected(c1, c2) == 0) && (vhc2gc.x[npost] < max_conn)) { + nc_append(i+ngcell+nbcell+nmcell, npost, dbr+4, 0.5e-3*scale_HC2GC_strength, 1.6, 10) + vhc2gc.x[npost] +=1 + if (debug_ == 1 ) {print i, npost, dbr+4} + } else { + j -= 1 + if (debug_ == 1 ) {print "HC2GC"} + } + } + + // HC -> BC + for j=0, 3 { + Gauz3 = rdhc2bc.repick() // [-2:2] + if (i+Gauz3 > nbcell-1) {npost = i+Gauz3-nbcell} + if (i+Gauz3 < 0) {npost = i+Gauz3+nbcell} + if ((i+Gauz3 > -1) && (i+Gauz3 < nbcell)) {npost = i+Gauz3} + dbr = rdsyna.repick() + + c1 = cells.object(npost+ngcell) + c2 = cells.object(i+ngcell+nbcell+nmcell) + if ((is_connected(c1, c2) == 0) && (vhc2bc.x[npost] < 5)) { // NOTE < 5 coded explicitly here! + nc_append(i+ngcell+nbcell+nmcell, npost+ngcell, dbr+10, 0.5e-3, 1.6, 10) // Hcell to Bcell + if (debug_ == 1 ) {print npost, dbr+10} + vhc2bc.x[npost] += 1 + } else { + j -= 1 + if (debug_ == 1 ) {print "HC2BC"} + } + } + + // HC -> MC + for j=0, 3 { + Gauz3 = rdhc2mc.repick() // [-2:2] + //print Gauz3 + if (i*2+2+Gauz3 > 14) {npost = i*2+2+Gauz3-15 } + if (i*2+2+Gauz3 < 0) {npost = i*2+2+Gauz3+15} + if ((i*2+2+Gauz3 >-1) && (i*2+2+Gauz3 < 15)) {npost = i*2+2+Gauz3} + //print npost + dbr = rdsynb.repick() + // if ((is_connected(MossyCell[npost], HIPPCell[i]) == 0) && (vhc2mc.x[npost] < 2) && (killMC.contains(ngcell+nbcell+npost) == 0)) //use if killing MC + if (nhcell != 6) {max_conn = 4} else {max_conn = 2} + + c1 = cells.object(npost+ngcell+nbcell) + c2 = cells.object(i+ngcell+nbcell+nmcell) + if ((is_connected(c1, c2) == 0) && (vhc2mc.x[npost] < max_conn)) { // NOTE < 2 coded explicitly here!! => changed to 4 + nc_append(i+ngcell+nbcell+nmcell, npost+ngcell+nbcell, dbr+13, 1.5e-3, 1, 10) // Hcell to Mcell + //if (debug_ == 1 ) {print npost, dbr+13} + vhc2mc.x[npost] += 1 + } else { + j -= 1 + if (debug_ == 1 ) {print "HC2MC"} + } + // if (killMC.contains(ngcell+nbcell+npost) == 1){ j += 1 print "dead MC"} //use if killing MC + } + } + } + + // ######################################################################## + // FUNCTIONS FOR DATA STORAGE AND PRINTING + // ######################################################################## + proc update_voltage_trace_matrix() { local i + /* Procedure called in every time step to add voltage of + recorded cell at time t + */ + VmT.append(t) + for i=0, (ngcell+nbcell +nmcell +nhcell -1) { + VmMat[i].append( cells.object[i].soma.v(0.5)) + } + } + + proc convert_voltage_trace_to_spikes() { local i, j, max_id + /* Converts membrane potential vector to a spike train + */ + max_id = npp + if ( print_template == 1) { max_id = npp } + if ( print_Raster == 1) { max_id = ngcell+nbcell +nmcell +nhcell-1 } + + for i=0, (max_id) { + // Used to make a binary version of a spike train. + // Vector().spikebin(, ) + // is a vector of membrane potential. + // is the voltage threshold for spike detection. + // is set to all zeros except at the onset of spikes + // (the first dt which the spike crosses threshold) + Spike[i].spikebin(VmMat[i], 0) + } + } + + obfunc myprecell() { localobj nil + /* Read out the ID of pre or precell (by Ted Carnevale) + + Arguments: + $o1: NetCon + + */ + // nil points to NULLobject + if ($o1.precell() == nil) { + return $o1.pre() + } else { + return $o1.precell() + } + } + + proc save_network_connections(){ local i localobj dfile + /* Writes network adjacency list to file + */ + strdef DGNC_name_ + dfile = new File() + sprint (DGNC_name_, "%s%s-%d-%d-%s.%s", datadir, "DGNC", PP_box_startid_, PP_box_nr_, idname, suffix) + print "\tWriting network connections to file: ", DGNC_name_ + dfile.wopen(DGNC_name_) + for i=0, nclist.count-1 { + dfile.printf("%s\t%s\n", myprecell(nclist.object(i)), nclist.object(i).postcell) + } + dfile.close() + } + + proc initialize_voltage_trace_file(){ + /* Creates the header for the membrane voltage file + */ + strdef DGVt_name_ // Voltage file name + efile = new File() + + // [ TODO ] -- n_to_print here is not good. Needs to be an object parameter or something. + n_to_print = 500 + sprint (DGVt_name_,"%s%s-%d-%d-%s.%s", datadir, "DGVt", PP_box_startid_, PP_box_nr_, idname, suffix) + print "\tWriting membrane voltage traces of all cells to file: ", DGVt_name_ + efile.wopen(DGVt_name_) + efile.printf("# t\t") + efile.printf("\n") + efile.close(DGVt_name_) + } + + proc print_voltage_traces_to_file(){ local j + /* Prints voltage traces to file + + TODO: + - Have better way to implement n_to_print so that + we can print not only granule cells, but still do this neatly + */ + efile.aopen(DGVt_name_) + efile.printf("%f\t", t) + for i = 0, n_to_print-1 { + efile.printf("%f\t", cells.object[int(i*ngcell/n_to_print)].soma.v(0.5))} + efile.printf("\n") + efile.close(DGVt_name_) + } + + proc write_spike_times_to_file() { local i, j localobj spikefile + strdef DGsp_name_ + sprint (DGsp_name_,"%s%s-%d-%d-%s.%s", datadir, "DGsp", PP_box_startid_, PP_box_nr_, idname, suffix) + print "\tWriting spike times to file: ", DGsp_name_ + spikefile = new File() + spikefile.wopen(DGsp_name_) + k = 0 + while (k < VmT.size) { + for j = 0, (ngcell+nbcell+nmcell+nhcell-1) { + if(Spike[j].x[k] != 0) { + // Writes out time of spike and cell id! + spikefile.printf("%f\t%d\n", VmT.x[k], j) + } + } + k +=1 + } + spikefile.close(DGsp_name_) + } + + proc write_stimuli_to_file() {local i, j localobj stimfile, noisestimfile + /* Writes stimuli to file + */ + stimfile = new File() + strdef f_name_ + sprint (f_name_,"%s%s-%d-%d-%s.%s", datadir, "StimIn", PP_box_startid_, PP_box_nr_, idname, suffix) + + print "\tWriting input trains to file: ", f_name_ + + stimfile.wopen(f_name_) + for i=0,npp-1 { + for j=0,vec_stim[i].size()-1 { + stimfile.printf("%f\t%d\n",vec_stim[i].x[j],i) + } + } + stimfile.close(f_name_) + + // write out noise input... + if ((print_stim_noise == 1) && (spontaneous_activity_ == 1)) { + + noisestimfile = new File() + strdef f_name_noise_ + sprint (f_name_noise_,"%s%s-%d-%d-%s.%s", datadir, "NoiseStimIn", PP_box_startid_, PP_box_nr_, idname, suffix) + + noisestimfile.wopen(f_name_noise_) + for i=0,ngcell-1 { + for j=0,vec_stim_noise[i].size()-1 { + noisestimfile.printf("%f\t%d\n",vec_stim_noise[i].x[j],i) + } + } + noisestimfile.close(f_name_noise_) + } + + + } + + // ######################################################################## + // FUNCTIONS TO INITIALIZE AND RUN SIMULATION + // ######################################################################## + proc initialize_simulation_run() { local dtsav, temp, secsav + /* Initialize a simulation run + + Arguments: + - $1 : int : 100 [ TODO ] NEED TO UNDERSTAND THIS MORE. SOMETHING ABOUT SEED SETTING + + */ + + finitialize(v_init) + // [ TODO ] - Figure out what these three next lines are for + t = -1000 // negative t step to initialize to steady state + dtsav = dt // original simulation step size + dt= 10 // larger step size + + // if cvode is on, turn it off to do large fixed step + // [ TODO ] - Learn why + temp = cvode.active() + if ( temp != 0 ) { cvode.active(0) } + while( t < -100 ) { fadvance() } + + //restore cvode if reqd + if (temp!=0) { cvode.active(1) } + dt = dtsav + t = 0 + + if ( cvode.active() ){ + cvode.re_init() + } else { + fcurrent() + } + + // restart number generators for PPStim + t = 0 + finitialize(v_init) + trial_old = trial + + if ($1 == 100) { + // Reset all Poisson inputs with the correct seed + for i = 0, rslist.count()-1 rslist.o(i).start(trial) + } else { + // reset int((rslist.count()-1)*(1-$1/100.)) Poisson inputs with a different seed + trial = trial + 1 + for i = 0, int((rslist.count()-1)*(1-$1/100.)) rslist.o(i).start(trial) + + // reset rslist.count()-int((rslist.count()-1)*(1-$1/100.)) Poisson inputs with the correct seed + trial = trial_old + for i = int((rslist.count()-1)*(1-$1/100.)), rslist.count()-1 rslist.o(i).start(trial) + } + + // init noise Generators (with seed trial_noise ...) + // reset all Poisson noise inputs with the trial_noise seed + trial = trial_noise + for i = 0, rs_noiselist.count()-1 rs_noiselist.o(i).start(trial) + trial = trial_old + + VmT = new Vector() + for i=0, cells.count-1 { + VmMat[i] = new Vector() + } + + for i=0, (ngcell + nbcell + nmcell + nhcell -1) { + Spike[i] = new Vector() // vector of spikes for each cell + Spike_times[i] = new Vector() // vector of spikes for each cell + } + + // finalize at the end + // [ TODO ] - Figure out what's going on here + t = 0 + finitialize(v_init) + frecord_init() + + // Select input pattern + for i=PP_box_startid_, PP_box_startid_+PP_box_nr_-1 { + PPSt[i].pp.status = 1 + PPSt[i].pp.start = PP_box_start_ + PPSt[i].pp.forcestop = PP_box_stop_ + PPSt[i].pp.nspk = PP_box_nspk_ + } + + // Initialize the GC spontaneous activity generators + // [ TODO ] Need more general way to set the number that receive SA generators + // This needs to be based on the desired frequency of spont. activity. + if (spontaneous_activity_ == 1) { + // Rdsa_thresh is the number of granule cells that should be active + // with a single spike in the simulation window in order to achieve + // a spontaneous firing rate of `spontaneous_activity_rate` + rdsa_thresh = ngcell*(spontaneous_activity_rate/(1000/tstop)) + for i=0,ngcell-1 { + rdsa_rv = rdsa.repick() + if (rdsa_rv <= rdsa_thresh) { + PPSt_noise[i].pp.status = 1 + PPSt_noise[i].pp.start = 0. + PPSt_noise[i].pp.forcestop = tstop + PPSt_noise[i].pp.nspk = 1 + } + if (rdsa_rv > rdsa_thresh) { + PPSt_noise[i].pp.status = 0 + PPSt_noise[i].pp.start = 0. + PPSt_noise[i].pp.forcestop = tstop + PPSt_noise[i].pp.nspk = 0 + } + } + } + } + + proc run() {local i, j + // MAKE NETWORK + make_cells() + make_input_stimulation() + make_connections() + + // INITIALIZE OPTIMIZER + cvode = new CVode() + + // RE-INITIALIZE DATA CONTAINERS + objref Spike[ntotal-1], Spike_times[ntotal-1] + objref VmT, VmMat[cells.count] + + // RUN SIMULATION + print "Running simulation." + + // Write connections to file + save_network_connections() + for (PP_box_startid_ = 0; PP_box_startid_ <= n_patterns; PP_box_startid_ += 1 ) { + initialize_simulation_run(100) + + // File header voltage output file + if (print_Vtrace == 1) { + initialize_voltage_trace_file() + } + + while (t 1) { scale_na_conductances = $2 } + if (narg > 2) { scale_kdr_conductances = $3 } + if (narg > 3) { scale_ka_conductances = $4 } + if (narg > 4) { gbar_ht_ = $5 } + if (narg > 5) { gbar_lt_ = $6 } + if (narg > 6) { scale_size_ = $7 } + if (narg > 7) { scale_gpas_dg_ = $8 } + if (narg > 8) { scale_sk_dg_ = $9 } + if (narg > 9) { scale_gabaa_ = $10 } + if (narg > 10) { scale_kir_ = $11 } + + // Run actual initialization + pre_list = new List() + subsets() + gctemp() + synapse() + } + + proc subsets(){ local i + all = new SectionList() + soma all.append() + for i=0, 3 gcdend1 [i] all.append() + for i=0, 3 gcdend2 [i] all.append() + + gcldend = new SectionList() + gcdend1 [0] gcldend.append() + gcdend2 [0] gcldend.append() + + pdend = new SectionList() + gcdend1 [1] pdend.append() + gcdend2 [1] pdend.append() + + mdend = new SectionList() + gcdend1 [2] mdend.append() + gcdend2 [2] mdend.append() + + ddend = new SectionList() + gcdend1 [3] ddend.append() + gcdend2 [3] ddend.append() + } + proc gctemp() { + + scale_area = 1./1.13 * scale_size_ + + // ********** Parameters for reversal potentials (assigned below) ********* + e_gabaa_ = -70. // reversal potential GABAA + + // ***************** Parameters + g_pas_fit_ = 1.44e-05 * scale_gpas_dg_ + gkbar_kir_fit_ = 1.44e-05 * scale_kir_ + ggabaabar_fit_ = 0.722e-05 * scale_gabaa_ + + // *********************** PAS ****************************************** + cm_fit_ = 1. + Ra_fit_ = 184. // fitted + + // *********************** KIR ***************************************** + vhalfl_kir_fit_ = -98.923594 // for Botzman I/V curve, fitted + kl_kir_fit_ = 10.888538 // for Botzman I/V curve, fitted + q10_kir_fit_ = 1. // temperature factor, set to 1 + vhalft_kir_fit_ = 67.0828 // 3 values for tau func from Stegen et al. 2011 + at_kir_fit_ = 0.00610779 + bt_kir_fit_ = 0.0817741 + + // ********************* Neuron Morphology etc *************************** + LJP_ = -10. // Liquid junction potential [mV] + V_rest = -68.16+LJP_ // resting potential [mV] + V_init = -68.16+LJP_ // initial potential [mV] + + // ******************** GABAA ******************** + e_pas_fit_ = -83.8 + e_pas_fit_Dend = -81.74 + + soma {nseg=1 L=16.8*scale_area diam=16.8*scale_area} // changed L & diam + + gcdend1 [0] {nseg=1 L=50*scale_area diam=3*scale_area} + for i = 1, 3 gcdend1 [i] {nseg=1 L=150*scale_area diam=3*scale_area} + + gcdend2 [0] {nseg=1 L=50*scale_area diam=3*scale_area} + for i = 1, 3 gcdend2 [i] {nseg=1 L=150*scale_area diam=3*scale_area} + + forsec all { + insert ccanl + catau_ccanl = 10 + caiinf_ccanl = 5.e-6 + insert HT + gbar_HT = gbar_ht_ + kan_HT = 0.5 + kbn_HT = 0.3 + insert LT + gbar_LT = gbar_lt_ + Ra=Ra_fit_ + } + + soma {insert ichan2 + gnatbar_ichan2=0.12 * scale_na_conductances // value Aradi & Holmes 1999 + gkfbar_ichan2=0.016 * scale_kdr_conductances + gksbar_ichan2=0.006 * scale_kdr_conductances + gl_ichan2 = g_pas_fit_ + el_ichan2 = e_pas_fit_ // set leak reversal poti to gain Vrest of cell + insert ka + gkabar_ka=0.012 * scale_ka_conductances + insert km + gbar_km = 0.001 * scale_kdr_conductances + insert nca + gncabar_nca=0.002 + insert lca + glcabar_lca=0.005 + insert tca + gcatbar_tca=0.000037 + insert sk + gskbar_sk=0.001 * scale_sk_dg_ + insert bk + gkbar_bk=0.0006 + cm=cm_fit_ + } + + forsec gcldend {insert ichan2 + gnatbar_ichan2=0.018 * scale_na_conductances // value Aradi & Holmes 1999 + gkfbar_ichan2=0.004 * scale_kdr_conductances + gksbar_ichan2=0.006 * scale_kdr_conductances + gl_ichan2 = g_pas_fit_ + el_ichan2 = e_pas_fit_ // set leak reversal poti to gain Vrest of cell + insert nca // HAV-N- Ca channel + gncabar_nca=0.003 // value Aradi & Holmes 1999 + insert lca + glcabar_lca=0.0075 + insert tca + gcatbar_tca=0.000075 + insert sk + gskbar_sk=0.0004 * scale_sk_dg_ + insert bk + gkbar_bk=0.0006 + cm=cm_fit_ + } + + forsec pdend {insert ichan2 + gnatbar_ichan2=0.013 * scale_na_conductances // value Aradi & Holmes 1999 + gkfbar_ichan2=0.004 * scale_kdr_conductances + gksbar_ichan2=0.006 * scale_kdr_conductances + gl_ichan2 = g_pas_fit_ * (0.000063/0.00004) + el_ichan2 = e_pas_fit_Dend // see comment above + insert nca // HAV-N- Ca channel + gncabar_nca=0.001 // value Aradi & Holmes 1999 + insert lca + glcabar_lca=0.0075 + insert tca + gcatbar_tca=0.00025 + insert sk + gskbar_sk=0.0002 * scale_sk_dg_ + insert bk + gkbar_bk=0.001 + cm=cm_fit_*1.6 + } + + forsec mdend {insert ichan2 + gnatbar_ichan2=0.008 * scale_na_conductances // value Aradi & Holmes 1999 + gkfbar_ichan2=0.001 * scale_kdr_conductances + gksbar_ichan2=0.006 * scale_kdr_conductances + gl_ichan2 = g_pas_fit_ * (0.000063/0.00004) + el_ichan2 = e_pas_fit_Dend // see comment above + insert nca + gncabar_nca=0.001 // value Aradi & Holmes 1999 + insert lca + glcabar_lca=0.0005 + insert tca + gcatbar_tca=0.0005 + insert sk + gskbar_sk=0.0 * scale_sk_dg_ + insert bk + gkbar_bk=0.0024 + cm=cm_fit_*1.6 + } + + forsec ddend {insert ichan2 + gnatbar_ichan2=0.0 * scale_na_conductances + gkfbar_ichan2=0.001 * scale_kdr_conductances + gksbar_ichan2=0.008 * scale_kdr_conductances + gl_ichan2 = g_pas_fit_ * (0.000063/0.00004) + el_ichan2 = e_pas_fit_Dend // see comment above + insert nca + gncabar_nca=0.001 // value Aradi & Holmes 1999 + insert lca + glcabar_lca=0.0 + insert tca + gcatbar_tca=0.001 + insert sk + gskbar_sk=0.0 * scale_sk_dg_ + insert bk + gkbar_bk=0.0024 + cm=cm_fit_*1.6 + } + + + connect gcdend1[0](0), soma(1) + connect gcdend2[0](0), soma(1) + for i=1,3 { + connect gcdend1[i](0), gcdend1[i-1](1) + } + for i=1,3 { + connect gcdend2[i](0), gcdend2[i-1](1) + } + + forsec all { + insert kir // kir conductance added in Yim et al. 2015, note that eK=-90mV is used instead of -105mV as reported in the paper + gkbar_kir = gkbar_kir_fit_ + vhalfl_kir = vhalfl_kir_fit_ + kl_kir = kl_kir_fit_ + vhalft_kir = vhalft_kir_fit_ + at_kir = at_kir_fit_ + bt_kir = bt_kir_fit_ + ggabaa_ichan2 = ggabaabar_fit_ // added GabaA in Yim et al. 2015 + egabaa_ichan2 = e_gabaa_ // reversal potential GABAA added in Yim et al. 2015 + ena = 50 // ena was unified from enat=55 (BC, HIPP, MC) and enat=45 (GC) in Santhakumar et al. (2005) + ek = -90 // simplified ekf=eks=ek=esk; note the eK was erroneously reported as -105mV in the Yim et al. 2015 + cao_ccanl = 2 } +} // end of gctemp() + + // Retrieval of objref arguments uses the syntax: $o1, $o2, ..., $oi. + // http://web.mit.edu/neuron_v7.1/doc/help/neuron/general/ocsyntax.html#arguments + proc connect_pre() { + soma $o2 = new NetCon (&v(1), $o1) + } + + + // Define synapses on to GCs using + //- an Exp2Syn object (parameters tau1 -rise, tau2 -decay, + // time constant [ms] and e - rev potential [mV] + // delay [ms] and weight -variablr betw 0 and 1 [1 corresponding to 1 'S] + + proc synapse() { + gcdend1[3] syn = new Exp2Syn(0.5) // PP syn based on data from Greg Hollrigel and Kevin Staley NOTE: both synapses are identical! + syn.tau1 = 1.5 syn.tau2 = 5.5 syn.e = 0 + pre_list.append(syn) + + gcdend2[3] syn = new Exp2Syn(0.5) // PP syn based on Greg and Staley + syn.tau1 = 1.5 syn.tau2 = 5.5 syn.e = 0 + pre_list.append(syn) + + gcdend1[1] syn = new Exp2Syn(0.5) // MC syn *** Estimated + syn.tau1 = 1.5 syn.tau2 = 5.5 syn.e = 0 + pre_list.append(syn) + + gcdend2[1] syn = new Exp2Syn(0.5) // MC syn *** Estimated + syn.tau1 = 1.5 syn.tau2 = 5.5 syn.e = 0 + pre_list.append(syn) + + gcdend1[3] syn = new Exp2Syn(0.5) // HIPP syn based on Harney and Jones corrected for temp + syn.tau1 = 0.5 syn.tau2 = 6 syn.e = -70 + pre_list.append(syn) + + gcdend2[3] syn = new Exp2Syn(0.5) // HIPP syn based on Harney and Jones corrected for temp + syn.tau1 = 0.5 syn.tau2 = 6 syn.e = -70 + pre_list.append(syn) + + soma syn = new Exp2Syn(0.5) // BC syn based on Bartos + syn.tau1 = 0.26 syn.tau2 = 5.5 syn.e = -70 + pre_list.append(syn) + + gcdend1[1] syn = new Exp2Syn(0.5) // NOTE: SPROUTED SYNAPSE based on Molnar and Nadler + syn.tau1 = 1.5 syn.tau2 = 5.5 syn.e = 0 + pre_list.append(syn) + + gcdend2[1] syn = new Exp2Syn(0.5) // NOTE: SPROUTED SYNAPSE + syn.tau1 = 1.5 syn.tau2 = 5.5 syn.e = 0 + pre_list.append(syn) + + // Total of 7 synapses per GC 0,1 PP; 2,3 MC; 4,5 HIPP and 6 BC 7,8 Sprout + } + + func is_art() { return 0 } + +endtemplate GranuleCell diff --git a/HIPP.hoc b/objects/HIPP.hoc similarity index 99% rename from HIPP.hoc rename to objects/HIPP.hoc index 65a9268..6ae1dbe 100644 --- a/HIPP.hoc +++ b/objects/HIPP.hoc @@ -11,11 +11,10 @@ // http://onlinelibrary.wiley.com/doi/10.1002/hipo.22373/abstract // modified by +// Abraham Nunes / 2022 // Man Yi Yim / 2015 // Alexander Hanuschkin / 2011 -objref Hcell[nhcell] - begintemplate HIPPCell ndend1=3 diff --git a/MC.hoc b/objects/MC.hoc similarity index 99% rename from MC.hoc rename to objects/MC.hoc index 389c94f..6c67954 100644 --- a/MC.hoc +++ b/objects/MC.hoc @@ -12,11 +12,10 @@ // http://onlinelibrary.wiley.com/doi/10.1002/hipo.22373/abstract // modified by +// Abraham Nunes / 2022 // Man Yi Yim / 2015 // Alexander Hanuschkin / 2011 -objref Mcell[nmcell] - begintemplate MossyCell ndend1=4 diff --git a/PP.hoc b/objects/PP.hoc similarity index 78% rename from PP.hoc rename to objects/PP.hoc index 16ed86a..0ede0c1 100644 --- a/PP.hoc +++ b/objects/PP.hoc @@ -17,8 +17,8 @@ begintemplate PPstim -external PP_rate_, tstop -external PP_box_start_,PP_box_stop_ // stimulation window instead of Poisson input ... +//external PP_rate_, tstop +//external PP_box_start_,PP_box_stop_ // stimulation window instead of Poisson input ... public pp, connect_pre, is_art, acell create acell @@ -26,7 +26,12 @@ objref pp proc init() { - actemp() + actemp() + pp_index = $1 + PP_rate_ = $2 // Is this even used?? + tstop = $3 + PP_box_start_ = $4 + PP_box_stop_ = $5 } proc actemp() { diff --git a/objects/ranstream.hoc b/objects/ranstream.hoc new file mode 100644 index 0000000..086f7c6 --- /dev/null +++ b/objects/ranstream.hoc @@ -0,0 +1,31 @@ +// Modified by Abraham Nunes / 2022 +// Removed dependency on an external global variables to be +// better contained. +begintemplate RandomStream + public r, repick, start, stream + objref r + proc init() { + /* Initializes random stream + + Arguments: + $1 : int : Stream ID + $2 : int : Seed + */ + stream = $1 + random_stream_offset_ = 1000 + r = new Random() + start($2) + } + func start() { + /* Resets the RNG + + Arguments: + $1 : int : Seed + + */ + return r.MCellRan4(stream*random_stream_offset_ + 1 + $1) + } + func repick() { + return r.repick() + } +endtemplate RandomStream diff --git a/plot-fi-curves.jl b/plot-fi-curves.jl new file mode 100644 index 0000000..9bd602e --- /dev/null +++ b/plot-fi-curves.jl @@ -0,0 +1,29 @@ +using Plots, DataFrames, CSV, LaTeXStrings + +df = CSV.read("ficurve.txt", delim="\t", header=0, DataFrame) +k_unique = unique(df[:,4]) + +labels = ["HC", "LR", "NR"] +ls = [:solid, :dash, :dot] +p = plot(xlabel=L"pA", ylabel=L"Hz", legend=:topleft) +for i ∈ 1:3 + p = plot!(1000 .* df[df[:,4] .== k_unique[i],1], df[df[:,4] .== k_unique[i], 2], + c=:black, ls=ls[i], label=labels[i]) +end +p + +V = DataFrame() +for i ∈ 0:2 + v0 = CSV.read("vm-$i.000000.txt", delim="\t", header=0, DataFrame) + v0[:,"ID"] .= i + V = [V; v0] +end + +p = plot() +for i ∈ 0:2 + p = plot!(V[V.ID .== i,1], V[V.ID .== i, 10], label="$i") +end +p + + +4.5 * (16.8^2 * π \ No newline at end of file diff --git a/plot_DG_all.py b/plot_DG_all.py deleted file mode 100644 index 108d949..0000000 --- a/plot_DG_all.py +++ /dev/null @@ -1,401 +0,0 @@ -### Analysis of DG network data ### -# This Python code plots DG neurons' activity. -# Enter the idname, sid and nid. - -# ModelDB file along with publication: -# Yim MY, Hanuschkin A, Wolfart J (2015) Hippocampus 25:297-308. -# http://onlinelibrary.wiley.com/doi/10.1002/hipo.22373/abstract - -# modified and augmented by -# Man Yi Yim / 2015 -# Alexander Hanuschkin / 2011 - -import pylab -import numpy -import matplotlib as mpl - -# Font size in the figure -font_size = 12 # 20 -mpl.rcParams['axes.titlesize'] = font_size+10 -mpl.rcParams['xtick.labelsize'] = font_size+6 -mpl.rcParams['ytick.labelsize'] = font_size+6 -mpl.rcParams['axes.labelsize'] = font_size+8 -mpl.rcParams['legend.fontsize'] = font_size -mpl.rcParams['font.size'] = font_size+10 - -# Parameters -tstop = 200 # duration of the simulation to be analysed -is_all_Vt = 1 # 1 if the dgvtfile contains all GCs -selected_Vt = numpy.array([135,157,185,224,239]) # used if is_all_Vt is 0; put 5 IDs of GCs to print -npp = 100 # number of PP inputs -sid = 2 # PP_box_start_ -nid = 6 # number of active PPs -#idname = "-pp10-gaba1-kir1-st0" -print 'duration = ' + str(tstop) -print 'npp = ' + str(npp) -print 'idname = ' + idname - -# Import data -#bginfile = 'BGIn-'+str(sid)+'-'+str(nid)+idname+'.txt' -stiminfile = 'StimIn-'+str(sid)+'-'+str(nid)+idname+'.txt' -dgspfile = 'DGsp-'+str(sid)+'-'+str(nid)+idname+'.txt' -dgvtfile = 'DGVt-'+str(sid)+'-'+str(nid)+idname+'.txt' -con_file = 'DGNC-'+str(sid)+'-'+str(nid)+idname+'.txt' -con_file_id = 'id_DGNC.txt' - -# Load files -#f1 = open(bginfile,'r') -#if f1.readline() == '': -# bgin = numpy.empty(0) -#else: -# bgin = numpy.loadtxt(bginfile) -f2 = open(stiminfile,'r') -if f2.readline() == '': - stimin = numpy.empty(0) -else: - stimin = numpy.loadtxt(stiminfile) -f3 = open(dgspfile,'r') -if f3.readline() == '': - dgsp = numpy.empty(0) -else: - dgsp = numpy.loadtxt(dgspfile) -dgvt = numpy.loadtxt(dgvtfile) - -def extract_conn(): - """extract_conn() extracts connectivity and prints the global IDs for Santhakumar model by default. - Modify the number in this code if you change the network size or ordering. - ID 0-499: GC - ID 500-505: BC - ID 506-520: MC - ID 521-526: HC - ID 527-626: PP - ID >=527: BG """ - f_read = open(con_file) - f_write = open(con_file_id,'w') - for line in f_read: - if line[0] == 'G': - j = 13 - while line[j] != ']': - j += 1 - pre = int(line[12:j]) - elif line[0] == 'B': - j = 12 - while line[j] != ']': - j += 1 - pre = 500 + int(line[11:j]) - elif line[0] == 'M': - j = 11 - while line[j] != ']': - j += 1 - pre = 506 + int(line[10:j]) - elif line[0] == 'H': - j = 10 - while line[j] != ']': - j += 1 - pre = 521 + int(line[9:j]) - elif line[0] == 'N': - j = 9 #12 - while line[j] != ']': - j += 1 - if line[7] == 'B': # correlated PP - pre = 527 + int(line[11:j]) - elif line[7] == '1': # uncorrelated PP - pre = 527 + int(line[11:j]) - elif line[7] == '[': # BG - pre = 527 + npp + int(line[8:j]) - else: - print 'Please check the file: ' + con_file - j += 2 - if line[j] == 'G': - j += 13 - k = j - while line[j] != ']': - j += 1 - post = int(line[k-1:j]) - elif line[j] == 'B': - j += 12 - k = j - while line[j] != ']': - j += 1 - post = 500 + int(line[k-1:j]) - elif line[j] == 'M': - j += 11 - k = j - while line[j] != ']': - j += 1 - post = 506 + int(line[k-1:j]) - elif line[j] == 'H': - j += 10 - k = j - while line[j] != ']': - j += 1 - post = 521 + int(line[k-1:j]) - else: - print 'Please check the file: ' + con_file - f_write.write(str(pre)+'\t'+str(post)+'\n') - f_write.close() - -def conn_mat(is_weighted = 1): - """conn_mat(is_weighted = 1) prints out the connectivity matrix according to con_file_id - if is_weighted = 0, the matrix element is either connected (white) or unconnected (black) - if is_weighted = 1, the matrix element is the connection strength that has to be entered manually. - """ - conn = numpy.loadtxt(con_file_id) - ncell = 527 - npp = 100 - # connection strength - gc2gc = 0 - gc2bc = 4.7e-3 - gc2mc = 0.2e-3 - gc2hc = 0.5e-3 - bc2gc = 1.6e-3 - bc2bc = 7.6e-3 - bc2mc = 1.5e-3 - bc2hc = 0 - mc2gc = 0.3e-3 - mc2bc = 0.3e-3 - mc2mc = 0.5e-3 - mc2hc = 0.2e-3 - hc2gc = 0.5e-3 - hc2bc = 0.5e-3 - hc2mc = 1.5e-3 - hc2hc = 0 - pp2gc = 0.02 - pp2bc = 0 - if is_weighted: - cmat = numpy.zeros((ncell+npp,ncell)) - for j in range(conn.size/2): - if conn[j,0] < 500: - if conn[j,1] < 500: - cmat[ncell+npp-conn[j,0]-1,conn[j,1]] = gc2gc - elif conn[j,1] < 506: - cmat[ncell+npp-conn[j,0]-1,conn[j,1]] = gc2bc - elif conn[j,1] < 521: - cmat[ncell+npp-conn[j,0]-1,conn[j,1]] = gc2mc - elif conn[j,1] < 527: - cmat[ncell+npp-conn[j,0]-1,conn[j,1]] = gc2hc - elif conn[j,0] < 506: - if conn[j,1] < 500: - cmat[ncell+npp-conn[j,0]-1,conn[j,1]] = bc2gc - elif conn[j,1] < 506: - cmat[ncell+npp-conn[j,0]-1,conn[j,1]] = bc2bc - elif conn[j,1] < 521: - cmat[ncell+npp-conn[j,0]-1,conn[j,1]] = bc2mc - elif conn[j,1] < 527: - cmat[ncell+npp-conn[j,0]-1,conn[j,1]] = bc2hc - elif conn[j,0] < 521: - if conn[j,1] < 500: - cmat[ncell+npp-conn[j,0]-1,conn[j,1]] = mc2gc - elif conn[j,1] < 506: - cmat[ncell+npp-conn[j,0]-1,conn[j,1]] = mc2bc - elif conn[j,1] < 521: - cmat[ncell+npp-conn[j,0]-1,conn[j,1]] = mc2mc - elif conn[j,1] < 527: - cmat[ncell+npp-conn[j,0]-1,conn[j,1]] = mc2hc - elif conn[j,0] < 527: - if conn[j,1] < 500: - cmat[ncell+npp-conn[j,0]-1,conn[j,1]] = hc2gc - elif conn[j,1] < 506: - cmat[ncell+npp-conn[j,0]-1,conn[j,1]] = hc2bc - elif conn[j,1] < 521: - cmat[ncell+npp-conn[j,0]-1,conn[j,1]] = hc2mc - elif conn[j,1] < 527: - cmat[ncell+npp-conn[j,0]-1,conn[j,1]] = hc2hc - elif conn[j,0] < 527 + npp: - if conn[j,1] < 500: - cmat[ncell+npp-conn[j,0]-1,conn[j,1]] = pp2gc - elif conn[j,1] < 506: - cmat[ncell+npp-conn[j,0]-1,conn[j,1]] = pp2bc - # The first figure is the connectivity matrix of projections from HC, MC and BC to all - pylab.figure() - pylab.pcolor(cmat[npp:ncell+npp-500,:],cmap=pylab.cm.spectral) - pylab.colorbar() - pylab.title('Connectivity matrix') - pylab.axis([0,ncell,0,27]) - pylab.xlabel('Post') - pylab.ylabel('Pre') - pylab.yticks([3,13.5,24],('HC','MC','BC')) - pylab.xticks([250,503,513.5,523],('GC','BC','MC','HC')) - pylab.savefig('120412_conn_mat_to_all.png') - # The second figure is the connectivity matrix of projections from all to HC, MC and BC - pylab.figure() - pylab.pcolor(cmat[npp:,500:],cmap=pylab.cm.spectral) - pylab.colorbar() - pylab.title('Connectivity matrix') - pylab.axis([0,27,0,ncell]) - pylab.xlabel('Post') - pylab.ylabel('Pre') - pylab.yticks([3,13.5,24,177],('HC','MC','BC','GC')) - pylab.xticks([3,13.5,23],('BC','MC','HC')) - pylab.savefig('120412_conn_mat_from_all.png') - else: - cmat = numpy.ones((ncell+npp,ncell)) - for j in range(conn.size/2): - if conn[j,0] < ncell+npp: # background input is not considered - cmat[ncell+npp-conn[j,0]-1,conn[j,1]] = 0 - # This figure is the connectivity matrix of projections from PP and all neurons to all neurons - pylab.figure(figsize=(16,12)) - pylab.pcolor(cmat,cmap=pylab.cm.gray) - pylab.title('Connectivity matrix: white means connected') - pylab.axis([0,ncell,0,ncell+npp]) - pylab.xlabel('Post') - pylab.ylabel('Pre') - pylab.yticks([50,103,113.5,124,377],('PP','HC','MC','BC','GC')) - pylab.xticks([250,503,513.5,523],('GC','B','M','H')) - pylab.savefig('120412_conn_mat.png') - -# Load connectivity -extract_conn() # extract the connectivity profile 'con_file_id' for further analysis -# conn_mat(1) # print out the connectivity matrix figure -conn = numpy.loadtxt(con_file_id) -pylab.figure(figsize=(10,15)) - -# Plot PP activity -pylab.subplot(411) -pylab.plot([-1],[-1],linewidth=2) -pylab.plot([-1],[-1],linewidth=2) -pylab.plot([-1],[-1],linewidth=2) -pylab.plot([-1],[-1],linewidth=2) -pylab.plot([-1],[-1],'k',linewidth=2) -pylab.legend(['GC','BC','MC','HC','PP']) -for j in range(stimin.size/2): - if stimin.size == 2: - pylab.plot(stimin[0]*numpy.ones(2),numpy.array([stimin[1]-0.4,stimin[1]+0.4]),'k',linewidth=2) - else: - pylab.plot(stimin[j,0]*numpy.ones(2),numpy.array([stimin[j,1]-0.4,stimin[j,1]+0.4]),'k',linewidth=2) -#pylab.xticks(range(0,201,25)) -pylab.axis([0,tstop,-0.5,20-0.5]) -pylab.ylabel('PP') -pylab.title(idname) - -# Plot DG activity -pylab.subplot(412) -if dgsp.size > 0: - for j in range(dgsp.size/2): - if dgsp.size == 2: - if dgsp[1] < 500: - color = 'b' - elif dgsp[1] < 506: - color = 'g' - elif dgsp[1] < 521: - color = 'r' - elif dgsp[1] < 527: - color = 'c' - pylab.plot(dgsp[0]*numpy.ones(2),numpy.array([dgsp[1]-0.4,dgsp[1]+0.4]),color,linewidth=2) - else: - if dgsp[j,1] < 500: - color = 'b' - elif dgsp[j,1] < 506: - color = 'g' - elif dgsp[j,1] < 521: - color = 'r' - elif dgsp[j,1] < 527: - color = 'c' - pylab.plot(dgsp[j,0]*numpy.ones(2),numpy.array([dgsp[j,1]-0.4,dgsp[j,1]+0.4]),color,linewidth=2) -#pylab.xticks(range(0,201,25)) -pylab.axis([0,tstop,-0.5,527-0.5]) -pylab.ylabel('DG') - -# Plot spiking GCs and their inputs -nspk = 0 -if dgsp.size == 2 and dgsp[1] < 500: - nspk = 1 -elif dgsp.size > 2: - nspk = numpy.unique(dgsp[dgsp[:,1]<500,1]).size -if is_all_Vt: - spkgc = numpy.arange(5) - if dgsp.size > 0: - if dgsp.size == 2 and dgsp[1] < 500: - spkgc = dgsp[1] - elif dgsp.size > 2: - spkgc = numpy.unique(dgsp[dgsp[:,1]<500,1]) - if spkgc.size > 5: - spkgc = spkgc[:5] - elif spkgc.size < 5: - for j in range(5-spkgc.size): - spkgc = numpy.append(spkgc,[j]) # append a GC -else: - spkgc = selected_Vt -spkgc.sort() -print 'selected neurons = ' + str(spkgc) - -# Rasterplot of a few GCs and their inputs -pylab.subplot(413) -for j in range(dgsp.size/2): - if dgsp.size == 2: - if any(spkgc == dgsp[1]): - n = pylab.argwhere(spkgc==dgsp[1])[0] - pylab.plot(dgsp[0]*numpy.ones(2),numpy.array([n*2-0.4,n*2+0.4])+1,'b',linewidth=2) - else: - if any(spkgc == dgsp[j,1]): - n = pylab.argwhere(spkgc==dgsp[j,1])[0] - pylab.plot(dgsp[j,0]*numpy.ones(2),numpy.array([n*2-0.4,n*2+0.4])+1,'b',linewidth=2) -for j in range(spkgc.size): - ind = conn[conn[:,1]==spkgc[j],0] - for k in ind: - if k < 527: - for l in range(dgsp.size/2): - if dgsp.size == 2: - if dgsp[1] == k: - if k < 500: - color = 'b' - elif k < 506: - color = 'g' - elif k < 521: - color = 'r' - elif k < 527: - color = 'c' - pylab.plot(dgsp[0]*numpy.ones(2),numpy.array([j*2-0.4,j*2+0.4])+2,color,linewidth=2) - else: - if dgsp[l,1] == k: - if k < 500: - color = 'b' - elif k < 506: - color = 'g' - elif k < 521: - color = 'r' - elif k < 527: - color = 'c' - pylab.plot(dgsp[l,0]*numpy.ones(2),numpy.array([j*2-0.4,j*2+0.4])+2,color,linewidth=2) - else: - for l in range(stimin.size/2): - if stimin.size == 2: - if stimin[1] == k-527: - pylab.plot(stimin[0]*numpy.ones(2),numpy.array([j*2-0.4,j*2+0.4])+2,'k',linewidth=2) - else: - if stimin[l,1] == k-527: - pylab.plot(stimin[l,0]*numpy.ones(2),numpy.array([j*2-0.4,j*2+0.4])+2,'k',linewidth=2) -# for l in range(bgin.size/2): -# if bgin.size == 2: -# if bgin[1] == k-527-npp: -# pylab.plot(bgin[0]*numpy.ones(2),numpy.array([j*2-0.4,j*2+0.4])+2,'0.5',linewidth=2) -# else: -# if bgin[l,1] == k-527-npp: -# pylab.plot(bgin[l,0]*numpy.ones(2),numpy.array([j*2-0.4,j*2+0.4])+2,'0.5',linewidth=2) -pylab.plot([0,tstop],[2.5,2.5],'k') -pylab.plot([0,tstop],[4.5,4.5],'k') -pylab.plot([0,tstop],[6.5,6.5],'k') -pylab.plot([0,tstop],[8.5,8.5],'k') -pylab.yticks(numpy.arange(1.5,10,2),spkgc.astype(int)) -#pylab.xticks(range(0,201,25)) -pylab.axis([0,tstop,0.5,spkgc.size*2+0.5]) -pylab.ylabel('GCs & input') - -# Plot membrane potentials -pylab.subplot(414) -for j in range(spkgc.size): -# if spkgc[j] > 5: - pylab.plot(dgvt[:,0],dgvt[:,spkgc[j]+1],'b',linewidth=2) -# else: -# pylab.plot(dgvt[:,0],dgvt[:,spkgc[j]+1],'m',linewidth=2) -# if any(dgsp[:,1] == spkgc[j]): -# pylab.plot(dgvt[:,0],dgvt[:,spkgc[j]+1],'b',linewidth=2) -# else: -# pylab.plot(dgvt[:,0],dgvt[:,spkgc[j]+1],'m',linewidth=2) -pylab.ylabel('Potential (mV)') -#pylab.xticks(range(0,201,25)) -pylab.xlabel('Time (ms) - '+str(nspk)+' GCs spike') -#pylab.xlabel('Time (ms)') -pylab.xlim(0,tstop) -pylab.savefig('DG-'+str(sid)+'-'+str(nid)+idname+'.eps') -pylab.show() diff --git a/printfile.hoc b/printfile.hoc deleted file mode 100644 index 484d0bb..0000000 --- a/printfile.hoc +++ /dev/null @@ -1,170 +0,0 @@ -//********************** PRINT FILES ************************* - -// ModelDB file along with publication: -// Yim MY, Hanuschkin A, Wolfart J (2015) Hippocampus 25:297-308. -// http://onlinelibrary.wiley.com/doi/10.1002/hipo.22373/abstract - -// Man Yi Yim / 2015 - -// Print out the output spikes, membrane potential and network connections in the main code into files -// This code contains the following proc's -// proc saveNet() -// proc sMatrix_init() -// proc sMatrix() -// proc SpkMx() -// proc SpkMx_template() -// proc write_customstate() -// proc write_stimIn() - -// Read out the ID of pre or precell (by Ted Carnevale) -// argument is a NetCon - -obfunc myprecell() { localobj nil - // nil points to NULLobject - if ($o1.precell() == nil) { - return $o1.pre() - } else { - return $o1.precell() - } -} - -//***Network connections*** -strdef DGNC_name_ -objref dfile -dfile = new File() -proc saveNet(){ local i - print "write out Net cons file ..." - sprint (DGNC_name_,"%s-%d-%d%s.%s", "DGNC", PP_box_startid_, PP_box_nr_, idname, suffix) - dfile.wopen(DGNC_name_) - for i=0, nclist.count-1 { - dfile.printf("%s\t%s\n", myprecell(nclist.object(i)), nclist.object(i).postcell)} -// vbc2gc.x[i], vmc2gc.x[i], vhc2gc.x[i], vgc2gc.x[i])} ?? - dfile.close() -} - -//***Header of the voltage file*** -strdef strmat -objref efile -efile = new File() -strdef DGVt_name_ -proc sMatrix_init(){ // here the header info for sMatrix is written out... - n_to_print = 500 - print "Output memb voltage traces of every ", ngcell/n_to_print, " GC and all other cells to file" - sprint (DGVt_name_,"%s-%d-%d%s.%s", "DGVt", PP_box_startid_, PP_box_nr_, idname, suffix) - efile.wopen(DGVt_name_) - efile.printf("# t\t") - efile.printf("\n") - efile.close(DGVt_name_) -} - -proc sMatrix(){ local j - efile.aopen(DGVt_name_) - efile.printf("%f\t", t) - for i = 0, n_to_print-1 { - efile.printf("%f\t", cells.object[int(i*ngcell/n_to_print)].soma.v(0.5))} - efile.printf("\n") - efile.close(DGVt_name_) -} - -proc SpkMx() { local i, j - strdef DGsp_name_ - sprint (DGsp_name_,"%s-%d-%d%s.%s", "DGsp", PP_box_startid_, PP_box_nr_, idname, suffix) - print "Create: ", DGsp_name_ // output filename - dfile.wopen(DGsp_name_) - k = 0 - while (k < VmT.size) { - for j = 0, (ngcell+nbcell+nmcell+nhcell-1) { - if(Spike[j].x[k] != 0) { - dfile.printf("%f\t%d\n", VmT.x[k], j) // write out time of spike and cell id! - // Spike_times[j].append(VmT.x[k]) // generate vector of spike times for each cell. Moved to seperate function - } - } - k +=1 - } - dfile.close(DGsp_name_) -} - -proc SpkMx_template() { local i, j - strdef DGsp_name_ - - sprint (DGsp_name_,"%s-%d-%d%s.%s", "GCsp", PP_box_startid_, PP_box_nr_, idname, suffix) - - print "Create: ",DGsp_name_ // output filename - - dfile.wopen(DGsp_name_) - - k = 0 - while(k < VmT.size) { - for j = 0, ngcell-1 { - if(Spike[j].x[k] != 0) { - dfile.printf("%f\t%d\n", VmT.x[k], j) // write out time of spike and cell id! - } - } - k +=1 - } - dfile.close(DGsp_name_) -} - -proc write_customstate() { local i, j - strdef f_name_ - sprint (f_name_,"%s-%d-%d%s.%s", "NetworkState_OUT_", PP_box_startid_, PP_box_nr_, idname, suffix) - - print "Create: ",f_name_ // output filename - - dfile.wopen(f_name_) - - for i = 0, ngcell -1 { - forsec Gcell[i].all { // iterate over all sections of this cell - for (x,0) { // iterate over internal nodes - dfile.printf("%f\n", v) // write out time of spike and cell id! - } - } - } - for i = 0, nbcell -1 { - forsec Bcell[i].all { // iterate over all sections of this cell - for (x,0) { // iterate over internal nodes - dfile.printf("%f\n", v) // write out time of spike and cell id! - } - } - } - for i = 0, nmcell -1 { - forsec Mcell[i].all { // iterate over all sections of this cell - for (x,0) { // iterate over internal nodes - dfile.printf("%f\n", v) // write out time of spike and cell id! - } - } - } - for i = 0, nhcell -1 { - forsec Hcell[i].all { // iterate over all sections of this cell - for (x,0) { // iterate over internal nodes - dfile.printf("%f\n", v) // write out time of spike and cell id! - } - } - } - dfile.close(f_name_) -} - -proc write_stimIn() { - strdef f_name_ - sprint (f_name_,"%s-%d-%d%s.%s", "StimIn", PP_box_startid_, PP_box_nr_, idname, suffix) - - print "Create: ",f_name_ // output filename - - dfile.wopen(f_name_) - for i=0,npp-1 { - for j=0,vec_stim[i].size()-1 { - dfile.printf("%f\t%d\n",vec_stim[i].x[j],i) - } - } - - // write out noise input... - if (print_stim_noise == 1){ - if (debug_ == 2) {print "write out noise input..." } - for i=0,npp-1 { - for j=0,vec_stim_noise[i].size()-1 { - dfile.printf("%f\t%d\n",vec_stim_noise[i].x[j],i+npp) - } - } - } - dfile.close(f_name_) -} \ No newline at end of file diff --git a/ranstream.hoc b/ranstream.hoc deleted file mode 100644 index fa0cb95..0000000 --- a/ranstream.hoc +++ /dev/null @@ -1,23 +0,0 @@ -random_stream_offset_ = 1000 - -begintemplate RandomStream -public r, repick, start, stream -external random_stream_offset_ -external trial -objref r -proc init() { - stream = $1 - // print "init RandomStream instance" - r = new Random() - start() -} -func start() { - if (stream<5) { - // print "start Seed for Poisson Generator ",stream," with ",(stream*random_stream_offset_ + 1) - } - return r.MCellRan4(stream*random_stream_offset_ + 1 + trial) -} -func repick() { - return r.repick() -} -endtemplate RandomStream diff --git a/readme.html b/readme.html deleted file mode 100644 index 0b83817..0000000 --- a/readme.html +++ /dev/null @@ -1,129 +0,0 @@ - - -

Conductance-based computer model of the DG network realizing a spatiotemporal pattern separation task employing either physiological or leaky GC phenotype.

See the original paper for details:

-

Yim MY, Hanuschkin A, Wolfart J (2015) Intrinsic rescaling of granule cells restores pattern separation ability of a dentate gyrus network model during epileptic hyperexcitability. Hippocampus 25:297-308.

-

http://onlinelibrary.wiley.com/doi/10.1002/hipo.22373/abstract

-

 

-

Dr. Man Yi Yim / 2015

-

Dr. Alexander Hanuschkin / 2011

-

 

- - - -

To run the NEURON simulation and data analysis under Unix system:

-

 

-

1. Compile the mod files using the command

-

> nrnivmodl

-

 

-

2. Run the simulation (select the figure you want to simulate by setting fig in main.hoc before running)

-

> ./x86_64/special main.hoc

-

if your computer is running the 64-bit version, or

-

> ./i686/special main.hoc

-

for the 32-bit.

-

 

-

3. Open ipython or other command cells for Python, and run the data analysis

-

> ipython

-

> run fig1.py

-

 

-

Alternatively, you can set the idname of the following python codes and run the codes separately.

-

 

-

a) To plot the network activity in a trial (e.g. Fig 1C,D), run the python code plot_DG_all.py

-

fig1.jpg

-

 

-

b) To plot the activity of a neuron (e.g. Fig 1B), run the python code GCinput.py

-

fig2.jpg

-

 

-

c) To plot the network input and GC output (Fig 1E), run the python code inout_pattern.py

-

 

-

4. To make a scatter plot of similarity scores and fit the data (Fig 1E) , run the python code sim_score.py and then the matlab code FitSimScore_ForallFigs.m

-

fig3.jpg

-

 

- -

File description

-

Main code: run this code for the simulation

-

main.hoc

-

 

-

Printing code: format of the file output

-

printfile.hoc

-

 

-

Neuron models: morphology, conductances, ion channels and neuronal properties

-

GC.hoc

-

BC.hoc

-

MC.hoc

-

HIPP.hoc

-

 

-

Input models: properties of the inputs

-

PP.hoc

-

ranstream.hoc

-

 

-

Conductances: dynamics and properties of conductances

-

BK.mod

-

CaL.mod

-

CaN.mod

-

CaT.mod

-

ccanl.mod

-

HCN.mod

-

ichan2.mod

-

Ka.mod

-

Kir.mod

-

SK.mod

-

 

-

Spike generators:

-

netstimbox.mod

-

netstim125.mod

-

 

-

Python-Matlab-Analysis:

-

FitSimScore_ForallFigs.m fits the sim score data points by the method of least

-

squares.

-

plot_DG_all.py plots DG neurons' activity.

-

GCinput.py extracts and plots the inputs to a selected GC.

-

inout_pattern.py plots the inputs and GC outputs.

-

sim_score.py creates a scatter plot of output vs input sim scores. 

-

 

- - - -

Introduced changes in Mod files compared to the original DG model of Santhakumar et al. 2005

- -

In our scripts, the previously existing different potassium equilibrium potentials (Ekf, Eks, Ek..) were reduced to a single common Ek (e.g. GC.hoc, ichan2.mod, ....)). -

- -

CaL.mod
-CaN.mod
-CaT.mod
-These are new mod files for L-, N- and T-type calcium channels written by A. Hanuschkin following the description in Ca ion & L/T/N-Ca channels model of
-Aradi I, Holmes WR (1999) J Comput Neurosci 6:215-35.
-Note that eCa is calculated during simulation by ccanl.mod (see below). ecat, ecal values set in Santhakumar are not used in our model scripts.

- -

ccanl.mod
-Warning by Ted Carnevale 2015:
-The expression that this mechanism uses to calculate the contribution of ica to the rate of change of calcium concentration in the shell is
--ica*(1e7)/(depth*FARADAY)
-but it should really be
--ica*(1e7)/(depth*2*FARADAY)
-because the valence of ca is 2. The result of this omission is that the mechanism behaves as if the shell is only 1/2 as thick as the value specified by the depth parameter.

- -

ichan2.mod
-- added a tonic (leak) GABAA conductance to be modified during epilepsy (see Young CC, Stegen M, Bernard R, Muller M, Bischofberger J, Veh RW, Haas CA, Wolfart J (2009) J Physiol 587:4213-4233
-http://onlinelibrary.wiley.com/doi/10.1113/jphysiol.2009.170746/abstract)

- -

Kir.mod
-New Mod file
-Added an inward rectifier potassium conductance to be modied during epilepsy (see Young CC, Stegen M, Bernard R, Muller M, Bischofberger J, Veh RW, Haas CA, Wolfart J (2009) J Physiol 587:4213-4233)
-Channel description and parameters from:
-Stegen M, Kirchheim F, Hanuschkin A, Staszewski O, Veh R, and Wolfart J. Cerebral Cortex, 22:9, 2087-2101, 2012.
-http://cercor.oxfordjournals.org/content/22/9/2087.long

- -

SK.mod
-Correction: use of correct dynamics (see rate() lines: 95-101)

- - -

Other remarks

-

BK.mod
-Please note that cai was not assiged here in the original Santhakumar et al. (2005) version (which we used). It should be cai = ncai + lcai + tcai, as noted by
-Morgan RJ, Santhakumar V, Soltesz I (2007) Prog Brain Res 163:639-58
-The bug was fixed to make the channel properly dependent on the current calcium concentration. See
- https://senselab.med.yale.edu/modeldb/showModel.cshtml?model=124513&file=/dentate_gyrus/CaBK.mod

- - - \ No newline at end of file diff --git a/run.sh b/run.sh new file mode 100644 index 0000000..0ad2872 --- /dev/null +++ b/run.sh @@ -0,0 +1,24 @@ +# Clean up anything from previous runs +sh cleanup.sh + +# Make directory structure for outputs +mkdir data +mkdir data/dgnetwork +mkdir data/fi-curves + +mkdir figures +mkdir figures/raster-plots +mkdir figures/pattern-separation +mkdir figures/voltage-tracings +mkdir figures/fi-curves + +# Get FI Curves + +# Run DG Model +nrnivmodl mods +x86_64/special main.hoc + +# Plot Network Structure + +# Analyze Data +julia analysis.jl diff --git a/sim_score.py b/sim_score.py deleted file mode 100644 index 203185e..0000000 --- a/sim_score.py +++ /dev/null @@ -1,91 +0,0 @@ -### Analysis of DG network data ### -# This Python code creates a scatter plot of output vs input sim scores. -# Enter the idname - -# ModelDB file along with publication: -# Yim MY, Hanuschkin A, Wolfart J (2015) Hippocampus 25:297-308. -# http://onlinelibrary.wiley.com/doi/10.1002/hipo.22373/abstract - -# modified and augmented by -# Man Yi Yim / 2015 -# Alexander Hanuschkin / 2011 - -import pylab -import numpy -import matplotlib as mpl -font_size = 12 # 20 -mpl.rcParams['axes.titlesize'] = font_size+10 -mpl.rcParams['xtick.labelsize'] = font_size+6 -mpl.rcParams['ytick.labelsize'] = font_size+6 -mpl.rcParams['axes.labelsize'] = font_size+8 -mpl.rcParams['legend.fontsize'] = font_size -mpl.rcParams['font.size'] = font_size+10 - -sid = 0 -sidm = 12 -nid = 6 -step = 2 -delta = numpy.array([20])/2. -dur = 200. # duration of the simulation -bin = 1. -ninput = 100 # number of input sources -ncell = 500 # number of neurons (GCs) -#idname ='-pp10-gaba1-kir1-st0' - -def corr_score(file1,file2,delta,bin=1.,dur=100.,ncell=500): - """Similarity score by correlation coefficient. The spike trains are convolved with a triangular kernel.""" - d1 = numpy.loadtxt(file1) - d2 = numpy.loadtxt(file2) - x = numpy.zeros(int(ncell*dur/bin)) - y = numpy.zeros(int(ncell*dur/bin)) - for j in range(ncell): - if d1.size == 2: - s1 = numpy.array(d1[0]*(d1[1]==j)) - else: - s1 = d1[d1[:,1]==j,0] - if d2.size == 2: - s2 = numpy.array(d2[0]*(d2[1]==j)) - else: - s2 = d2[d2[:,1]==j,0] - kern = numpy.append(numpy.arange(delta/bin),numpy.arange(delta/bin,-1,-1)) - ts1,dump = pylab.histogram(s1,numpy.arange(0.,dur+bin,bin)) - ts2,dump = pylab.histogram(s2,numpy.arange(0.,dur+bin,bin)) - x[j*dur/bin:(j+1)*dur/bin] = numpy.convolve(ts1,kern,'same') - y[j*dur/bin:(j+1)*dur/bin] = numpy.convolve(ts2,kern,'same') - x = x - pylab.mean(x) - y = y - pylab.mean(y) - cor = sum(x*y)/(len(x)*pylab.std(x)*pylab.std(y)) - return cor - -pylab.figure() -fin1 = 'StimIn-'+str(sid)+'-'+str(nid)+idname+'.txt' -fout1 = 'GCsp-'+str(sid)+'-'+str(nid)+idname+'.txt' -inscore = numpy.empty((delta.size,sidm*(sidm+1)/2+1)) -outscore = numpy.empty((delta.size,sidm*(sidm+1)/2+1)) -for k in range(delta.size): - inscore[k,0] = corr_score(fin1,fin1,delta[k],bin,dur,ninput) - outscore[k,0] = corr_score(fout1,fout1,delta[k],bin,dur,ncell) -count = 1 -for j in range(0,sidm): - fin1 = 'StimIn-'+str(j)+'-'+str(nid)+idname+'.txt' - fout1 = 'GCsp-'+str(j)+'-'+str(nid)+idname+'.txt' - for m in range(j+1,sidm+1): - fin2 = 'StimIn-'+str(m)+'-'+str(nid)+idname+'.txt' - fout2 = 'GCsp-'+str(m)+'-'+str(nid)+idname+'.txt' - for k in range(delta.size): - inscore[k,count] = corr_score(fin1,fin2,delta[k],bin,dur,ninput) - outscore[k,count] = corr_score(fout1,fout2,delta[k],bin,dur,ncell) - count += 1 - -for k in range(delta.size): - pylab.plot(inscore[k,:],outscore[k,:],'ko',label='width='+str(2*delta[k])) -pylab.plot([0,1],[0,1],'k--',linewidth=2) -pylab.axis([-0.05,1.05,-0.05,1.05]) -pylab.xlabel('Input similarity score') -pylab.ylabel('Output similarity score') -pylab.savefig('sim_score_'+idname+'.eps') -f_write = open('sim_score'+idname+'.txt','w') -for j in range(inscore.size): - f_write.write(str(inscore[0][j])+'\t'+str(outscore[0][j])+'\n') -f_write.close() -pylab.show() diff --git a/utilities.jl b/utilities.jl new file mode 100644 index 0000000..c166029 --- /dev/null +++ b/utilities.jl @@ -0,0 +1,208 @@ +####################################################################### +# UTILITIES.JL +# Functions for analysis of data, I/O, plotting, etc. +# +####################################################################### +using DataFrames, CSV, Plots, PyCall, LsqFit, LaTeXStrings, Statistics +np = pyimport("numpy") +ss = pyimport("scipy.stats") + +####################################################################### +# I/O Functions +####################################################################### +""" + create_directories_if_not_exist() + +Creates necessary directories if they don't yet exist. +""" +function create_directories_if_not_exist() + if !("figures" ∈ readdir()) + mkdir("figures") + end + + for d ∈ ["raster-plots", "pattern-separation", "voltage-tracings", "fi-curves"] + if !(d ∈ readdir("figures")) + mkdir("figures/"*d) + end + end +end + +""" + load_spike_files(patterns, model_label, neuron_ids; neurons_per_pattern, data_dir) + +Loads spike time files and concatenates them. File names from simulation are structured as follows: + - "DGsp-{p}-{neurons_per_pattern}-{model_label}.txt" + - "StimIn-{p}-{neurons_per_pattern}-{model_label}.txt" + +The `neuron_ids` object is a `Dict` which includes lower and upper bound ranges on neuron IDs: + ``` + neuron_ids = Dict( + "DG" => [0, 500], + "BC" => [500,506], + "MC" => [506, 521], + "HIPP" => [521, 527] + ) + ``` +""" +function load_spike_files( + patterns::Union{UnitRange{Int64}, Vector{Int64}}, + model_label::String, + neuron_ids::Dict; + neurons_per_pattern::Int64=6, + data_dir::String="data/dgnetwork/") + + df = DataFrame() + for p ∈ patterns + spike_fname = data_dir*"DGsp-"*string(p)*"-"*string(neurons_per_pattern)*"-"*model_label*".txt" + stims_fname = data_dir*"StimIn-"*string(p)*"-"*string(neurons_per_pattern)*"-"*model_label*".txt" + spikes = CSV.read(spike_fname, delim="\t", header=0, DataFrame) + stimin = CSV.read(stims_fname, delim="\t", header=0, DataFrame) + + if size(spikes, 1) > 0 + rename!(spikes, ["Time", "Neuron"]) + spikes[:,"Population"] .= "" + spikes[:,"Pattern"] .= p + spikes[:,"NeuronsPerPattern"] .= neurons_per_pattern + spikes[:,"Model"] .= model_label + for k ∈ keys(neuron_ids) + lb, ub = neuron_ids[k] + spikes[lb .<= spikes[:, "Neuron"] .< ub, "Population"] .= k + end + + df = [df; spikes] + end + + if size(stimin, 1) > 0 + rename!(stimin, ["Time", "Neuron"]) + stimin[:,"Population"] .= "PP" + stimin[:,"Pattern"] .= p + stimin[:,"NeuronsPerPattern"] .= neurons_per_pattern + stimin[:,"Model"] .= model_label + df = [df; stimin] + end + end + + return df +end + +####################################################################### +# FUNCTIONS FOR STATISTICAL ANALYSIS +####################################################################### +""" + spike_train_correlation(spike_train1, spike_train2, n_neurons, delta, n_bins, duration) + +Computes the similarity score for two spike trains using the Pearson + correlation coefficient. This function was copied from that in + Yim et al. (2015), but modified in a few ways. We added documentation, + switched to using DataFrames (for clarity/brevity of code), and computed + the correlation coefficient using the `scipy.pearsonr` function. +""" +function spike_train_correlation( + spike_train1::DataFrame, + spike_train2::DataFrame, + n_neurons::Int64; + delta::Int64=20, + n_bins::Int64=1, + duration::Float64=200.0) + + # Create kernel + tri_rng = collect(1:round(Int64, delta/n_bins)) + triangular_kernel = [tri_rng; reverse(tri_rng)[2:end]]' + + # Create bins for histograms + bins = collect(1:n_bins:(duration + n_bins)) + + x = zeros(round(Int64, n_neurons * duration / n_bins)) + y = zeros(round(Int64, n_neurons * duration / n_bins)) + for i ∈ 1:n_neurons + s1 = spike_train1[spike_train1.Neuron .== (i-1), "Time"] + s2 = spike_train2[spike_train2.Neuron .== (i-1), "Time"] + timeseries1, _ = np.histogram(s1, bins) + timeseries2, _ = np.histogram(s2, bins) + idx_low = round(Int64, (i-1)*duration/n_bins + 1) + idx_high = round(Int64, i*duration/n_bins) + x[idx_low:idx_high] = np.convolve(timeseries1, triangular_kernel, "same") + y[idx_low:idx_high] = np.convolve(timeseries2, triangular_kernel, "same") + end + return ss.pearsonr(x, y) +end + + +""" + pattern_separation_curve(spike_times, n_pp_cells, n_granule_cells; delta, n_bins, duration) + +Computes pairwise correlations between input and output patterns +""" +function pattern_separation_curve( + spike_times::DataFrame, + n_pp_cells::Int64, + n_granule_cells::Int64; + delta::Int64=20, + n_bins::Int64=1, + duration::Float64=200.0) + + out = DataFrame() + n_patterns = unique(spike_times.Pattern) |> length + for i ∈ 0:(n_patterns - 2) + for j ∈ i:(n_patterns-1) + pp_a = spike_times[(spike_times.Population .== "PP") .& (spike_times.Pattern .== i), ["Neuron", "Time"]] + pp_b = spike_times[(spike_times.Population .== "PP") .& (spike_times.Pattern .== j), ["Neuron", "Time"]] + r_in, p_in = spike_train_correlation(pp_a, pp_b, n_pp_cells; delta=delta, n_bins=n_bins, duration=duration) + + gc_a = spike_times[(spike_times.Population .== "DG") .& (spike_times.Pattern .== i), ["Neuron", "Time"]] + gc_b = spike_times[(spike_times.Population .== "DG") .& (spike_times.Pattern .== j), ["Neuron", "Time"]] + r_out, p_out = spike_train_correlation(gc_a, gc_b, n_granule_cells; delta=delta, n_bins=n_bins, duration=duration) + + out_ij = DataFrame( + "Pattern A"=>[i], + "Pattern B"=>[j], + "Input Correlation"=>[r_in], + "Output Correlation"=>[r_out], + "Input Correlation p-Value"=>[p_in], + "Output Correlation p-Value"=>[p_out], + ) + out = [out; out_ij] + end + end + return out +end + +""" + fit_power_law(x, y) + +Fits `a + (1-a)*(x^b)` to the pattern separation data. +""" +function fit_power_law(x::Vector, y::Vector) + # Clip for tractability + x = max.(x, 0.0001) + y = max.(y, 0.0001) + + @. model(x, w) = w[1] + (1-w[1])*(x^w[2]) + res = curve_fit(model, x, y, [0., 1.]) + a, b = res.param + f(x) = a + (1-a)*(x^b) + return f +end + +####################################################################### +# FUNCTIONS FOR PLOTTING +####################################################################### + +""" + raster_plot(spikes; xlims, xlab, ylab) + +Returns a raster plot +""" +function raster_plot( + spikes::DataFrame; + xlims::Vector=[0,200], + xlab::String="Time (ms)", + ylab::String="Neuron ID" + ) + + p = plot(xlab=xlab, ylab=ylab, xlims=xlims, yticks=false) + p = scatter!(spikes[:,1], spikes[:,2], markershape=:vline, + label=nothing, c=:black, msc=:black, msa=1, ma=1, + markerstrokewidth=1) + return p +end \ No newline at end of file