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
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
                            
% Acker et al 2003 model variant neuron network with random connectivity
% with integrate and fire postsynaptic cell
% eric zilli - oct 10, 2009
%
% release version 1.0. check modeldb for updates.
%
% this source code is released into the public domain
%
% The Acker et al model is a biophysical model of an entorhinal cortex
% layer II stellate. The model uses voltage-shifted Hodghkin-Huxley loligo
% Na and K channels for spiking, plus a non-inactivating persistent-sodium
% current and a two-component (fast and slow) H-current which are modified
% by Acker et al. from original current models by Erik Fransen. We modify the model by
% scaling down all the conductance densities by increasing the membrane
% capacitance. This slows the minimum firing frequency of the model down so
% that the noise level can be scaled to match our biological data (ISI mean
% 0.2 s, ISI std. dev. 0.036 s).
%
% Interesting questions to continue this line of research: Can we calculate
% these statistics analytically? Approximately? Can we predict the
% statistics of a single cell in a network based on numerical measurements
% of a single uncoupled cell?
%
% 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

simdur = 4000; % (ms)
dt = .01; % ms

nruns = 1;

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

useNoise = 1;
useFrozenNoise = 0; % requires more memory
spikeThreshold = 20; % mV

% if true, use a slower version of the H current
slowslow = 0;

% if true, analyze only the last half of the ISIs
lastHalf = 1;

useGapJunctions = 0;
if ~slowslow
  uniqueNoiseSTD = 3.44*useNoise;
else
  uniqueNoiseSTD = 10*useNoise;
end

% approximate probability any two cells are connected (autoconnections will
% be removed)
pcons = 1;

% type = 1 n=150 uncoupled, set I and noise to match biology
% type = 2 n=250, synaptic, varying g (to find coupling for Figure 10)
% type = 3 n=250, gap-junction, varying g
for type=[3]
  if type==1
    useGapJunctions = 0;
    allncells = 150;
    pcons = 0;
    gs = 0;
    I = -2.37;
    simdur = 4000; % ms
    useNoise = 1;
    uniqueNoiseSTD = 3.44*useNoise;
  elseif type==2
    useGapJunctions = 0;
    allncells = 250;
    pcons = 1;
    gs = 2.^[3:.5:8];
    I = -2.37;
    simdur = 4000; % ms
    useNoise = 1;
    uniqueNoiseSTD = 3.44*useNoise;
  elseif type==3
    useGapJunctions = 1;
    allncells = 250;
    pcons = 1;
    gs = 2.^[-6:.5:4];
    I = -2.37;
    simdur = 4000; % ms
    useNoise = 1;
    uniqueNoiseSTD = 3.44*useNoise;
  end
end

fileOut = 0;
if ~fileOut
  'no fileout'
end
outfilename = 'acker_model_';
if ncells>1
  if useGapJunctions
    outfilename = [outfilename 'g'];
  else
    outfilename = [outfilename 's'];
  end
  endif useNoise
  outfilename = [outfilename 'n'];
end
outfilename = [outfilename sprintf('_vary_gncellspcon_%gnoise_%s.txt',useNoise,datestr(now,'mmmmdd'))];

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

Cm = 5; 1.5; % uF/cm^2
gNa = 52; % mS/cm^2
gNap = 0.5; % mS/cm^2
gH = 1.5; % mS/cm^2
gK = 11; % mS/cm^2
gLeak = 0.5; % mS/cm^2
EH = -20; % (mv)
ELeak = -65; % (mv)
EK = -90; % (mv)
ENa = 55; % (mv)

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);
  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

        % I&F time constant
        tau = 4.5; % (ms)
        membraneDecay = exp(-dt/tau);
        postWeight=1.2/ncells;

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

        state = zeros(7,round(simdur/dt));

        if useFrozenNoise
          randn('seed',3);
          frozenNoise = randn(ncells,round(simdur/dt)+1);
        end

        % magnitude of noise added to all cells on each step
        commonNoiseSTD = 0;

        t = 0; % current time in simulation
        tind = 1;

        % initial values
        v = -55*ones(1,ncells);
        mNa = zeros(1,ncells);
        hNa = zeros(1,ncells);
        mNap = zeros(1,ncells);
        nK = zeros(1,ncells);
        mHfast = zeros(1,ncells);
        mHslow = zeros(1,ncells);
        state(:,1) = [v(1); mNa(1); hNa(1); mNap(1); nK(1); mHfast(1); mHslow(1)];

        % 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 = tind+1; 1+round(t/dt);

          % 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;
          alphamNa = -0.1*(oldv+23)./(exp(-0.1*(oldv+23))-1);
          betamNa = 4*exp(-(oldv+48)/18);
          mNa = mNa + dt*(alphamNa.*(1-mNa) - betamNa.*mNa);

          alphahNa = 0.07*exp(-(oldv+37)/20);
          betahNa = 1./(exp(-0.1*(oldv+7))+1);
          hNa = hNa + dt*(alphahNa.*(1-hNa) - betahNa.*hNa);

          alphanK = -0.01*(oldv+27)./(exp(-0.1*(oldv+27))-1);
          betanK = 0.125*exp(-(oldv+37)/80);
          nK = nK + dt*(alphanK.*(1-nK) - betanK.*nK);

          alphamNap = 1./(0.15*(1+exp(-(oldv+38)/6.5)));
          betamNap = exp(-(oldv+38)/6.5)./(0.15*(1+exp(-(oldv+38)/6.5)));
          mNap = mNap + dt*(alphamNap.*(1-mNap) - betamNap.*mNap);

          minfHfast = 1./(1+exp((oldv+79.2)/9.78));
          mtauHfast = 0.51./(exp((oldv-1.7)/10) + exp(-(oldv+340)/52)) + 1;
          mHfast = mHfast + dt*((minfHfast-mHfast)./mtauHfast);

          minfHslow = 1./(1+exp((oldv+71.3)/7.9));
          if ~slowslow
            mtauHslow = 5.6./(exp((oldv-1.7)/14) + exp(-(oldv+260)/43)) + 1;
          else
            mtauHslow = 65.5813 + 248.0469*exp(-(-79.2190-oldv).^2/33.5178^2); % orig
          end
          mHslow = mHslow + dt*((minfHslow-mHslow)./mtauHslow);

          if useFrozenNoise
            vnoise = frozenNoise(:,tind);
          else
            vnoise = uniqueNoiseSTD*randn(1,ncells);
          end
          if useGapJunctions
            if pcon==0
              v = v + dt*(vnoise - gH*(0.65*mHfast + 0.35*mHslow).*(v-EH) - (gNap*mNap+gNa*mNa.^3.*hNa).*(v-ENa) - gK*nK.^4.*(v-EK) - gLeak*(v-ELeak) + I)/Cm;
            elseif pcon==1
              v = v + dt*(vnoise - gH*(0.65*mHfast + 0.35*mHslow).*(v-EH) - (gNap*mNap+gNa*mNa.^3.*hNa).*(v-ENa) - gK*nK.^4.*(v-EK) - gLeak*(v-ELeak) + I + g*(sum(v)-ncells*v))/Cm;
            else
              vdiffs = repmat(v',ncells,1) - repmat(v,1,ncells);
              v = v + dt*(vnoise - gH*(0.65*mHfast + 0.35*mHslow).*(v-EH) - (gNap*mNap+gNa*mNa.^3.*hNa).*(v-ENa) - gK*nK.^4.*(v-EK) - gLeak*(v-ELeak) + I + g*sum(C*vdiffs,2))/Cm;
            end
          else
            if pcon==0
              v = v + dt*(vnoise - gH*(0.65*mHfast + 0.35*mHslow).*(v-EH) - (gNap*mNap+gNa*mNa.^3.*hNa).*(v-ENa) - gK*nK.^4.*(v-EK) - gLeak*(v-ELeak) + I)/Cm;
            elseif pcon==1
              v = v + dt*(vnoise - gH*(0.65*mHfast + 0.35*mHslow).*(v-EH) - (gNap*mNap+gNa*mNa.^3.*hNa).*(v-ENa) - gK*nK.^4.*(v-EK) - gLeak*(v-ELeak) + I + g*(sum(spikes(:,tind-1),1)-spikes(:,tind-1)'))/Cm;
            else
              v = v + dt*(vnoise - gH*(0.65*mHfast + 0.35*mHslow).*(v-EH) - (gNap*mNap+gNa*mNa.^3.*hNa).*(v-ENa) - gK*nK.^4.*(v-EK) - gLeak*(v-ELeak) + I + g*(C*spikes(:,tind-1))')/Cm;
            end
          end
          state(:,tind) = [v(1); mNa(1); hNa(1); mNap(1); nK(1); mHfast(1); mHslow(1)];

          % spikes if upward crossing spikeThreshold mV
          spikes(:,tind) = (oldv<spikeThreshold).*(v>=spikeThreshold);

          post(:,tind) = post(:,tind-1)*membraneDecay + tau*postWeight*sum(spikes(:,tind));
          % faster(?) if we don't need postspikes:
          post(post(:,tind)>1,tind) = -1e-10;
          % slower way if postspikes is needed:
          %           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

          % instead of dropping, 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 lastHalf && length(ISIs)
            ISIs = ISIs(end/2: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^2/(eps+indivstds(ce))^2;
          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);
          spikeTimes = find(post(ce,(startAnalysisTime/dt):(endAnalysisTime/dt))==-1e-10);
          ISIs = dt*diff(spikeTimes)/1000; % convert to seconds

          % instead of dropping, 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;

          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%d\n', ...
            ncells,pcon,I,gs(gi),simdur/1000, useNoise*uniqueNoiseSTD,tau,postWeight,...
            spikeThreshold,...
            mean(1./ISIs),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
        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%d\n', ...
          ncells,pcon,I,gs(gi),simdur/1000, useNoise*uniqueNoiseSTD,tau,postWeight,...
          spikeThreshold,...
          mean(1./ISIs),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)));

        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,...
            mean(1./ISIs),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

warning on

Loading data, please wait...