% 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