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
% 2d grid simulation
% eric zilli - october 13, 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).
%
% This script uses experimentally recorded rat trajectory (given by the
% hafting_trajectory function) as input to multiple VCOs which can
% interfere to produce spatial firing of a postsynaptic cell which,
% hopefully, is arrayed in the classical hexagonal grid pattern.
%
% Running this script requires that an FI curve of the network has been
% generated (e.g. by SI_acker_model_FI_relation.m) for a cell or network
% with identical parameters to those used here.
%
% This is generally way slow: behavioral timescales are big and neural
% timescales are tiny. Hints: using pcon=1 lets the script bypass the
% connectivity matrix which really speeds things up. If the resonant postsynaptic
% cell starts spiking too much the script will slow way down because the
% vector of spiketimes of the postsynaptic cell is a pre-allocated sparse
% vector. When adjusting postsynaptic cell parameters, start with
% short simulation to act as a sort of sanity check.
% For 2 VCO, 160 s simulations, my desktop takes about
% 16 hours. It is easily possible there are simple speedups that would make
% this more wieldy and the interested reader might learn a lot about the
% model by trying to find them!
%
% Interesting questions to continue this line of research:
% Can we analytically find appropriate parameters for the
% postsynaptic models given requirements such as: the field diameter should
% be about half the field spacing? (Here's a start: that means on the edge
% of a field between two neighboring fields,
% all the VCOs are at phase zero except for one which is at pi/2,
% corresponding to a fixed time given a baseline frequency, which means the
% decayed sum of the inputs from all but one VCO (decayed over the time corresponding
% to pi/2 radians, or one quarter the duration of the oscillator period)
% plus the input from one VCO should equal the firing threshold [but a tau
% too high will carry over activity from the previous cycle!]). If noise is
% correlated between the VCO inputs, what is the effect on the model and on
% the stability of the oscillations between the VCOs? Zilli et al 2009 only
% examined this question for abstract oscillators, not for these more
% complex models where the effects may be worse.
%
% More possibilities are allowed than the single set of preset parameters here,
% but postsynaptic parameters might need to be adjusted.
% Note that parameters will need to change if the baseline frequency
% changes.
%
% 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.
%
% This generates Figure 10 and has just that one "preset".

clear all;
% grab trajectory: this is a bit slow so best to do it once then comment it
% out

simdur = 240*1e3; % ms
trajdt = 0.1; % sampling rate for trajectory; ms
[dxs,dys,fdxs,fdys] = hafting_trajectory(simdur,trajdt);
clear dxs dys

dt = .01; % ms

useFilteredTrajectory = 1;

% only store 1 out of every storeskip timesteps
storeskip = 20;

% to extend a run, enter the total previous run duration and increase
% simdur above to the new end-of-simulation-time (and comment out the "clear
% all" line!)
extendPreviousRun = 0;
previousDur = 12000; % (ms)
if extendPreviousRun
  %% this may not work anymore! try on a really short simulation first if
  %% you want to use this
  'extending previous run instead of overwriting!'
end

nruns = 1;
nVCOs = 2;

runAbstract = 1;
runNetwork = 1;

% if ==1, abstract output will be plotted, ==0 means resonator output
% plotted in its place, ==-1 means nothing in those spaces
showAbstractFig = -1;

saveFigures = 0;

% live plot of network spatial activity
runningPlot = 0;

if runningPlot
  runh = figure;
end

useGapJunctions = 0;
slowslow = 0;
ncells = 250;
pcon = 1;
useNoise = 1;
commonNoiseScale = 0;
uniqueNoiseScale = 3.44*useNoise;
g = 80; % syn
I = -2.35;
spikeThreshold = 20; % mV

load Acker_sn_FI_n250.mat;

% prefix used for saving figures
simprefix = sprintf('filttraj_n%d_%gnoise_%ds_%s',ncells,useNoise,round(simdur/1000),datestr(now,'mmmmdd'));
if ~useFilteredTrajectory
  simprefix = ['un' simprefix];
end

% shared grid variables: freq = baselineFreq+beta*speed*(cos(prefHD-curHD));
baselineFreq = freqs(round(length(freqs)/2));
beta = 2; % Hz/(m/s)

baselineMult = 1;

% non-gated lif params:
tau = 40; % ms
membraneDecay = exp(-dt/tau);
posthWeight=0.15/ncells;
baseWeight = 0.8/ncells;

% gated lif params:
% I&F params
basegateDur = 40; % ms
tau = 50; % (msec)
gatedmembraneDecay = exp(-dt/tau);
weightMult = 1.2;
gatedposthWeight = weightMult/ncells/nVCOs;

% non-gated res params:
weightMult = 100;
resposthWeight=weightMult*1/ncells/(1+nVCOs);
resbaseWeight = 2*posthWeight;
posthc = -0.005;
posthw = baselineFreq*2*pi/1000; % Hz (times 2pi for angular freq, /1000 to convert to Hz

%% can downsample FI curve to test the needed level of precision
%   freqs = freqs(1:50:end);
%   currents = currents(1:50:end);

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)

% connectivity matrix
if pcon~=1
  C = (rand(ncells)<pcon)>0;
  C = sparse(C.*(1-eye(ncells))); % remove autoconnections
else
  C = 0;
end
stabilities = zeros(1,nruns);
stabilities2 = zeros(1,nruns);
for run=1
  % activity of posthsynaptic I&F
  npost = 3;
  if ~extendPreviousRun
    posth = zeros(npost, round(simdur/dt/storeskip)+1);
    post = zeros(npost,1);

    % spike times of posthsynaptic I&F
    posthspikes = spalloc(npost, round(simdur/dt/storeskip)+1, 6*simdur); % expect firing at 6 Hz

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

    % initial values
    v = -55*ones(ncells,nVCOs);
    mNa = zeros(ncells,nVCOs);
    hNa = zeros(ncells,nVCOs);
    mNap = zeros(ncells,nVCOs);
    nK = zeros(ncells,nVCOs);
    mHfast = zeros(ncells,nVCOs);
    mHslow = zeros(ncells,nVCOs);

    vb = -55*ones(ncells,1);
    mNab = zeros(ncells,1);
    hNab = zeros(ncells,1);
    mNapb = zeros(ncells,1);
    nKb = zeros(ncells,1);
    mHfastb = zeros(ncells,1);
    mHslowb = zeros(ncells,1);

    % save voltage of one cell over simulation:
    vhist = zeros(nVCOs, round(simdur/dt/storeskip)+1);
    vhistb = zeros(1, round(simdur/dt/storeskip)+1);
    % binary array indicating which cells fire on each time steps
    spikes = zeros(ncells, nVCOs);
    spikesb = zeros(ncells, 1);

    % not keeping full array anymore, but still want to know times any cell
    % spiked:
    clear VCOSpikeTimes;
    VCOSpikeTimes{nVCOs} = []; % implicitly create this struct/class/whatever thing
    BaseSpikeTimes = [];
    VCOinds = zeros(1,nVCOs);
    Baseind = 0;

    %% abstract grid model variables:
    % dbasePhi/dt = baselineFreq
    % dVCOPhi/dt = (baselineFreq+beta*speed*(cos(prefHD-curHD)))
    basePhi = 0; % baseline oscillator phase variable
    VCOPhi = zeros(1,nVCOs); % VCO phase variable(s)
    prefHDs = [0 2*pi/3]; % VCO preferred direction(s), (radians)
    abstractGrid = zeros(1, round(simdur/dt/storeskip)+1);
    basehist = zeros(1,round(simdur/dt));
    vcohist = zeros(nVCOs,round(simdur/dt));

    x = 0; % m
    dx = 0; % m
    y = 0; % m
    dy = 0; % m

    pos = zeros(2, round(simdur/dt/storeskip)+1);
    baseI = interp1(freqs,currents,baselineFreq,'linear');
    I = baseI;
    gatedI = baseI*ones(1,nVCOs);
    ISIs = zeros(1,nVCOs);

    speed = 0; % m/s
    gatedspeeds = zeros(1,nVCOs);

    % external path integration to hold I constant:
    XatLastSpike = zeros(1,nVCOs);
    YatLastSpike = zeros(1,nVCOs);

    basegateCount = 0;
  else
    posth(npost,round(simdur/dt)+1) = 0; % implicit extension
    abstractGrid(1,round(simdur/dt/storeskip)+1) = 0; % implicit extension
    vhist(nVCOs,round(simdur/dt/storeskip)+1) = 0; % implicit extension
    vhistb(nVCOs,round(simdur/dt/storeskip)+1) = 0; % implicit extension
    pos(2,round(simdur/dt/storeskip)+1) = 0; % implicit extension
  end

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

    if mod(t,round(simdur)/100)<=dt
      disp(sprintf('t = %d, %g elapsed',t,toc));
    end

    %%% filtered, interpolated trajectory (scale dx/dy by mismatch between
    %%% time steps to properly scale speed)
    if useFilteredTrajectory
      dx = fdxs(ceil(t/trajdt))*(dt/trajdt);
      dy = fdys(ceil(t/trajdt))*(dt/trajdt);
    else
      dx = dxs(ceil(t/trajdt))*(dt/trajdt);
      dy = dys(ceil(t/trajdt))*(dt/trajdt);
    end

    x = x+dx;
    y = y+dy;

    speed = abs(sqrt((dx^2 + dy^2)/(dt/1000)^2)); % (m/s)
    curHD = atan2(dy,dx);
    if mod(tind,storeskip)==1
      storeind=storeind+1;
      pos(:,storeind) = [x; y];
    end

    if runAbstract
      %% abstract grid model: (divide dt by 1000 because freqs are in Hz
      %% but dt is in ms)
      basePhi = basePhi + dt/1000*2*pi*baselineFreq;
      VCOAngularFreqs = 2*pi*(baselineFreq+beta*speed*(cos(prefHDs-curHD)));
      VCOPhi = VCOPhi + dt/1000*VCOAngularFreqs;
      basehist(tind) = basePhi;
      vcohist(:,tind) = VCOPhi;
    end
    if mod(tind,storeskip)==1
      abstractGrid(storeind) = nVCOs*cos(basePhi)+sum(cos(VCOPhi));
    end

    if runNetwork
      %% active VCOs
      for vco=1:nVCOs
        oldv = v(:,vco);
        alphamNa = -0.1*(oldv+23)./(exp(-0.1*(oldv+23))-1);
        betamNa = 4*exp(-(oldv+48)/18);
        mNa(:,vco) = mNa(:,vco) + dt*(alphamNa.*(1-mNa(:,vco)) - betamNa.*mNa(:,vco));

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

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

        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(:,vco) = mNap(:,vco) + dt*(alphamNap.*(1-mNap(:,vco)) - betamNap.*mNap(:,vco));

        minfHfast = 1./(1+exp((oldv+79.2)/9.78));
        mtauHfast = 0.51./(exp((oldv-1.7)/10) + exp(-(oldv+340)/52)) + 1;
        mHfast(:,vco) = mHfast(:,vco) + dt*((minfHfast-mHfast(:,vco))./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(:,vco) = mHslow(:,vco) + dt*((minfHslow-mHslow(:,vco))./mtauHslow);

        desiredFreq = baselineFreq+beta*speed*(cos(prefHDs(vco)-curHD));
        lowind = max(find(freqs<desiredFreq));
        highind = min(find(freqs>desiredFreq));
        proportion = (desiredFreq-freqs(lowind))/(freqs(highind)-freqs(lowind));
        I = currents(lowind) + proportion*(currents(highind)-currents(lowind));

        vnoise = uniqueNoiseScale*randn(ncells,1);
        if useGapJunctions
          if pcon==0
            v(:,vco) = v(:,vco) + dt*(vnoise - gH*(0.65*mHfast(:,vco) + 0.35*mHslow(:,vco)).*(v(:,vco)-EH) - (gNap*mNap(:,vco)+gNa*mNa(:,vco).^3.*hNa(:,vco)).*(v(:,vco)-ENa) - gK*nK(:,vco).^4.*(v(:,vco)-EK) - gLeak*(v(:,vco)-ELeak) + I)/Cm;
          elseif pcon==1
            v(:,vco) = v(:,vco) + dt*(vnoise - gH*(0.65*mHfast(:,vco) + 0.35*mHslow(:,vco)).*(v(:,vco)-EH) - (gNap*mNap(:,vco)+gNa*mNa(:,vco).^3.*hNa(:,vco)).*(v(:,vco)-ENa) - gK*nK(:,vco).^4.*(v(:,vco)-EK) - gLeak*(v(:,vco)-ELeak) + I + g*(sum(v(:,vco))-ncells*v(:,vco)))/Cm;
          else
            vdiffs = repmat(v(:,vco)',ncells,1) - repmat(v(:,vco),1,ncells);
            v(:,vco) = v(:,vco) + dt*(vnoise - gH*(0.65*mHfast(:,vco) + 0.35*mHslow(:,vco)).*(v(:,vco)-EH) - (gNap*mNap(:,vco)+gNa*mNa(:,vco).^3.*hNa(:,vco)).*(v(:,vco)-ENa) - gK*nK(:,vco).^4.*(v(:,vco)-EK) - gLeak*(v(:,vco)-ELeak) + I + g*sum(C*vdiffs,2))/Cm;
          end
        else
          if pcon==0
            v(:,vco) = v(:,vco) + dt*(vnoise - gH*(0.65*mHfast(:,vco) + 0.35*mHslow(:,vco)).*(v(:,vco)-EH) - (gNap*mNap(:,vco)+gNa*mNa(:,vco).^3.*hNa(:,vco)).*(v(:,vco)-ENa) - gK*nK(:,vco).^4.*(v(:,vco)-EK) - gLeak*(v(:,vco)-ELeak) + I)/Cm;
          elseif pcon==1
            v(:,vco) = v(:,vco) + dt*(vnoise - gH*(0.65*mHfast(:,vco) + 0.35*mHslow(:,vco)).*(v(:,vco)-EH) - (gNap*mNap(:,vco)+gNa*mNa(:,vco).^3.*hNa(:,vco)).*(v(:,vco)-ENa) - gK*nK(:,vco).^4.*(v(:,vco)-EK) - gLeak*(v(:,vco)-ELeak) + I + g*(sum(spikes(:,vco))-spikes(:,vco)')')/Cm;
          else
            v(:,vco) = v(:,vco) + dt*(vnoise - gH*(0.65*mHfast(:,vco) + 0.35*mHslow(:,vco)).*(v(:,vco)-EH) - (gNap*mNap(:,vco)+gNa*mNa(:,vco).^3.*hNa(:,vco)).*(v(:,vco)-ENa) - gK*nK(:,vco).^4.*(v(:,vco)-EK) - gLeak*(v(:,vco)-ELeak) + I + g*(C*spikes(:,vco))')/Cm;
          end
        end
        % spikes if upward crossing spikeThreshold mV
        spikes(:,vco) = (oldv<spikeThreshold).*(v(:,vco)>=spikeThreshold);

        if any(spikes(:,vco))
          % faster I hope:
          if mod(VCOinds(vco),1000)==0
            VCOSpikeTimes{vco}(length(VCOSpikeTimes{vco})+1000) = 0;
          end
          VCOinds(vco) = VCOinds(vco)+1;
          VCOSpikeTimes{vco}(VCOinds(vco)) = t;

          % now save current position for end of this cycle:
          XatLastSpike(vco) = x;
          YatLastSpike(vco) = y;
        end
      end

      if mod(tind,storeskip)==1
        vhist(:,storeind) = v(1,:);
        if any(isnan(v(1,:)))
          error('it is bad again')
        end
      end

      %% baseline oscillator
      oldvb = vb;
      alphamNa = -0.1*(oldvb+23)./(exp(-0.1*(oldvb+23))-1);
      betamNa = 4*exp(-(oldvb+48)/18);
      mNab = mNab + dt*(alphamNa.*(1-mNab) - betamNa.*mNab);

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

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

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

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

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

      vnoise = uniqueNoiseScale*randn(ncells,1);
      if useGapJunctions
        if pcon==0
          vb = vb + dt*(vnoise - gH*(0.65*mHfastb + 0.35*mHslowb).*(vb-EH) - (gNap*mNapb+gNa*mNab.^3.*hNab).*(vb-ENa) - gK*nKb.^4.*(vb-EK) - gLeak*(vb-ELeak) + baseI)/Cm;
        elseif pcon==1
          vb = vb + dt*(vnoise - gH*(0.65*mHfastb + 0.35*mHslowb).*(vb-EH) - (gNap*mNapb+gNa*mNab.^3.*hNab).*(vb-ENa) - gK*nKb.^4.*(vb-EK) - gLeak*(vb-ELeak) + baseI + g*(sum(vb)-ncells*vb))/Cm;
        else
          vdiffs = repmat(vb',ncells,1) - repmat(vb,1,ncells);
          vb = vb + dt*(vnoise - gH*(0.65*mHfastb + 0.35*mHslowb).*(vb-EH) - (gNap*mNapb+gNa*mNab.^3.*hNab).*(vb-ENa) - gK*nKb.^4.*(vb-EK) - gLeak*(vb-ELeak) + baseI + g*sum(C*vdiffs,2))/Cm;
        end
      else
        if pcon==0
          vb = vb + dt*(vnoise - gH*(0.65*mHfastb + 0.35*mHslowb).*(vb-EH) - (gNap*mNapb+gNa*mNab.^3.*hNab).*(vb-ENa) - gK*nKb.^4.*(vb-EK) - gLeak*(vb-ELeak) + baseI)/Cm;
        elseif pcon==1
          vb = vb + dt*(vnoise - gH*(0.65*mHfastb + 0.35*mHslowb).*(vb-EH) - (gNap*mNapb+gNa*mNab.^3.*hNab).*(vb-ENa) - gK*nKb.^4.*(vb-EK) - gLeak*(vb-ELeak) + baseI + g*(sum(spikesb,1)-spikesb')')/Cm;
        else
          vb = vb + dt*(vnoise - gH*(0.65*mHfastb + 0.35*mHslowb).*(vb-EH) - (gNap*mNapb+gNa*mNab.^3.*hNab).*(vb-ENa) - gK*nKb.^4.*(vb-EK) - gLeak*(vb-ELeak) + baseI + g*(C*spikesb)')/Cm;
        end
      end

      % spikes if upward crossing spikeThreshold mV
      spikesb = (oldvb<spikeThreshold).*(vb>=spikeThreshold);

      if any(spikesb)
        if mod(Baseind,1000)==0
          BaseSpikeTimes(length(BaseSpikeTimes)+1000) = 0;
        end
        Baseind = Baseind+1;
        BaseSpikeTimes(Baseind) = t;

        % if baseline spikes, active VCOs can influence posth for
        % basegateCount steps
        basegateCount = basegateDur/dt;
      else
        basegateCount = basegateCount-1;
      end


      if mod(tind,storeskip)==1
        vhistb(:,storeind) = vb(1);
      end

      % while basegateCount>0, the gate is open
      if basegateCount>0
        basegate = 1;
      else
        basegate = 0;
      end

      spikesum = sum(sum(spikes));
      spikebsum = sum(spikesb);
      % posth cell 1 = non-gated lif
      post(1) = post(1)*membraneDecay + posthWeight*spikesum+baselineMult*baseWeight*spikebsum;
      if post(1)>1
        posthspikes(1,storeind) = posthspikes(1,storeind) + post(1)>1;
        post(1) = 0;
      end
      % posth cell 2 = gated lif
      post(2) = post(2)*gatedmembraneDecay + basegate*gatedposthWeight*spikesum;
      if post(2)>1
        posthspikes(2,storeind) = posthspikes(2,storeind) + post(2)>1;
        post(2) = 0;
      end
      % posth cell 3 = non-gated resonant
      post(3) = post(3) + dt*((posthc+posthw*i)*post(3) + resposthWeight*spikesum+baselineMult*resbaseWeight*spikebsum);
      %       if mod(tind,storeskip)==1
      %% since this doesn't spike a lot, let's add up the spikes instead
      %% of periodically sampling whether the cell is spiking
      posthspikes(3,storeind) = posthspikes(3,storeind) + real(post(3))>1;
      %       end
      if real(post(3))>1
        post(3) = i;
      end

      if mod(tind,storeskip)==1
        posth(:,storeind) = post;
      end

      if runningPlot
        figure(runh);
        plot(pos(1,:),pos(2,:),'.',pos(1,find(posthspikes>0)),pos(2,find(posthspikes>0)),'R.'), title('network VCO spatial firing'); xlabel('position (m)'), ylabel('position (m)');
        set(gca,'xlim',[-1 1],'ylim',[-1 1]);
        drawnow
      end
    end
  end
  toc

  for vco=1:nVCOs
    VCOSpikeTimes{vco}(VCOSpikeTimes{vco}==0) = [];
  end
  BaseSpikeTimes(BaseSpikeTimes==0) = [];

  clear C

  if runNetwork
    for ce=1:npost
      spikeTimes = find(posthspikes(ce,:)>0);

      nspikes(ce) = length(spikeTimes);

      ISIs = dt*diff(spikeTimes)/1000;

      % 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
  end
end

for vco=1:nVCOs
  VCOSpikeTimes{vco}(VCOSpikeTimes{vco}==0) = [];
  if length(VCOSpikeTimes{vco})>40000
    VCOSpikeTimes{vco} = VCOSpikeTimes{vco}(1:10:end);
  end
end
BaseSpikeTimes(BaseSpikeTimes==0) = [];
if length(BaseSpikeTimes)>40000
  BaseSpikeTimes = BaseSpikeTimes(1:10:end);
end

if simdur>4000
  plotduration = 4000; % ms
else
  plotduration  = simdur;
end
figure; plot(vhist(:,1:(plotduration/dt/storeskip))');
hold on; plot(vhistb(1:(plotduration/dt/storeskip))','r');
set(gcf,'position',[82 404 3*257 2/3*302])
set(gca,'fontsize',12,'box','off')
xlabel('Time (ms)','fontsize',12); ylabel('Voltage (mV)','fontsize',12)
if saveFigures && exist(['fig_Acker' typePrefix '_' simprefix '_2D_vcos_trace_0a.eps'])~=2
  print_eps(['fig_Acker' typePrefix '_' simprefix '_2D_vcos_trace_0a.eps'])
  saveas(gcf,['fig_Acker' typePrefix '_' simprefix '_2D_vcos_trace_0a.fig'])
  close
end

% plot non-gated lif posth:
figure; plot(dt:dt:plotduration/storeskip,posth(1,1:(plotduration/dt/storeskip)));
if any(posthspikes(1,1:(plotduration/dt/storeskip))>0)
  hold on; plot(dt*find(posthspikes(1,1:(plotduration/dt/storeskip))>0),1,'r.','MarkerSize',16);
end
set(gcf,'position',[82 404 3*257 2/3*302])
set(gca,'fontsize',12,'box','off')
xlabel('Time (ms)','fontsize',12); ylabel('Voltage (mV)','fontsize',12)
set(gca,'ylim',[-0.1 1.1]);
if saveFigures && exist(['fig_Acker' typePrefix '_' simprefix '_2D_posth_nogatelif_trace_0a.eps'])~=2
  print_eps(['fig_Acker' typePrefix '_' simprefix '_2D_posth_nogatelif_trace_0a.eps'])
  saveas(gcf,['fig_Acker' typePrefix '_' simprefix '_2D_posth_nogatelif_trace_0a.fig'])
end

% plot gated lif posth:
figure; plot(dt:dt:plotduration/storeskip,posth(2,1:(plotduration/dt/storeskip)));
if any(posthspikes(2,1:(plotduration/dt/storeskip))>0)
  hold on; plot(dt*find(posthspikes(2,1:(plotduration/dt/storeskip))>0),1,'r.','MarkerSize',16);
end
set(gcf,'position',[82 404 3*257 2/3*302])
set(gca,'fontsize',12,'box','off')
xlabel('Time (ms)','fontsize',12); ylabel('Voltage (mV)','fontsize',12)
set(gca,'ylim',[-0.1 1.1]);
if saveFigures && exist(['fig_Acker' typePrefix '_' simprefix '_2D_posth_gating_trace_0a.eps'])~=2
  print_eps(['fig_Acker' typePrefix '_' simprefix '_2D_posth_gating_trace_0a.eps'])
  saveas(gcf,['fig_Acker' typePrefix '_' simprefix '_2D_posth_gating_trace_0a.fig'])
end

% plot non-gated res posth:
figure; plot(dt:dt:plotduration/storeskip,real(posth(3,1:(plotduration/dt/storeskip))));
if any(posthspikes(3,:))
  hold on; plot(dt*find(posthspikes(3,1:(plotduration/dt/storeskip))>0),1,'r.','MarkerSize',16);
end
set(gcf,'position',[82 404 3*257 2/3*302])
set(gca,'fontsize',12,'box','off')
xlabel('Time (ms)','fontsize',12); ylabel('Voltage (mV)','fontsize',12)
set(gca,'ylim',[-1.1 1.1]);
if saveFigures && exist(['fig_Acker' typePrefix '_' simprefix '_2D_posth_nogateres_trace_0a.eps'])~=2
  print_eps(['fig_Acker' typePrefix '_' simprefix '_2D_posth_nogateres_trace_0a.eps'])
  saveas(gcf,['fig_Acker' typePrefix '_' simprefix '_2D_posth_nogateres_trace_0a.fig'])
end

% these in normalized coords:
trajH = 3/4*1/2; % <1/2
trajW = 3/4*1/3; % <1/3
diffH = 2/3*1/3; % <1/3
diffW = trajW; % <1/3

col1L = 1/4*1/3; % <1/3-trajW
col2L = 1/3+1/5*1/3;
col3L = 2/3+1/5*1/3;
row2B = 1/6*1/2; % <1/3-trajH
diffL = 1/5*1/3;
diffB1 = .007+1/8*1/2;
diffB2 = diffB1+0.005;
%   diffB3 = trajB;

figure('position',[82 404 3*257 2*277])
if nVCOs==2
  % gated lif:
  subplot('position',[col1L+0 row2B+1/2 trajW trajH]);
  plot(pos(1,1:500:length(pos)),pos(2,1:500:length(pos)),'.','MarkerSize',8);
  hold on; plot(pos(1,round(find(posthspikes(2,:)>0))),pos(2,round(find(posthspikes(2,:)>0))),'R.','MarkerSize',15)
  title('Gated posthsynaptic cell','fontsize',12)
  xlabel('Position (m)','fontsize',12), ylabel('Position (m)','fontsize',12)
  set(gca,'fontsize',12,'box','off')
  set(gca,'xlim',[min(pos(1,:)) max(pos(1,:))],'ylim',[min(pos(2,:)) max(pos(2,:))])

  % gated lif xcorr:
  edges{1} = linspace(min(pos(1,:)),max(pos(1,:)),50);
  edges{2} = linspace(min(pos(2,:)),max(pos(2,:)),50);
  rate = hist3([pos(1,round(find(posthspikes(2,:)>0))); pos(2,round(find(posthspikes(2,:)>0)))]','Edges',edges);
  subplot('position',[col2L+0 row2B+1/2 trajW trajH]); imagesc(rot90(conv2(rate,rate,'same')))
  if showAbstractFig
    title('Network model autocorrelogram','fontsize',12)
  else
    title(['Spatial autocorrelogram'],'fontsize',12)
  end
  set(gca,'xtick',[],'ytick',[]);
  set(gca,'fontsize',12)
  axis tight

  if showAbstractFig==1
    % abstract:
    abstractThr = 1.8;
    subplot('position',[col1L+0 row2B+0 trajW trajH]);
    plot(pos(1,1:500:length(pos)),pos(2,1:500:length(pos)),'.','MarkerSize',8);
    hold on; plot(pos(1,(find(abstractGrid>abstractThr))),pos(2,(find(abstractGrid>abstractThr))),'R.','MarkerSize',15)
    title('Abstract model','fontsize',12)
    xlabel('Position (m)','fontsize',12), ylabel('Position (m)','fontsize',12)
    set(gca,'xlim',[min(pos(1,:)) max(pos(1,:))],'ylim',[min(pos(2,:)) max(pos(2,:))])
    set(gca,'fontsize',12,'box','off')
    
    rate = hist3([pos(1,find(abstractGrid>abstractThr)); pos(2,find(abstractGrid>abstractThr))]','Edges',edges);
    subplot('position',[col2L+0 row2B trajW trajH]); imagesc(rot90(conv2(rate,rate,'same')))
    title(['Abstract model autocorrelogram'],'fontsize',12)
    set(gca,'xtick',[],'ytick',[]);
    set(gca,'fontsize',12)
    axis tight
  elseif showAbstractFig==0
    % non-gated res:
    subplot('position',[col1L+0 row2B+0 trajW trajH]);
    plot(pos(1,1:500:length(pos)),pos(2,1:500:length(pos)),'.','MarkerSize',8);
    hold on; plot(pos(1,round(find(posthspikes(3,:)>0))),pos(2,round(find(posthspikes(3,:)>0))),'R.','MarkerSize',15)
    title('Resonator posthsynaptic cell','fontsize',12)
    xlabel('Position (m)','fontsize',12), ylabel('Position (m)','fontsize',12)
    set(gca,'xlim',[min(pos(1,:)) max(pos(1,:))],'ylim',[min(pos(2,:)) max(pos(2,:))])
    set(gca,'fontsize',12,'box','off')

    rate = hist3([pos(1,round(find(posthspikes(3,:)==1))); pos(2,round(find(posthspikes(3,:)==1)))]','Edges',edges);
    subplot('position',[col2L+0 row2B trajW trajH]); imagesc(rot90(conv2(rate,rate,'same')))
    title(['Resonator autocorrelogram'],'fontsize',12)
    set(gca,'xtick',[],'ytick',[]);
    set(gca,'fontsize',12)
    axis tight
  end

  % vco 1 phase diff:
  spikeskips = 25;
  subplot('position',[col3L diffB2+2/3 diffW diffH]);
  plot(mod(vcohist(1,round(VCOSpikeTimes{1}(1:spikeskips:end)/dt)),2*pi),'.')
  %   hold on; plot(2-1e4*dxs(round(VCOSpikeTimes{1}(1:spikeskips:end)/dt)),'k.');
  title('VCO 1 phase drift','fontsize',12)
  ylabel('Phase difference (rad)','fontsize',12)
  set(gca,'xlim',[0 length(round(VCOSpikeTimes{1}(1:spikeskips:end)))],'ylim',[0 2*pi],'xtick',[])
  set(gca,'fontsize',12,'box','off')
  
  subplot('position',[diffL+2/3 diffB2+1/3 diffW diffH]);
  plot(mod(vcohist(2,round(VCOSpikeTimes{2}(1:spikeskips:end)/dt)),2*pi),'.')
  title('VCO 2 phase drift','fontsize',12)
  ylabel('Phase difference (rad)','fontsize',12)
  set(gca,'xlim',[0 length(round(VCOSpikeTimes{2}(1:spikeskips:end)))],'ylim',[0 2*pi],'xtick',[])
  set(gca,'fontsize',12,'box','off')

  % base phase diff:
  subplot('position',[col3L diffB2 diffW diffH]);
  plot(mod(basehist(round(BaseSpikeTimes(1:spikeskips:end)/dt)),2*pi),'.')
  title('Baseline phase drift','fontsize',12)
  xlabel('Time','fontsize',12)
  ylabel('Phase difference (rad)','fontsize',12)
  set(gca,'xlim',[0 length(round(BaseSpikeTimes(1:spikeskips:end)))],'ylim',[0 2*pi],'xtick',[])
  set(gca,'fontsize',12,'box','off')
end
if saveFigures && exist(['fig_Acker' typePrefix '_' simprefix '_2D_all_0a.eps'])~=2
  save(['fig_Acker' typePrefix '_' simprefix '_2D_variables_0a.mat'])
  % free some memory:
  %   clear pos abstractGrid vhist vcohist basehist posth
  print_eps(['fig_Acker' typePrefix '_' simprefix '_2D_all_0a.eps'])
  saveas(gcf,['fig_Acker' typePrefix '_' simprefix '_2D_all_0a.fig'])
  close
end

Loading data, please wait...