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 time varying input with LIF-measured stability
% eric zilli - oct ??, 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 research: Setting slowslow=1
% decreases the adaptation's initial magnitude but the effect
% has an effect lasting many ISIs long. Is
% the corresponding time constant model realistic? Experiments looking at H
% current time constants at spiking voltages do not appear to exist, but
% gating variable "friction" in the extended HH model accounts for a
% possible non-zero base for time constants which might play a major role
% here. In-vitro experiments are needed.
% Under what conditions does the ISI period or frequency overshoot and when does it undershoot?
% What is the best way to quantify the history dependence on firing rate
% effect?
% Is it possible for the intrinsic properties of one neuron A to be exactly
% the inverse of B in the sense of this adaptation such that an
% input signal I to A would be distorted by the adaptation into an
% output I' which would then be un-distorted by B to recover a constant function
% of I? (E.g. I goes from driving at 7 Hz to 11 Hz, A's firing rate adapts from
% 11.5 Hz down to 11 Hz. Can there be a B receiving A's spikes that fires at 11 Hz
% given input at 11.5 Hz and which adapts at just the right rate that B continues to
% fire at 11 Hz as the input from A slows down slightly?)
%
% This script starts a network with a specified current level (or
% frequency) and after a specified number of spikes allows the current
% level to be changed. This can be done with a single starting and ending
% current level (Figure S2A-D), a single starting and multiple ending levels
% (Figure S2E-F), multiple starting and single ending level, or a grid of injected
% current levels (Figure S3). These correspond to three types of figures that
% might be of interest.
%
% If the FI curve of the cell or network is known, it can be used to find
% the appropriate current levels. Otherwise they must be set by hand. The
% FI curve needn't be terribly high in resolution, so generating them as
% needed should be quick and easy.
%
% 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;

% type=1 - (Figure S2B,D) class 2 excitable, n=1, noise-free, 1 start, 1 end
% type=2 - (Figure S2F) class 2 excitable, n=1, noise-free, 1 start, 5 end, low-to-high
% type=3 - class 2 excitable, n=1, noise-free, 7:0.5:11 Hz start, 7:0.5:11 Hz end
for type=[3]
  simdur = 10000; % (ms)
  dt = .01; % ms

  nruns =1;

  saveFigures = 0;

  useNoise = 0;
  spikeThreshold = 20; % mV
  slowslow = 0;

  numDropISIs = 10;
  numPreISIs = 10;
  numPostISIs = 10;

  if type==1
    excitationClass = 2;
    useGapJunctions = 1;
    g = 0;
    uniqueNoiseSTD = 0*useNoise;
    ncells = 1;
    pcon = 0;
    load Acker_FI_n1.mat;
    inputMags = interp1(freqs,currents,7,'cubic');
    inputMags2 = interp1(freqs,currents,11,'cubic');
  elseif type==2
    excitationClass = 2;
    useGapJunctions = 1;
    g = 0;
    uniqueNoiseSTD = 0*useNoise;
    ncells = 1;
    pcon = 0;
    load Acker_FI_n1.mat;
    desiredFreqs = 7:11;
    inputMags = interp1(freqs,currents,7,'cubic');
    inputMags2 = interp1(freqs,currents,desiredFreqs,'cubic');
  elseif type==3
    excitationClass = 2;
    useGapJunctions = 1;
    g = 0;
    uniqueNoiseSTD = 0*useNoise;
    ncells = 1;
    pcon = 0;
    load Acker_FI_n1.mat;
    desiredFreqs = [7 11];
    inputMags = interp1(freqs,currents,desiredFreqs,'cubic');
    inputMags2 = interp1(freqs,currents,desiredFreqs,'cubic');
  end

  fileOut = 0;
  typePrefix = [sprintf('%d',excitationClass)];
  if ncells>1
    if useGapJunctions
      typePrefix = [typePrefix 'g'];
    else
      typePrefix = [typePrefix 's'];
    end
  end
  if useNoise
    typePrefix = [typePrefix 'n'];
  end
  outfilename = sprintf('acker_model_%s_history_dependence_%dn_%gnoise_%s.txt',typePrefix,ncells,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)

  slopes = zeros(length(inputMags));
  if length(inputMags)==1 && length(inputMags2)>1
    figISIs = zeros(length(inputMags2),numPreISIs+numPreISIs);
  elseif length(inputMags)>1 && length(inputMags2)==1
    figISIs = zeros(length(inputMags),numPreISIs+numPreISIs);
  elseif length(inputMags)>1 && length(inputMags2)>1
    figISIs = zeros(length(inputMags),length(inputMags2));
  end

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

  for inputi=1:length(inputMags)
    for inputi2=1:length(inputMags2)

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

        % spike times of postsynaptic I&F
        postspikes = 0;

        % I&F params
        tau = 2.5; % (sec)
        membraneDecay = exp(-dt/tau);
        postWeight=4/ncells;

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

        numISIs = -1;

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

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

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

        I = inputMags(inputi);
        Is = zeros(1,round(simdur/dt));

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

          %% done below now:
          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);
          end
          mHslow = mHslow + dt*((minfHslow-mHslow)./mtauHslow);

          vnoise = uniqueNoiseSTD*randn(1,ncells);
          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) mHslow];

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

          %% do numPreISIs then change I and do numPostISIs
          if any(spikes(:,tind))
            numISIs = numISIs + 1; % starts equal to -1 so second spike is first ISI
          end
          if numISIs<(numPreISIs+numDropISIs)
            I = inputMags(inputi);
          else
            I = inputMags2(inputi2);
          end
          if numISIs==(numPreISIs+numPostISIs+numDropISIs)
            break
          end

          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

        % calculate mean and variance of ISIs for each postsynaptic cell
        npost = 1; %size(post,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);
        ISIs = dt*diff(find(spikes(1,:)>0))/1000;

        if length(inputMags)==1 && length(inputMags2)>1 % add so we can take the average over runs
          figISIs(inputi2,:) = figISIs(inputi2,:) + ISIs((numDropISIs+1):end);
        elseif length(inputMags)>1 && length(inputMags2)==1
          figISIs(inputi,:) = figISIs(inputi,:) + ISIs((numDropISIs+1):end);
        elseif length(inputMags)>1 && length(inputMags2)>1
          % using a complex number to hold two values in one matrix:
          figISIs(inputi2,inputi) = figISIs(inputi2,inputi) + ISIs(numDropISIs+numPreISIs+1) + ISIs(end-1)*i;
        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%d\n', ...
            ncells,pcon,inputMags(inputi),inputMags2(inputi2),g,simdur/1000, useNoise*uniqueNoiseSTD,tau,postWeight,...
            mean(ISIs(2:end-numPostISIs-2)), ISIs(end-numPostISIs-2), ISIs(end-numPostISIs+1), ISIs(end-numPostISIs+2),ISIs(end),length(ISIs),...
            (inputMags2(inputi2)-inputMags(inputi)), ISIs(end-numPostISIs+1)-ISIs(end-numPostISIs+2),ISIs(end-numPostISIs+1)-ISIs(end));
          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%d\n', ...
          ncells,pcon,inputMags(inputi),inputMags2(inputi2),g,simdur/1000, useNoise*uniqueNoiseSTD,tau,postWeight,...
          mean(ISIs(2:end-numPostISIs-2)), ISIs(end-numPostISIs-2), ISIs(end-numPostISIs+1), ISIs(end-numPostISIs+2),ISIs(end),length(ISIs),...
          (inputMags2(inputi2)-inputMags(inputi)), ISIs(end-numPostISIs+1)-ISIs(end-numPostISIs+2),ISIs(end-numPostISIs+1)-ISIs(end));

        slopes(inputi,inputi2) = (ISIs(end-numPostISIs-1)-ISIs(end-numPostISIs))/(inputMags2(inputi2)-inputMags(inputi));
      end
    end

    % average ISIs over multiple runs if noisy network is in use:
    if length(inputMags)>1 || length(inputMags2)>1
      figISIs = figISIs./nruns;
    end
  end
end

state = state(:,1:tind);

if length(inputMags)==1 && length(inputMags2)==1
  start=260000;stop=360000;
  figure;
  plot((start:stop)*dt,state(1,start:stop)');
  set(gca,'xlim',[start*dt stop*dt]);
  set(gcf,'position',[82 404 257 185]);
  xlabel('Time (ms)','fontsize',12); ylabel('Voltage (mV)','fontsize',12);
  set(gca,'fontsize',12,'box','off')
  if slowslow
    if saveFigures && exist(['fig_Ackerslow_history_dependence_vtrace_0a.eps'])~=2
      print_eps(['fig_Ackerslow_history_dependence_vtrace_0a.eps'])
      saveas(gcf,['fig_Ackerslow_history_dependence_vtrace_0a.fig'])
    end
  else
    if saveFigures && exist(['fig_Acker_history_dependence_vtrace_0a.eps'])~=2
      print_eps(['fig_Acker_history_dependence_vtrace_0a.eps'])
      saveas(gcf,['fig_Acker_history_dependence_vtrace_0a.fig'])
    end
  end
  figure;
  plot((start:stop)*dt,state(2,start:stop)');
  set(gcf,'position',[82 404 257 227]);
  set(gca,'ygrid','on','xlim',[start*dt stop*dt]);
  xlabel('Time (ms)','fontsize',12); ylabel('Slow H','fontsize',12);
  set(gca,'fontsize',12,'box','off')
  if slowslow
    if saveFigures && exist(['fig_Ackerslow_history_dependence_htrace_0a.eps'])~=2
      print_eps(['fig_Ackerslow_history_dependence_htrace_0a.eps'])
      saveas(gcf,['fig_Ackerslow_history_dependence_htrace_0a.fig'])
    end
  else
    if saveFigures && exist(['fig_Acker_history_dependence_htrace_0a.eps'])~=2
      print_eps(['fig_Acker_history_dependence_htrace_0a.eps'])
      saveas(gcf,['fig_Acker_history_dependence_htrace_0a.fig'])
    end
  end
elseif (length(inputMags)==1 && length(inputMags2)>1) || (length(inputMags)>1 && length(inputMags2)==1)
  figure; plot(1./figISIs','b.-','Markersize',16);
  xlabel('Interspike interval #','fontsize',12); ylabel('1/ISI duration (Hz)','fontsize',12)
  set(gcf,'position',[82 404 257 302])
  set(gca,'ylim',[min(desiredFreqs)-0.1 max(desiredFreqs)+0.5]);
  set(gca,'fontsize',12,'box','off')
  if slowslow
    if saveFigures && exist(['fig_Ackerslow_history_dependence_multi_0a.eps'])~=2
      print_eps(['fig_Ackerslow_history_dependence_multi_0a.eps'])
      saveas(gcf,['fig_Ackerslow_history_dependence_multi_0a.fig'])
    end
  else
    if saveFigures && exist(['fig_Acker_history_dependence_multi_0a.eps'])~=2
      print_eps(['fig_Acker_history_dependence_multi_0a.eps'])
      saveas(gcf,['fig_Acker_history_dependence_multi_0a.fig'])
    end
  end
elseif length(inputMags)>1 && length(inputMags2)>1
  figure; imagesc(desiredFreqs,desiredFreqs,abs(real(figISIs)-imag(figISIs)));
  title('ISI overshoot magnitude (ms)','fontsize',12)
  xlabel('Starting frequency (Hz)','fontsize',12); ylabel('Ending frequency (Hz)','fontsize',12)
  set(gcf,'position',[82 404 322 302])
  set(gca,'fontsize',12,'box','off')
  if saveFigures && exist(['fig_RS' typePrefix '_histdep_ISI_colorplot.eps'])~=2
    print_eps(['fig_RS' typePrefix '_histdep_ISI_colorplot.eps'])
    saveas(gcf,['fig_RS' typePrefix '_histdep_ISI_colorplot.fig'])
  end

  figure; imagesc(desiredFreqs,desiredFreqs,abs(1./real(figISIs)-1./imag(figISIs)));
  title('Frequency overshoot magnitude (Hz)','fontsize',12)
  xlabel('Starting frequency (Hz)','fontsize',12); ylabel('Ending frequency (Hz)','fontsize',12)
  set(gcf,'position',[82 404 322 302])
  set(gca,'fontsize',12,'box','off')
  if saveFigures && exist(['fig_RS' typePrefix '_histdep_freq_colorplot.eps'])~=2
    print_eps(['fig_RS' typePrefix '_histdep_freq_colorplot.eps'])
    saveas(gcf,['fig_RS' typePrefix '_histdep_freq_colorplot.fig'])
  end
end

Loading data, please wait...