Grid cell oscillatory interference with noisy network oscillators (Zilli and Hasselmo 2010)

 Download zip file 
Help downloading and running models
Accession:128812
To examine whether an oscillatory interference model of grid cell activity could work if the oscillators were noisy neurons, we implemented these simulations. Here the oscillators are networks (either synaptically- or gap-junction--coupled) of one or more noisy neurons (either Izhikevich's simple model or a Hodgkin-Huxley--type biophysical model) which drive a postsynaptic cell (which may be integrate-and-fire, resonate-and-fire, or the simple model) which should fire spatially as a grid cell if the simulation is successful.
Reference:
1 . Zilli EA, Hasselmo ME (2010) Coupled Noisy Spiking Neurons as Velocity-Controlled Oscillators in a Model of Grid Cell Spatial Firing J. Neurosci. 30(41):13850-13860
Citations  Citation Browser
Model Information (Click on a link to find other models with that property)
Model Type: Realistic Network; Neuron or other electrically excitable cell;
Brain Region(s)/Organism:
Cell Type(s): Neocortex spiny stellate cell; Abstract integrate-and-fire leaky neuron;
Channel(s): I Na,p; I Na,t; I K; I K,leak; I h;
Gap Junctions: Gap junctions;
Receptor(s):
Gene(s):
Transmitter(s):
Simulation Environment: MATLAB;
Model Concept(s): Oscillations; Synchronization; Simplified Models; Spike Frequency Adaptation; Grid cell;
Implementer(s): Zilli, Eric [zilli at bu.edu];
Search NeuronDB for information about:  I Na,p; I Na,t; I K; I K,leak; I h;
/
ZilliHasselmo2010
README.txt
Acker_sn_FI_n250.mat
fig_RS2sn_filthaftingtraj_n5000_p01_1noise_02res_July13_2D_spikes_0a.txt
fig_RS2sn_filthaftingtraj_n5000_p01_1noise_02res_July13_C_0a.mat
FigS4ABC_izh_simple0b_RS4sn_LIF_oct30_vary_gandncells_1noise_newpost_draft1.txt
FigS4DEF_izh_simple0b_RS4gn_LIF_oct30_vary_gandncells_1noise_newpost_draft1.txt
FigS8_izh_simple0b_RS2gn_LIF_nov08_vary_gandncells_0.125noise_newpost_draft1.txt
FigS8_izh_simple0b_RS2gn_LIF_nov08_vary_gandncells_0.25noise_newpost_draft1.txt
FigS8_izh_simple0b_RS2gn_LIF_nov08_vary_gandncells_0.5noise_newpost_draft1.txt
FigS8_izh_simple0b_RS2gn_LIF_nov08_vary_gandncells_1noise_newpost_draft1.txt
FigS8_izh_simple0b_RS2gn_LIF_nov08_vary_gandncells_2noise_newpost_draft1.txt
FigS8_izh_simple0b_RS2gn_LIF_nov08_vary_gandncells_4noise_newpost_draft1.txt
hafting_trajectory.m
rat_10925.mat
SI_acker_model_2d_grid.m
SI_acker_model_FI_history_dependence.m
SI_acker_model_FI_relations.m
SI_acker_model_stability_vs_params.m
SI_simple_model_2d_grid.m
SI_simple_model_FI_history_dependence.m
SI_simple_model_FI_relation.m
SI_simple_model_stability_vs_params.m
simple_model_RS1_FI_Jan09_n1.mat
simple_model_RS1gn_FI_Jan09_n250.mat
simple_model_RS1n_FI_Jan09_n1.mat
simple_model_RS1sn_FI_Jan09_n250.mat
simple_model_RS2_ext_FI_Jan09_n1.mat
simple_model_RS2_FI_Jan09_n1.mat *
simple_model_RS2gn_FI_Jan09_n250.mat
simple_model_RS2n_FI_Jan09_n1.mat
simple_model_RS2sn_FI_Jan01_n5000.mat
simple_model_RS2sn_FI_Jan01_n5000.txt
simple_model_RS2sn_FI_Jan09_n250.mat *
simple_model_RS2sn_FI_Jul10_n5000.mat
simple_model_RS2sn_FI_Jul10_n5000.txt
simple_model_RS2sn_FI_Jul11_n5000.mat
simple_model_RS2sn_FI_Jul11_n5000.txt
                            
% izhikevich simple model neuron network with random connectivity
% with integrate and fire postsynaptic cell
% eric zilli - june 27, 2009 using simple model MATLAB code from
% Izhikevich's book modified to be a network
%
% release version 1.0. check modeldb for updates.
%
% this source code is released into the public domain
%
% This script runs a network of optionally-noisy cells while the injected
% current is held constant. It allows the number of cells, the coupling
% type and strength, the connectivity probability, amount of noise, etc. to
% be changed.
%
% Note: this is a cleaned up version of the script used to generate the
% paper figures--it is possible errors were introduced during the cleaning!
% Feel free to contact me if there is any difficulty reproducing any
% results from the manuscript.

clear all
warning off

% type=1 - class 1 excitable, n=100, synaptic, varying g
% type=2 - class 1 excitable, n=100, gap-junction, varying g
% type=3 - (Fig S4A-C) class 2 excitable, n=250, synaptic, varying g
% type=4 - (Fig S4D-F) class 2 excitable, n=250, gap-junction, varying g
% type=5 - (Fig 5A,D) class 2 excitable, n = 1, ISI histograms
% type=6 - (Fig 5B,E) class 2 excitable, n = 250, gap-junction, ISI histograms
% type=7 - (Fig 5C,F) class 2 excitable, n = 250, synaptic, ISI histograms
% type=8 - (Fig S8) class 2 excitable, varying n and noise, gap-junction g=0.1
% type=9 - (Discussion) class 2 excitable, n=250, gap-junctions, 100% indep. noise, 0 corr. noise
% type=10 - (Discussion) class 2 excitable, n=250, gap-junctions, 95% indep. noise var, 5% corr. noise var
% type=11 - (Discussion) class 2 excitable, n=1, 5% realistic noise var
% type=12 - (Discussion) class 2 excitable, n=2500, gap-junctions, 5% realistic corr. noise var
% type=13 - class 2 excitable, n=250, inhibitory synaptic, varying g
% type=14 - (Fig 8D,H) class 2 excitable, n = 5000, p=0.01, synaptic, ISI histograms
for type=[1]
  % Note: some parameters may be overwritten in the 'if type' block:
  simdur = 8000; % (ms)
  dt = .1; % ms

  nruns = 1;

  % for statistics about each run
  fileOut = 0;
    
  % for ISI histograms:
  saveFigures = 0;
  numISIHists = 0;

  % these allow initial transients and such to be ignored:
  startAnalysisTime = dt; % (ms)
  endAnalysisTime = simdur; % (ms)

  % sets whether or not noise is used
  useNoise = 1;

  % approximate probability any two cells are connected (autoconnections will
  % be removed)
  % set to 0 (or set gs=0) to uncouple cells to, e.g., find noise level and
  % input level such that the median matches the biological data
  pcons = 1;

  weightmult = 1;
  
  % std. dev. of per-step noise term correlated among all network cells:
  commonNoiseSTD = 0*useNoise;
  commonNoise = 0;

  if type==1
    excitationClass = 1;
    useGapJunctions = 0;
    uniqueNoiseSTD = 130*useNoise;
    I = 62;
    % rough range of good gs = 150:1180
    gs = 150:50:1300;
    allncells = 100;
    nruns = 10;
  elseif type==2
    excitationClass = 1;
    useGapJunctions = 1;
    uniqueNoiseSTD = 130*useNoise;
    I = 62;
    % rough range of good gs = 0.022:16
    gs = 2.^[-5.5:0.5:4];
    allncells = 100;
    nruns = 10;
  elseif type==3
    excitationClass = 2;
    useGapJunctions = 0;
    uniqueNoiseSTD = 100*useNoise;
    I = 95.8;
    gs = 0:25:450;
    allncells = 250;
    nruns = 10;
  elseif type==4
    excitationClass = 2;
    useGapJunctions = 1;
    uniqueNoiseSTD = 100*useNoise;
    I = 95.8;
    gs = 2.^[-6:0.5:4];
    allncells = 250;
    nruns = 10;
  elseif type==5
    excitationClass = 2;
    useGapJunctions = 0;
    uniqueNoiseSTD = 100*useNoise;
    I = 95.8;
    gs = 0;
    allncells = 1;
    nruns = 1;
    simdur = 60000; % ms
    numISIHists = 1;
  elseif type==6
    excitationClass = 2;
    useGapJunctions = 1;
    uniqueNoiseSTD = 100*useNoise;
    I = 95.8;
    gs = 0.1;
    allncells = 250;
    nruns = 1;
    simdur = 60000; % ms
    numISIHists = 1;
  elseif type==7
    excitationClass = 2;
    useGapJunctions = 0;
    uniqueNoiseSTD = 100*useNoise;
    I = 95.8;
    gs = 256;
    allncells = 250;
    nruns = 1;
    simdur = 60000; % ms
    numISIHists = 1;
  elseif type==8
    excitationClass = 2;
    useGapJunctions = 1;
    %% run one at a time:
    uniqueNoiseSTD = 400*useNoise;
%     uniqueNoiseSTD = 200*useNoise;
%     uniqueNoiseSTD = 100*useNoise; % original level
%     uniqueNoiseSTD = 50*useNoise;
%     uniqueNoiseSTD = 25*useNoise;
%     uniqueNoiseSTD = 12.5*useNoise;
    I = 125;
    allncells = [16 32 64 128 256 512];
    gs = 0.1;
  elseif type==9
    excitationClass = 2;
    useGapJunctions = 1;
    uniqueNoiseSTD = 100*useNoise;
    commonNoiseSTD = 0*useNoise;
    I = 95.8;
    gs = 0.1;
    allncells = 250;
  elseif type==10
    excitationClass = 2;
    useGapJunctions = 1;
    uniqueNoiseSTD = sqrt(0.95)*100*useNoise;
    commonNoiseSTD = sqrt(0.05)*100*useNoise;
    I = 95.8;
    gs = 0.1;
    allncells = 250;
  elseif type==11
    excitationClass = 2;
    useGapJunctions = 1;
    uniqueNoiseSTD = 0*useNoise;
    commonNoiseSTD = sqrt(0.05)*100*useNoise;
    I = 95.8;
    gs = 0;
    allncells = 1;
  elseif type==12
    excitationClass = 2;
    useGapJunctions = 1;
    uniqueNoiseSTD = 0*useNoise;
    commonNoiseSTD = sqrt(0.05)*100*useNoise;
    I = 95.8;
    gs = 0.1;
    allncells = 2500;
  elseif type==13
    excitationClass = 2;
    useGapJunctions = 0;
    useNoise = 1;
    uniqueNoiseSTD = 100*useNoise;
    I = 110;
    gs = -(2.^[9:15]); % inhibitory synapses
    allncells = 1000;
    nruns = 1;
    simdur = 15000;
    weightmult = 8;
    endAnalysisTime=simdur;
elseif type==14
    excitationClass = 2;
    useGapJunctions = 0;
    uniqueNoiseSTD = 100*useNoise;
    I = 95.8;
    gs = 256;
    allncells = 5000;
    nruns = 1;
    simdur = 6000;0; % ms
    pcon = 0.01;
    numISIHists = 1;
end

  % file to save statistics for every individual cell in network, don't
  % save if indivname = ''
  indivname = '';
  % indivname = [strtok(outfilename,'.') '_individualcells.txt'];

  Cf=100; vr=-60; vt=-40; k=0.7;  % parameters used for RS
  a=0.03; c=-50; d=100;           % neocortical pyramidal neurons
  if excitationClass==1
    b=-2; % integrator
  elseif excitationClass==2
    b = 2; % resonator
  end
  vpeak=35; % spike cutoff

  % store min/med/max stability time of individual cells in each simulation
  mins = zeros(length(allncells),length(pcons),length(gs));
  medians = zeros(length(allncells),length(pcons),length(gs));
  maxs = zeros(length(allncells),length(pcons),length(gs));

  for ncellsi=1:length(allncells)
    ncells = allncells(ncellsi);
    
    typeString = '';
    if ncells>1
      if useGapJunctions
        typeString = [typeString 'g'];
      else
        typeString = [typeString 's'];
      end
    end
    if useNoise
      typeString = [typeString 'n'];
    end
    simprefix = sprintf('simple_model_RS%s_n%d_%s',typeString,ncells,datestr(now,'mmmdd'))

    if ~fileOut
      disp('no fileout')
    end
    outfilename = ['simple_model_RS' sprintf('%d',excitationClass) typeString];
    outfilename = [outfilename sprintf('_LIF_%s_vary_gandnoise_%gnoise.txt',datestr(now,'mmmdd'),useNoise)];

    for gi=1:length(gs)
      g = gs(gi);
      for pind=1:length(pcons)
        pcon = pcons(pind);

        stabilities = zeros(1,nruns);
        ISIstabilities = zeros(1,nruns);
        for run=1:nruns
          % activity of postsynaptic I&F
          post = zeros(1, round(simdur/dt));

          % spike times of postsynaptic I&F
          postspikes = spalloc(1, round(simdur/dt), 6*simdur); % expect firing at 6 Hz

          % LIF time constant
          tau = 4.5; % (ms)
          membraneDecay = exp(-dt/tau);
          postWeight=weightmult*3/ncells;

          % connectivity matrix
          if pcon>0
            C = (rand(ncells)<pcon)>0;
            C = sparse(C.*(1-eye(ncells))); % remove autoconnections
          else
            C = spalloc(ncells,ncells,1);
          end

          vhist = zeros(1,round(simdur/dt));

          t = 0; % current time in simulation

          v = vr*ones(ncells,1); % vr*rand(ncells,1);
          u = 0*v; % initial values

          % save state of one cell over simulation:
          state = zeros(2,simdur/dt);

          % binary array indicating which cells fire on each time steps
          spikes = sparse(ncells, simdur/dt,ncells*simdur); % sparse to save memory in big simulations

          tic
          while t<simdur
            t = t+dt; % advance to next time step
            tind = 1+round(t/dt);

            if commonNoiseSTD
              commonNoise = commonNoiseSTD*randn;
            end

            % if our LIF neuron is spiking constantly, let's just skip this
            % one (there is no apparent synchrony, which may either reflect the network
            % itself or the LIF parameters)
            if tind==500 && sum(postspikes)>(497)
              %             postspikes = NaN;
              break
            end

            oldv = v;
            oldu = u;
            if useGapJunctions
              if pcon==1
                % if all-to-all connected, the gap-junctions are equivalent
                % to sum(v)-ncells*v
                v = oldv + dt*(k*(oldv-vr).*(oldv-vt) - oldu + I + commonNoise + uniqueNoiseSTD*randn(ncells,1) + g*(sum(v)-ncells*v))/Cf;
              else
                % make matrix of differences in membrane potential
                vdiffs = C.*(repmat(oldv',ncells,1) - repmat(oldv,1,ncells));
                v = oldv + dt*(k*(oldv-vr).*(oldv-vt) - oldu + I + commonNoise + uniqueNoiseSTD*randn(ncells,1) + g*sum(vdiffs,2))/Cf;
                %% slowest way:
                % v = oldv + dt*(k*(oldv-vr).*(oldv-vt) - oldu + I + commonNoise + uniqueNoiseSTD*randn(ncells,1) + g*sum(C.*vdiffs,2))/Cf;
              end
            else
              if pcon==1
                v = oldv + dt*(k*(oldv-vr).*(oldv-vt) - oldu + I + commonNoise + uniqueNoiseSTD*randn(ncells,1) + g*(sum(spikes(:,tind-1),1)-spikes(:,tind-1)))/Cf; 
              else
                v = oldv + dt*(k*(oldv-vr).*(oldv-vt) - oldu + I + commonNoise + uniqueNoiseSTD*randn(ncells,1) + g*(C*spikes(:,tind-1)))/Cf;
              end
            end
            u = oldu + dt*a*(b*(oldv-vr)-oldu);

            vhist(tind) = v(1);

            % save and reset spikes when v>=vpeak
            spikes(:,tind) = v>=vpeak;
            u(v>=vpeak) = u(v>=vpeak)+d;
            oldv(v>=vpeak) = vpeak;
            v(v>=vpeak) = c;

            state(:,tind) = [oldv(1); oldu(1)];

            post(:,tind) = post(:,tind-1)*membraneDecay + tau*postWeight*sum(spikes(:,tind));
            postspikes(:,tind) = post(:,tind)>1;
            post(postspikes(:,tind)>1,tind) = 0; % reset to 0
          end
          toc

          clear C

          % calculate mean and variance of ISIs for each postsynaptic cell
          npost = 1;
          ISImeans = zeros(1,npost);
          ISIstds = zeros(1,npost);
          ISIstabilities = zeros(1,nruns);

          %% calculate ISI statistics for each individual unit in the
          %% network:
          indivmeans = zeros(1,ncells);
          indivstds = zeros(1,ncells);
          indivstabilities = zeros(1,ncells);
          for ce=1:ncells % only run this once because we only save the last value in stabilities
            spikeTimes = find(spikes(ce,(startAnalysisTime/dt):(endAnalysisTime/dt))>0);
            ISIs = dt*diff(spikeTimes)/1000; % convert to seconds

            % count ISIs as occuring between initial
            % spikes of bursts (and additional spikes as occuring in a given
            % burst if they are < 50 ms apart);
            fISIs = [];
            csum = 0;
            for is=1:length(ISIs)
              csum = csum + ISIs(is);
              if ISIs(is)>.050
                fISIs = [fISIs csum];
                csum=0;
              end;
            end
            ISIs = fISIs;

            if ce<=numISIHists
              figure; hist(ISIs,linspace(.1,.3,200))
              title(['VCO Cell #'  num2str(ce) ' ISIs'],'fontsize',12)
              xlabel('ISI duration (s)','fontsize',12)
              set(gca,'fontsize',12);
              set(gca,'box','off');
              set(gcf,'position',[82 404 257 185]);
              fname = [simprefix '_vco_isihist_' num2str(ce)];
              if saveFigures && exist([fname '.eps'])~=2
                print_eps([fname '.eps'])
                saveas(gcf,[fname '.fig'])
              end
            end


            % if we have 3 or fewer ISIs, let's just trash this cell:
            if length(ISIs)<=3
              indivmeans(ce) = NaN;
              indivstds(ce) = NaN;
              indivstabilities(ce) = NaN;
            else
              indivmeans(ce) = mean(ISIs);
              indivstds(ce) = std(ISIs);
              indivstabilities(ce) = 5*indivmeans(ce)^3/16/pi/pi/(eps+indivstds(ce))/(eps+indivstds(ce));
            end
          end
          % drop NaNs and sort:
          indivmeans(isnan(indivmeans)) = [];
          indivstds(isnan(indivstds)) = [];
          indivstabilities(isnan(indivstabilities)) = [];
          [indivstabilities,ix] = sort(indivstabilities);
          indivmeans = indivmeans(ix);
          indivstds = indivstds(ix);


          %% calculate ISI statistics for the postsynaptic LIF neurons
          for ce=1:npost % only run this once because we only save the last value in stabilities
            spikeTimes = find(postspikes(ce,(startAnalysisTime/dt):(endAnalysisTime/dt))>0);
            ISIs = dt*diff(spikeTimes)/1000; % convert to seconds

            % count ISIs as occuring between initial
            % spikes of bursts (and additional spikes as occuring in a given
            % burst if they are < 50 ms apart);
            fISIs = [];
            csum = 0;
            for is=1:length(ISIs)
              csum = csum + ISIs(is);
              if ISIs(is)>.050
                fISIs = [fISIs csum];
                csum=0;
              end;
            end
            ISIs = fISIs;

            if ce==1 && numISIHists
              figure; hist(ISIs,linspace(.1,.3,200))
              title('Postsynaptic LIF ISIs','fontsize',12)
              xlabel('ISI duration (s)','fontsize',12)
              set(gca,'fontsize',12)
              set(gca,'box','off')
              set(gcf,'position',[82 404 257 185]);
              fname = [simprefix '_post_isihist'];
              if saveFigures && exist([fname '.eps'])~=2
                print_eps([fname '.eps'])
                saveas(gcf,[fname '.fig'])
              end
            end

            ISImeans(ce) = mean(ISIs);
            ISIstds(ce) = std(ISIs);
            ISIstabilities(ce) = 5*ISImeans(ce)^3/16/pi^2/(eps+ISIstds(ce))^2;
          end

          % save to file:
          if fileOut
            fid = fopen(outfilename,'a');
            fprintf(fid,'%d\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%d\n', ...
              ncells,pcon,I,gs(gi),simdur/1000, useNoise*uniqueNoiseSTD,tau,postWeight,...
              Cf,vr,vt,k,a,b,c,d,vpeak,...
              1/ISImeans(1),ISImeans(1),ISIstds(1),ISIstabilities(1),length(ISIs),...
              min(indivmeans),median(indivmeans),max(indivmeans),...
              min(indivstds),median(indivstds),max(indivstds),...
              min(indivstabilities),median(indivstabilities),max(indivstabilities),...
              length(find(postspikes>0)));
            fclose(fid);
          end
          % print to screen:
          fprintf('%d\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%d\n', ...
            ncells,pcon,I,gs(gi),simdur/1000, useNoise*uniqueNoiseSTD,tau,postWeight,...
            Cf,vr,vt,k,a,b,c,d,vpeak,...
            1/ISImeans(1),ISImeans(1),ISIstds(1),ISIstabilities(1),length(ISIs), ...
            min(indivmeans),median(indivmeans),max(indivmeans),...
            min(indivstds),median(indivstds),max(indivstds),...
            min(indivstabilities),median(indivstabilities),max(indivstabilities),...
            length(find(postspikes>0)));

          % optinally save individual cell statistics to a file:
          if ~isempty(indivname) && fileOut
            fid = fopen(indivname,'a');
            fprintf(fid,'%d\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%d', ...
              ncells,pcon,I,gs(gi),simdur/1000, useNoise*uniqueNoiseSTD,tau,postWeight,...
              1/ISImeans(1),ISImeans(1),ISIstds(1),ISIstabilities(1),length(ISIs),...
              length(find(postspikes>0)));
            for ce=1:length(indivmeans)
              fprintf(fid,'\t%g',indivmeans(ce));
            end
            for ce=1:ncells-length(indivmeans) % for alignment in ouput, put back the NaNs we removed
              fprintf(fid,'\tNaN');
            end
            fprintf(fid,'\t');
            for ce=1:length(indivmeans)
              fprintf(fid,'\t%g',indivstds(ce));
            end
            for ce=1:ncells-length(indivmeans)
              fprintf(fid,'\tNaN');
            end
            fprintf(fid,'\t');
            for ce=1:length(indivmeans)
              fprintf(fid,'\t%g',indivstabilities(ce));
            end
            for ce=1:ncells-length(indivmeans)
              fprintf(fid,'\tNaN');
            end
            fprintf(fid,'\n');
            fclose(fid);
          end
        end

        medians(ncellsi,pind,gi) = median(stabilities);
        maxs(ncellsi,pind,gi) = max(stabilities);
        mins(ncellsi,pind,gi) = min(stabilities);
      end
    end
  end % end ncells loop
end % end run loop

warning on
return

%% Figure 1
data = textread('simple_model_RS2sn_vary_gandncells_1noise.txt');
figure; plot(data(:,4),data(:,18),'.');
set(gca,'xlim',[0 500])
set(gca,'ylim',[0 13])
set(gca,'fontsize',12,'box','off')
xlabel('Synaptic conductance','fontsize',12); ylabel('Network frequency (Hz)','fontsize',12)
set(gcf,'position',[82 404 257 302])
if exist('fig_freq_varyg.eps')~=2
  print_eps('fig_freq_varyg.eps')
  saveas(gcf,'fig_freq_varyg.fig')
end
figure; plot(data(:,4),data(:,21),'.'); xlabel('Synaptic conductance','fontsize',12); ylabel('Estimated stability time (s)','fontsize',12)
set(gcf,'position',[82 404 257 302])
set(gca,'xlim',[0 500])
set(gca,'ylim',[0 350])
set(gca,'fontsize',12,'box','off')
if exist('fig_stability_varyg.eps')~=2
  print_eps('fig_stability_varyg.eps')
  saveas(gcf,'fig_stability_varyg.fig')
end
figure; plot(data(:,4),data(:,20),'.'); xlabel('Synaptic conductance','fontsize',12); ylabel('Network period std. dev. (s)','fontsize',12)
set(gcf,'position',[82 404 257 302])
set(gca,'xlim',[0 500])
set(gca,'yscale','log')
set(gca,'fontsize',12,'box','off')
if exist('fig_stdev_varyg.eps')~=2
  print_eps('fig_stdev_varyg.eps')
  saveas(gcf,'fig_stdev_varyg.fig')
end
data = textread('simple_model_RS2gn_vary_gandncells_1noise.txt');
figure; plot(data(:,4),data(:,18),'.');
set(gca,'xlim',[1e-2 10])
set(gca,'xtick',[.01 .1 1 10])
set(gca,'xscale','log')
set(gca,'fontsize',12,'box','off')
% set(gca,'ylim',[0 9.1])
xlabel('Gap-junction conductance','fontsize',12); ylabel('Network frequency (Hz)','fontsize',12)
set(gcf,'position',[82 404 257 302])
if exist('fig_freq_varyg.eps')~=2
  print_eps('fig_freq_varyg.eps')
  saveas(gcf,'fig_freq_varyg.fig')
end
figure; plot(data(:,4),data(:,21),'.'); xlabel('Gap-junction conductance','fontsize',12); ylabel('Estimated stability time (s)','fontsize',12)
set(gcf,'position',[82 404 257 302])
set(gca,'xlim',[1e-2 20])
set(gca,'xtick',[.01 .1 1 10])
set(gca,'xscale','log')
set(gca,'fontsize',12,'box','off')
if exist('fig_stability_varyg.eps')~=2
  print_eps('fig_stability_varyg.eps')
  saveas(gcf,'fig_stability_varyg.fig')
end
figure; plot(data(:,4),data(:,20),'.'); xlabel('Gap-junction conductance','fontsize',12); ylabel('Network period std. dev. (s)','fontsize',12)
set(gcf,'position',[82 404 257 302])
set(gca,'xlim',[1e-2 20])
set(gca,'xtick',[.01 .1 1 10])
set(gca,'xscale','log')
set(gca,'yscale','log')
set(gca,'fontsize',12,'box','off')
if exist('fig_stdev_varyg.eps')~=2
  print_eps('fig_stdev_varyg.eps')
  saveas(gcf,'fig_stdev_varyg.fig')
end


%% Figure 2
% histogram plots are buried in the script above, search for numISIHists


%% Figure 3
% first run type=8 for each standard deviation of the noise, saving output
% we expect all files have the same number of rows
nfiles = 6; % i.e. number of noise levels
startrow = 0;
endrow = 6;

ISIfreqs = zeros(endrow-startrow,nfiles);
ISImeans = zeros(endrow-startrow,nfiles);
ISIstds = zeros(endrow-startrow,nfiles);
ISIstability = zeros(endrow-startrow,nfiles);

data = textread('simple_model_RS2gn_LIF_vary_gandnoise_4noise.txt');
ISIfreqs(:,1) = data((startrow+1):endrow,18);
ISImeans(:,1) = data((startrow+1):endrow,19);
ISIstds(:,1) = data((startrow+1):endrow,20);
ISIstability(:,1) = data((startrow+1):endrow,21);
data = textread('simple_model_RS2gn_LIF_vary_gandnoise_2noise.txt');
ISIfreqs(:,2) = data((startrow+1):endrow,18);
ISImeans(:,2) = data((startrow+1):endrow,19);
ISIstds(:,2) = data((startrow+1):endrow,20);
ISIstability(:,2) = data((startrow+1):endrow,21);
data = textread('simple_model_RS2gn_LIF_vary_gandnoise_1noise.txt');
ISIfreqs(:,3) = data((startrow+1):endrow,18);
ISImeans(:,3) = data((startrow+1):endrow,19);
ISIstds(:,3) = data((startrow+1):endrow,20);
ISIstability(:,3) = data((startrow+1):endrow,21);
data = textread('simple_model_RS2gn_LIF_vary_gandnoise_0.5noise.txt');
ISIfreqs(:,4) = data((startrow+1):endrow,18);
ISImeans(:,4) = data((startrow+1):endrow,19);
ISIstds(:,4) = data((startrow+1):endrow,20);
ISIstability(:,4) = data((startrow+1):endrow,21);
data = textread('simple_model_RS2gn_LIF_vary_gandnoise_0.25noise.txt');
ISIfreqs(:,5) = data((startrow+1):endrow,18);
ISImeans(:,5) = data((startrow+1):endrow,19);
ISIstds(:,5) = data((startrow+1):endrow,20);
ISIstability(:,5) = data((startrow+1):endrow,21);
data = textread('simple_model_RS2gn_LIF_vary_gandnoise_0.125noise.txt');
ISIfreqs(:,6) = data((startrow+1):endrow,18);
ISImeans(:,6) = data((startrow+1):endrow,19);
ISIstds(:,6) = data((startrow+1):endrow,20);
ISIstability(:,6) = data((startrow+1):endrow,21);

ncellsvect = data((startrow+1):endrow,1);

figure; imagesc(log10(ISIstability(1:end,:))');
set(gca,'yticklabel',[400 200 100 50 25 12.5],'xtick',[1:6],'xticklabel',[16 32 64 128 256 512]);
xlabel('Number of cells','fontsize',12); ylabel('Noise standard deviation','fontsize',12);
title('Estimated stability time (log s)','fontsize',12)
set(gca,'fontsize',12)
h=colorbar; set(h,'fontsize',12);
colormap(flipud(gray));
set(gcf,'position',[82 404 322 275])
if exist('fig_RS2gn_varynnoise.eps')~=2
  print_eps('fig_RS2gn_varyn.eps')
  saveas(gcf,'fig_RS2gn_varyn.fig')
end