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
                            
% 2d grid simulation using izhikevich simple model neuron network
% eric zilli - july 2, 2009 using simple model MATLAB code from
% Izhikevich's book modified to be a network
%
% release version 1.0. check modeldb for updates.
%
% this source code is released into the public domain
%
% This script 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_simple_model_FI_relation.m) for a cell or network
% with identical parameters to those used here.
%
% This is generally 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, 320 s simulations, my slow laptop takes about
% 1 hour. 160 s simulations at 30 minutes wall time is a good place to
% tweak parameters. 4 s simulations at 1.5 min to ballpark params.
%
% Presets are given for figures in the paper which all Type 2 excitable,
% but Type 1 excitability is as easy as setting the variable and providing
% the proper FI curve. 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.

% clear all

simdur = 1*240e3; 320e3; % ms
trajdt = 0.1; % sampling rate for returned trajectory (will interp down to dt later); ms
if simdur>4e3
  [dxs,dys,fdxs,fdys] = hafting_trajectory(simdur,trajdt);
else
  [dxs,dys,fdxs,fdys] = hafting_trajectory(4e3,trajdt);
end

% % straight line trajectory:
% fdxs = 4e-5*ones(size(fdxs));
% fdys = 0*fdys;

% % bent line trajectory:
% speedmult = 0.1;
% fdxs(1:end/2) = speedmult*.0005*dt;
% fdys(1:end/2) = 0;
% fdxs(end/2:end) = -speedmult*-.0005*dt;
% fdys(end/2:end) = 0.00001*speedmult*.0004*dt;

% % elliptical trajectory:
% tv=(0:dt:simdur)*2*2*pi/8e3; a=0.001/4; b=0.00025*1.5;
% fdxs = -dt*(- a*cos(pi/2)*sin(tv) - b*cos(tv)*sin(pi/2));
% fdys = dt*(- b*cos(pi/2)*cos(tv) - a*sin(pi/2)*sin(tv));

% if you're trying to find good postsynaptic cell parameters, you're gonna
% keep some sanity by seeding the random number generator to compare runs
% with different params:
rand('seed',7);
randn('seed',7);

% type=1 (Figures S5, S9) 1 VCO, Class 2 excitable, n=1, noise=0, filt
% type=2 (Figure 3) 2 VCOs, Class 2 excitable, n=1, noise=0, filt
% type=3 (Figure S1) 2 VCOs, Class 2 excitable, n=1, noise=0, unfilt
% type=4 (Figure 4) 2 VCOs, Class 2 excitable, n=1, noise=full, filt
% type=5 (Figures 6, 7) 2 VCOs, Class 2 excitable, synaptic, n=250, noise=full, filt
% type=6 (Figures S6, S7) 2 VCOs, Class 2 excitable, gap-junction, n=250, noise=full, filt (the preset params for this one could use work!)
% type=7 3 VCOs, Class 2 excitable, synaptic, n=250, noise=full, filt
% type=8 3 VCOs, Class 2 excitable, gap-junction, n=250, noise=full, filt (doesn't work well)
% type=9 (Figure 9) 2 VCOs, Class 2 excitable inhibitory, gap-junction, n=250, noise=full, filt
% type=10 2 VCOs, Class 2 excitable, gap-junction, n=250, noise=full, filt, nPostIn=1
% type=11 2 VCOs, Class 2 excitable, synaptic, n=5000, p=0.01
% type=12 6 VCOs, Class 2 excitable, synaptic, n=250, noise=full, filt, cheatprecession
% type=13 6 VCOs, Class 2 excitable, n=1, noise=0, filt, cheatprecession
% type=14 6 VCOs, Class 2 excitable, gap-junction, n=250, noise=full, filt, nPostIn=25, cheatprecession (doesn't work, ninhib messes up FI curve)
% type=15 sandbox
for type=[7]
  dt = .1; % ms

  nruns = 1;
  nVCOs = 2; % might be overridden below

  % generally you will want to run both models so that you have the
  % respective VCO phase differences as a visual performance measure (in
  % addition to the spatial firing itself)
  runAbstract = 1;
  runNetwork = 1;

  % generally leave these = 1, unless debugging one or the other
  runNetworkVCOs = 1;
  runNetworkBaseline = 1;

  % if true, figures will be saved to disk if a figure of the same name does
  % not already exist
  saveFigures = 0;

  % plotType = 0 means gated LIF post on top, first quarter time network or
  % abstract on bottom (Figures 3, 4, S1)
  % plotType = 1 means gated LIF post on top, non-gated res on bottom
  % (Figures 6, 7, S6, S7)
  % plotType = 2 means all three post plus abstract spatial firing for
  % nVCOs=1 (overrides abstractThr) (Figure S5, S9)
  % plotType = 3 means all three post plus xcorrs plus abstract (useful for setting params)
  plotType = 3;

  % if true, plots traces of cell #1 from each VCO and the postsynaptic
  % cells (each one gets its own figure)
  % this makes Figures 7, 12, 14)
  plotVCOsAndPosts = 1;

  % if journalChargesALotForColor, gray-scale plots will be produced
  journalChargesALotForColor=1;

  % if 1, plots the velocity component on the VCO 1 phase diff plot for
  % plotType==0
  plotVelocity = 0;

  % for plotType = 0, this controls whether the bottom plot is the abstract
  % model's output or the first 1/4 of the simulation for the network model
  showAbstractFig = 0;

  % if =1 the postsynaptic cells (LIF) are quiescent at rest (I=0), if =0 the
  % first postsynaptic cell (simple model) actively (tonically) spikes at
  % I=0 and its inputs become inhibitory
  stablePost = 1;

  %% baseline synapses onto the post cell are this many times stronger than
  %% VCO syns (not counting that they are also nVCO times stronger)
  baselineMult = 1;

  % live plot of network spatial activity
  runningPlot = 0;

  % number of cells in each VCO the postsynaptic cell receives input from
  nPostIn = 1;
  
  % if cheatHDRectification=1, VCO inputs to the grid cell are
  % suppressed if the current head direction is not within 180 degrees of
  % the VCO's preferred direction
  cheatHDRectification=0;

  % if opponentVCOInhibition=1, cells 1:nInhib of each VCO get an input
  % current of magnitude opponentInhibAmp if the animal's direction is not
  % within 90 degrees of its preferred direction
  opponentVCOInhibition = 0;
  nInhib = 0;
  opponentInhibAmp = 0;
  
  % if non-empty, we'll save the VCO spike times to disk (for efficiently
  % testing postsynaptic parameters, can just generate the VCOs themselves
  % one and re-use the spikes)
    saveVCOSpikesFileName = '';
%   saveVCOSpikesFileName = '15s_res2_fi_spikes_dur4s.txt';

  % if non-empty, we'll load VCO spike times from disk from this file
  loadVCOSpikesFileName = '';
%     loadVCOSpikesFileName = 'testspikes.txt';

  % connectivity matrix
  % if not-empty, this filed will be loaded to provide the connectivity
  % matrix C
  loadC = '';
  %   % to make a C:
  %   C = (rand(ncells)<pcon)>0;
  %   C = sparse(C-eye(size(C))>0); % remove autoconnections
  %   save C_simple_nX_pX.mat C

  if runningPlot
    runh = figure;
  end

  %%% ANY OF THESE may be replaced in the type-unique parameters below

  % if baselineLIFGating = 1, the LIF only gets active VCO inputs and only during the
  % basegateDur milliseconds after each baseline VCO spike
  baselineLIFGating = 1;
  basegateDur = 10; % ms
  abstractThr=3;
  simprefix = '';
  ncells = 1;
  pcon = 1;
  useNoise = 0;
  commonNoiseSTD = 0;
  useFilteredTrajectory = 1;

  if type==1
    simprefix = sprintf('filthaftingtraj_n1_0noise_1vco_%s',datestr(now,'mmmmdd'));
    excitationClass = 2;
    useGapJunctions = 1;
    showAbstractFig = 1;
    plotType = 2;
    ncells = 1;
    nPostIn = ncells;
    pcon = 1;
    useNoise = 0;
    useFilteredTrajectory = 0;
    load simple_model_RS2_FI_Jan09_n1.mat;
    uniqueNoiseSTD = 100*useNoise;
    g = 0;
    nVCOs = 1;
    abstractThr=1.65*nVCOs;
    % shared grid variables: freq = baselineFreq+beta*speed*(cos(prefHD-curHD));
    baselineFreq = freqs(-1+round(length(freqs)/2)); % Hz
    beta = 2; % Hz/(m/s)

    % non-gated lif params:
    tau = 5; % ms
    membraneDecay = exp(-dt/tau);
    postWeight=0.91;
    baseWeight = 1*postWeight;

    % gated lif params:
    basegateDur = 25; % ms
    tau = 20; % (msec)
    gatedmembraneDecay = exp(-dt/tau);
    weightMult = 1.3;
    gatedpostWeight = weightMult/ncells/nVCOs;

    % non-gated res params:
    weightMult = 1;
    respostWeight=2;
    resbaseWeight = 6;
    postc = -0.01;
    postw = baselineFreq*2*pi/1000; % Hz

  elseif type==2
    simprefix = sprintf('filthaftingtraj_n1_0noise_01res%s',datestr(now,'mmmmdd'));
    excitationClass = 2;
    useGapJunctions = 0;
    showAbstractFig = 1;
    plotType = 0;
    plotVelocity = 1;
    ncells = 1;
    nPostIn = ncells;
    pcon = 0;
    useNoise = 0;
    uniqueNoiseSTD = 100*useNoise;
    useFilteredTrajectory = 1;
    load simple_model_RS2_FI_Jan09_n1.mat;
    g = 0;
    % shared grid variables: freq = baselineFreq+beta*speed*(cos(prefHD-curHD));
    baselineFreq = freqs(round(length(freqs)/2)); % Hz
    beta = 2; % Hz/(m/s)

    % non-gated lif params:
    tau = 40*7.8989/baselineFreq; 35; % ms
    membraneDecay = exp(-dt/tau);
    postWeight = 0.14;
    baseWeight = .8;

    % gated lif params:
    basegateDur = 40; % ms
    tau = 25; % (msec)
    gatedmembraneDecay = exp(-dt/tau);
    weightMult = 1.4;
    gatedpostWeight = weightMult/ncells/nVCOs;

    % non-gated res params:
    weightMult = 1;
    respostWeight=1.4; weightMult*1/ncells/(nVCOs);
    resbaseWeight = 6; 5; 2*postWeight;
    postc = -0.01;
    postw = baselineFreq*2*pi/1000; % Hz (times 2pi for angular freq, /1000 to convert to /s
  elseif type==3
    simprefix = sprintf('unfilthaftingtraj_n1_0noise_01res_%s',datestr(now,'mmmmdd'));
    excitationClass = 2;
    useGapJunctions = 0;
    showAbstractFig = 0;
    plotType = 0;
    ncells = 1;
    nPostIn = ncells;
    pcon = 0;
    useNoise = 0;
    useFilteredTrajectory = 0;
    load simple_model_RS2_FI_Jan09_n1.mat;
    uniqueNoiseSTD = 100*useNoise;
    g = 0;
    % shared grid variables: freq = baselineFreq+beta*speed*(cos(prefHD-curHD));
    baselineFreq = freqs(round(length(freqs)/2)); % Hz
    beta = 2; % Hz/(m/s)
    % non-gated lif params:
    tau = 35; % ms
    membraneDecay = exp(-dt/tau);
    postWeight = 0.16;
    baseWeight = .8;

    % gated lif params:
    basegateDur = 40; % ms
    tau = 25; % (msec)
    gatedmembraneDecay = exp(-dt/tau);
    weightMult = 1.4;
    gatedpostWeight = weightMult/ncells/nVCOs;

    % non-gated res params:
    weightMult = 1;
    respostWeight=1.5; weightMult*1/ncells/(nVCOs);
    resbaseWeight = 5; 2*postWeight;
    postc = -0.0075;
    postw = baselineFreq*2*pi/1000; % Hz
  elseif type==4
    simprefix = sprintf('filthaftingtraj_n1_1noise_01res_%s',datestr(now,'mmmmdd'));
    excitationClass = 2;
    useGapJunctions = 0;
    showAbstractFig = 0;
    plotType = 0;
    ncells = 1;
    nPostIn = ncells;
    pcon = 0;
    useNoise = 1;
    useFilteredTrajectory = 1;
    basegateDur = 40;
    load simple_model_RS2n_FI_Jan09_n1;
    uniqueNoiseSTD = 100*useNoise;
    % shared grid variables: freq = baselineFreq+beta*speed*(cos(prefHD-curHD));
    baselineFreq = freqs(round(length(freqs)/2)); % Hz
    beta = 2; % Hz/(m/s)
    %% Use whatever parameters you like here, you won't get a grid cell!
    % non-gated lif params:
    tau = 35; % ms
    membraneDecay = exp(-dt/tau);
    postWeight = 0.16;
    baseWeight = .8;

    % gated lif params:
    basegateDur = 40; % ms
    tau = 25; % (msec)
    gatedmembraneDecay = exp(-dt/tau);
    weightMult = 1.4;
    gatedpostWeight = weightMult/ncells/nVCOs;

    % non-gated res params:
    weightMult = 1;
    respostWeight=1.5; weightMult*1/ncells/(nVCOs);
    resbaseWeight = 5; 2*postWeight;
    postc = -0.0075;
    postw = baselineFreq*2*pi/1000; % Hz
  elseif type==5
    simprefix = sprintf('filthaftingtraj_n250_1noise_01res_%s',datestr(now,'mmmmdd'));
    excitationClass = 2;
    useGapJunctions = 0;
    showAbstractFig = 0;
    plotType = 1;
    ncells = 250;
    nPostIn = ncells;
    pcon = 1;
    useNoise = 1;
    useFilteredTrajectory = 1;
    basegateDur = 1;
    load simple_model_RS2sn_FI_Jan09_n250.mat;
    uniqueNoiseSTD = 100*useNoise;
    g = 256;
    baselineFreq = freqs(round(2+length(freqs)/2)); % Hz
    beta = 2; % Hz/(m/s)

    % non-gated lif params:
    tau = 5; % ms
    weightMult = 0.7;
    membraneDecay = exp(-dt/tau);
    postWeight=weightMult*1/ncells/(2+nVCOs);
    baseWeight = 2*postWeight;

    % gated lif params:
    basegateDur = 1; 50; % ms
    tau = 5; % (msec)
    gatedmembraneDecay = exp(-dt/tau);
    weightMult = 1.5; 1.3;
    gatedpostWeight = 0.0016;

    % non-gated res params:
    weightMult = 1;
    respostWeight=0.002;
    resbaseWeight = 0.0024;
    postc = -0.01;
    postw = baselineFreq*2*pi/1000; % Hz
  elseif type==6
    simprefix = sprintf('filthaftingtraj_n250_1noise_02res_%s',datestr(now,'mmmmdd'));
    excitationClass = 2;
    useGapJunctions = 1;
    showAbstractFig = 0;
    plotType = 1;
    ncells = 250;
    nPostIn = ncells;
    pcon = 1;
    useNoise = 1;
    useFilteredTrajectory = 1;
    load simple_model_RS2gn_FI_Jan09_n250.mat;
    uniqueNoiseSTD = 100*useNoise;
    g = 0.1;

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

    % non-gated lif params:
    tau = 20; % ms
    membraneDecay = exp(-dt/tau);
    postWeight=.19/250;
    baseWeight = .8/250;

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

    % non-gated res params:
    weightMult = 8;
    respostWeight=weightMult*1/ncells/(nVCOs);
    resbaseWeight = 16*postWeight;
    postc = -0.01;
    postw = baselineFreq*2*pi/1000; % Hz
  elseif type==7
    rand('seed',1);
    randn('seed',1);
    simprefix = sprintf('filthaftingtraj_n250_1noise_3vco_%s',datestr(now,'mmmmdd'));
    excitationClass = 2;
    useGapJunctions = 0;
    showAbstractFig = 0;
    plotType = 2;
    nVCOs = 3;
    ncells = 250;
    nPostIn = ncells;
    pcon = 1;
    useNoise = 1;
    useFilteredTrajectory = 1;
    load simple_model_RS2sn_FI_Jan09_n250.mat;
    uniqueNoiseSTD = 100*useNoise;
    g = 256;
    baselineFreq = freqs(round(-1+length(freqs)/2)); % Hz
    beta = 2; % Hz/(m/s)

    % non-gated lif params:
    tau = 5; % ms
    weightMult = 1; 0.8;
    membraneDecay = exp(-dt/tau);
    postWeight=weightMult*1/ncells/(2+nVCOs);
    baseWeight = 1.5*postWeight; 2*postWeight;

    % gated lif params:
    basegateDur = .05; 1; 50; % ms
    tau = 5; % (msec)
    gatedmembraneDecay = exp(-dt/tau);
    gatedpostWeight = 0.0014;

    % non-gated res params:
    weightMult = 7;
    respostWeight=0.0015; 0.002;
    resbaseWeight = 0.003; 0.0024;
    postc = -0.01;
    postw = baselineFreq*2*pi/1000; % Hz
  elseif type==8
    simprefix = sprintf('filthaftingtraj_n250_1noise_3vco_%s',datestr(now,'mmmmdd'));
    excitationClass = 2;
    useGapJunctions = 1;
    showAbstractFig = 0;
    plotType = 3;
    nVCOs = 3;
    ncells = 250;
    nPostIn = ncells;
    pcon = 1;
    useNoise = 1;
    useFilteredTrajectory = 1;
    load simple_model_RS2gn_FI_Jan09_n250.mat;
    uniqueNoiseSTD = 100*useNoise;
    g = 0.1;
    baselineFreq = freqs(1+round(length(freqs)/2)); % Hz
    beta = 2; % Hz/(m/s)

    % non-gated lif params:
    tau = 30; % ms
    membraneDecay = exp(-dt/tau);
    postWeight=.15/250;
    baseWeight = .8/250;

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

    % non-gated res params:
    weightMult = 9;
    respostWeight=weightMult*1/ncells/(nVCOs);
    resbaseWeight = 16*postWeight;
    postc = -0.01;
    postw = baselineFreq*2*pi/1000; % Hz
  elseif type==9
    simprefix = sprintf('filthaftingtraj_n250_inhib_1noise_02res_%s',datestr(now,'mmmmdd'));
    excitationClass = 2;
    useGapJunctions = 1;
    showAbstractFig = 0;
    plotType = 3;
    ncells = 250;
    nPostIn = ncells;
    stablePost = 0;
    pcon = 1;
    useNoise = 1;
    useFilteredTrajectory = 1;
    load simple_model_RS2gn_FI_Jan09_n250.mat;
    uniqueNoiseSTD = 100*useNoise;
    g = 0.1;

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

    % simple model postsynaptic params:
    tau = 20; % ms
    membraneDecay = exp(-dt/tau);
    postWeight=1000*-.19/250;
    baseWeight = 1000*-.8/250;

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

    % non-gated res params:
    weightMult = 8;
    respostWeight=weightMult*1/ncells/(nVCOs);
    resbaseWeight = 16*postWeight;
    postc = -0.01;
    postw = baselineFreq*2*pi/1000; % Hz
  elseif type==10
    simprefix = sprintf('filthaftingtraj_n250_1noise_02res_%s',datestr(now,'mmmmdd'));
    excitationClass = 2;
    useGapJunctions = 1;
    showAbstractFig = 0;
    plotType = 1;
    ncells = 250;
    nPostIn = 1;
    pcon = 1;
    useNoise = 1;
    useFilteredTrajectory = 1;
    load simple_model_RS2gn_FI_Jan09_n250.mat;
    uniqueNoiseSTD = 100*useNoise;
    g = 0.1;

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

    % non-gated lif params:
    tau = 40; % ms
    membraneDecay = exp(-dt/tau);
    postWeight = 0.14;
    baseWeight = .8;

    % gated lif params:
    basegateDur = 30; % ms
    tau = 20; % (msec)
    gatedmembraneDecay = exp(-dt/tau);
    %     weightMult = 1.4;
    gatedpostWeight = 0.9;% weightMult/nVCOs;

    % non-gated res params:
    weightMult = 1;
    respostWeight=1.5;
    resbaseWeight = 5.5;
    postc = -0.01;
    postw = baselineFreq*2*pi/1000; % Hz (times 2pi for angular freq, /1000 to convert to /s
  elseif type==11
    simprefix = sprintf('filthaftingtraj_n5000_p01_1noise_02res_%s',datestr(now,'mmmmdd'));
    loadVCOSpikesFileName = 'fig_RS2sn_filthaftingtraj_n5000_p01_1noise_02res_July13_2D_spikes_0a.txt';
    excitationClass = 2;
    useGapJunctions = 0;
    showAbstractFig = 0;
    plotType = 3;
    ncells = 5000;
    nPostIn = 5000;
    pcon = 0.01;
    useNoise = 1;
    useFilteredTrajectory = 1;
    load simple_model_RS2sn_FI_Jul11_n5000.mat
    uniqueNoiseSTD = 100*useNoise;
    g = 256;
    loadC = 'fig_RS2sn_filthaftingtraj_n5000_p01_1noise_02res_July13_C_0a.mat';

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

    % non-gated lif params:
    tau = 12; % ms
    membraneDecay = exp(-dt/tau);
    postWeight = 0.00018;
    baseWeight =  .00019;

    % gated lif params:
    basegateDur = 3; % ms
    tau = 5; % (msec)
    gatedmembraneDecay = exp(-dt/tau);
    %     weightMult = 1.4;
    gatedpostWeight = 0.00035;

    % non-gated res params:
    weightMult = .001;
    respostWeight=1e-3;
    resbaseWeight = 1e-3;
    postc = -0.01;
    postw = baselineFreq*2*pi/1000; % Hz (times 2pi for angular freq, /1000 to convert to /s
  elseif type==12
    rand('seed',1);
    randn('seed',1);
    simprefix = sprintf('filthaftingtraj_6vcoprec_n250_1noise_6vco_prec_%s',datestr(now,'mmmmdd'));
    excitationClass = 2;
    useGapJunctions = 0;
    showAbstractFig = 0;
    plotType = 2;
    nVCOs = 6;
    ncells = 250;
    nPostIn = ncells;
    cheatHDRectification = 1;
    pcon = 1;
    useNoise = 1;
    useFilteredTrajectory = 1;
    load simple_model_RS2sn_FI_Jan09_n250.mat;
    uniqueNoiseSTD = 100*useNoise;
    g = 256;
    baselineFreq = freqs(round(-1+length(freqs)/2)); % Hz
    beta = 2; % Hz/(m/s)
    
    phaseLocked = 0;
    if phaseLocked==0
      % post-stronger (precessing) params:
      % non-gated lif params:
      tau = 5; % ms
      weightMult = 1.5; 1; 0.8;
      membraneDecay = exp(-dt/tau);
      postWeight=1.4*weightMult*1/ncells/(2+nVCOs);
      baseWeight = 0.5*postWeight; 1.5*postWeight; 2*postWeight;

      % gated lif params:
      basegateDur = .05; 1; 50; % ms
      tau = 1.2*5; % (msec)
      gatedmembraneDecay = exp(-dt/tau);
      gatedpostWeight = 0.0015; 0.0014;

      % non-gated res params:
      weightMult = 7;
      respostWeight= .0035; 1.5*0.0015; 0.002;
      resbaseWeight = .002; 0.003; 0.0024;
      postc = -0.01;
      postw = baselineFreq*2*pi/1000; % Hz
    elseif phaseLocked==1
      % baseline-stronger (non-precessing) params:
      % non-gated lif params:
      tau = 5; % ms
      weightMult = 1.25; 1; 0.8;
      membraneDecay = exp(-dt/tau);
      postWeight=1*weightMult*1/ncells/(2+nVCOs);
      baseWeight = 2.5*postWeight;

      % gated lif params:
      basegateDur = .05; 1; 50; % ms
      tau = 1.2*5; % (msec)
      gatedmembraneDecay = exp(-dt/tau);
      gatedpostWeight = 1.2*0.0014;

      % non-gated res params:
      weightMult = 7;
      respostWeight= .3*0.003; 0.002;
      resbaseWeight = 1.6*0.003; 0.0024;
      postc = -0.01;
      postw = baselineFreq*2*pi/1000; % Hz
    end

    % % params for no-cheat test:
%     % non-gated lif params:
%     tau = 5; % ms
%     weightMult = .5*1.5; 1; 0.8;
%     membraneDecay = exp(-dt/tau);
%     postWeight=weightMult*1/ncells/(2+nVCOs);
%     baseWeight = 1.5*postWeight; 2*postWeight;
% 
%     % gated lif params:
%     basegateDur = .05; 1; 50; % ms
%     tau = 1.2*5; % (msec)
%     gatedmembraneDecay = exp(-dt/tau);
%     gatedpostWeight = .5*0.0014;
% 
%     % non-gated res params:
%     weightMult = 7;
%     respostWeight= .5*1.5*0.0015; 0.002;
%     resbaseWeight = 0.003; 0.0024;
%     postc = -0.01;
%     postw = baselineFreq*2*pi/1000; % Hz
  elseif type==13
    simprefix = sprintf('filthaftingtraj_6vcoprec_n1_0noise_01res%s',datestr(now,'mmmmdd'));
    excitationClass = 2;
    useGapJunctions = 0;
    showAbstractFig = 1;
    plotType = 2;
    plotVelocity = 1;
    nVCOs = 6;
    ncells = 1;
    nPostIn = ncells;
    cheatHDRectification = 1;
    pcon = 0;
    useNoise = 0;
    uniqueNoiseSTD = 100*useNoise;
    useFilteredTrajectory = 1;
    load simple_model_RS2_FI_Jan09_n1.mat;
    g = 0;
    % shared grid variables: freq = baselineFreq+beta*speed*(cos(prefHD-curHD));
    baselineFreq = freqs(round(length(freqs)/2)); % Hz
    beta = 2; % Hz/(m/s)

    % non-gated lif params:
    tau = 10; 40*7.8989/baselineFreq; 35; % ms
    membraneDecay = exp(-dt/tau);
    postWeight = .5; 0.14;
    baseWeight = .1; .2; .8;

    % gated lif params:
    basegateDur = 40; % ms
    tau = 25; % (msec)
    gatedmembraneDecay = exp(-dt/tau);
    weightMult = 3; 2;
    gatedpostWeight = weightMult/ncells/nVCOs;

    % non-gated res params:
    weightMult = 1;
    respostWeight= 3; 3; 4; 1.4; weightMult*1/ncells/(nVCOs);
    resbaseWeight = 1.5; 2; 6; 5; 2*postWeight;
    postc = -0.01;
    postw = baselineFreq*2*pi/1000; % Hz (times 2pi for angular freq, /1000 to convert to /s
  elseif type==14
    simprefix = sprintf('filthaftingtraj_6vcoprec_post25_n250_1noise_02res_%s',datestr(now,'mmmmdd'));
    excitationClass = 2;
    useGapJunctions = 1;
    showAbstractFig = 0;
    plotType = 2;
    ncells = 250;
    nVCOs = 6;
    nPostIn = 1;
    pcon = 1;
    useNoise = 1;
    useFilteredTrajectory = 1;
    load simple_model_RS2gn_FI_Jan09_n250.mat;
    
    opponentVCOInhibition = 1;
    nInhib = nPostIn;
    opponentInhibAmp = -max(currents);
    uniqueNoiseSTD = 100*useNoise;
    g = 0.1;

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

    % non-gated lif params:
    tau = 12; % ms
    membraneDecay = exp(-dt/tau);
    postWeight = 0.15/10;
    baseWeight = .07/10;

    % gated lif params:
    basegateDur = 30; % ms
    tau = 10; % (msec)
    gatedmembraneDecay = exp(-dt/tau);
    %     weightMult = 1.4;
    gatedpostWeight = 0.05;% weightMult/nVCOs;

    % non-gated res params:
    weightMult = 1;
    respostWeight=.14;
    resbaseWeight = 0.01;
    postc = -0.01;
    postw = baselineFreq*2*pi/1000; % Hz (times 2pi for angular freq, /1000 to convert to /s
  elseif type==15
    simprefix = sprintf('filthaftingtraj_n1_0noise_01res%s',datestr(now,'mmmmdd'));
    excitationClass = 2;
    useGapJunctions = 0;
    showAbstractFig = 0;
    plotType = 1;
    plotVelocity = 1;
    ncells = 1;
    nPostIn = ncells;
    pcon = 0;
    useNoise = 0;
    uniqueNoiseSTD = 100*useNoise;
    useFilteredTrajectory = 1;
    load simple_model_RS2_FI_Jan09_n1.mat;
    g = 0;
    % shared grid variables: freq = baselineFreq+beta*speed*(cos(prefHD-curHD));
    baselineFreq = freqs(round(length(freqs)/2)); % Hz
    beta = 2; % Hz/(m/s)

    % non-gated lif params:
    tau = 40*7.8989/baselineFreq; 35; % ms
    membraneDecay = exp(-dt/tau);
    postWeight = 1.2*0.14;
    baseWeight = .8;

    % gated lif params:
    basegateDur = 40; % ms
    tau = 25; % (msec)
    gatedmembraneDecay = exp(-dt/tau);
    weightMult = 1.4;
    gatedpostWeight = weightMult/ncells/nVCOs;

    % non-gated res params:
    weightMult = 1;
    respostWeight=3*1.4; weightMult*1/ncells/(nVCOs);
    resbaseWeight = .5*6; 5; 2*postWeight;
    postc = -0.01;
    postw = baselineFreq*2*pi/1000; % Hz (times 2pi for angular freq, /1000 to convert to /s
  end

  typePrefix = sprintf('%d',excitationClass);
  if ncells>1
    if useGapJunctions
      typePrefix = [typePrefix 'g'];
    else
      typePrefix = [typePrefix 's'];
    end
  end
  if useNoise
    typePrefix = [typePrefix 'n'];
  end

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

  if excitationClass<0
    d = 60;
  end

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

  stabilities = zeros(1,nruns);
  stabilities2 = zeros(1,nruns);
  for run=1:nruns
    % activity of postsynaptic I&F
    npost = 3;
    post = zeros(npost, round(simdur/dt)+1);
    if stablePost==0
      postu = zeros(npost, round(simdur/dt)+1);
      postinhib = zeros(npost, round(simdur/dt)+1);
    end
    vpost = zeros(nVCOs,round(simdur/dt)+1);
    vpostb = zeros(1,round(simdur/dt)+1);

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

    t = 0; % current time in simulation
    v = vr*ones(ncells,nVCOs); % vr*rand(ncells,1);
    u = 0*v; % initial values
    vb = vr*ones(ncells,1); % vr*rand(ncells,1);
    ub = 0*vb; % initial values

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

    % 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 4*pi/3 pi/3 pi 5*pi/3]; % VCO preferred direction(s), (radians)
    if nVCOs>length(prefHDs)
      prefHDs = repmat(prefHDs,1,ceil(nVCOs/length(prefHDs)));
    end
    prefHDs = prefHDs(1:nVCOs);

    abstractGrid = zeros(1,round(simdur/dt));
    basehist = zeros(1,round(simdur/dt));
    vcohist = zeros(nVCOs,round(simdur/dt));

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

    commonNoise = 0;

    pos = zeros(2,round(simdur/dt));
    baseI = currents(find(freqs==baselineFreq));
    I = baseI;
    ISIs = zeros(1,nVCOs);

    speed = 0; % m/s

    basegateCount = 0;

    buffer = zeros(1,nVCOs+1);

    if ~isempty(loadVCOSpikesFileName)
      loadedSpikes = textread(loadVCOSpikesFileName,'%d');
      loadedSpikes = reshape(loadedSpikes,nVCOs+1,[])';
    end

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

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

      % move virtual rat: pregenerated trajectory in variables dxs dys
      if useFilteredTrajectory
        %% filtered trajectory:
        dx = fdxs(tind);
        dy = fdys(tind);
      else
        %% unfiltered trajectory:
        dx = dxs(tind);
        dy = dys(tind);
      end
      x = x+dx;
      y = y+dy;


      speed = abs(sqrt((dx^2 + dy^2)/(dt/1000)^2)); % (m/s); abs because we include the direction via cosine later
      pos(:,tind) = [x; y];
      curHD = atan2(dy,dx);

      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;
        abstractGrid(tind) = nVCOs*cos(basePhi)+sum(cos(VCOPhi));
        basehist(tind) = basePhi;
        vcohist(:,tind) = VCOPhi;
      end

      if isempty(loadVCOSpikesFileName)
        if runNetwork
          if commonNoiseSTD
            commonNoise = commonNoiseSTD*randn;
          end
          %% active VCO networks
          if runNetworkVCOs
            for vco=1:nVCOs
              oldv = v(:,vco);
              oldu = u(:,vco);
              % much faster if we interpolate the FI curve ourself, though we do lose
              % the ability to extrapolate, which will cause an error
              % here with one of lowind or highind being empty
              desiredFreq = baselineFreq + beta*speed*(cos(prefHDs(vco)-curHD));
              lowind = find(freqs<desiredFreq,1,'last');
              highind = find(freqs>desiredFreq,1,'first');
              proportion = (desiredFreq-freqs(lowind))/(freqs(highind)-freqs(lowind));
              I = currents(lowind) + proportion*(currents(highind)-currents(lowind));
              
              inhibI = zeros(size(v(:,vco)));
              if opponentVCOInhibition
                inhibI(1:nInhib) = opponentInhibAmp.*(cos(prefHDs(vco)-curHD)<0);
              end
              
              if useGapJunctions
                if pcon==1
                  % if all-to-all connected, the gap-junctions are equivalent
                  % to sum(v(:,vco))-ncells*v(:,vco)
                  v(:,vco) = oldv + dt*(k*(oldv-vr).*(oldv-vt) - oldu + I + inhibI + commonNoise + uniqueNoiseSTD*randn(ncells,1) + g*(sum(v(:,vco))-ncells*v(:,vco)))/Cf;
                else
                  % make matrix of differences in membrane potential
                  vdiffs = C.*(repmat(oldv',ncells,1) - repmat(oldv,1,ncells));
                  v(:,vco) = oldv + dt*(k*(oldv-vr).*(oldv-vt) - oldu + I + inhibI + commonNoise + uniqueNoiseSTD*randn(ncells,1) + g*sum(vdiffs,2))/Cf;
                  %% slowest way:
                  % v(:,vco) = oldv + dt*(k*(oldv-vr).*(oldv-vt) - oldu + I + inhibI + commonNoise + uniqueNoiseSTD*randn(ncells,1) + g*sum(C.*vdiffs,2))/Cf;
                end
              else % using synaptic connections:
                if pcon==1
                  v(:,vco) = oldv + dt*(k*(oldv-vr).*(oldv-vt) - oldu + I + inhibI + commonNoise + uniqueNoiseSTD*randn(ncells,1) + g*(sum(spikes(:,vco),1)-spikes(:,vco)))/Cf;
                else
                  v(:,vco) = oldv + dt*(k*(oldv-vr).*(oldv-vt) - oldu + I + inhibI + commonNoise + uniqueNoiseSTD*randn(ncells,1) + g*(C*double(spikes(:,vco))))/Cf;
                end
              end
              u(:,vco) = oldu + dt*a*(b*(oldv-vr)-oldu);

              % save and reset spikes when v>=vpeak
              spikes(:,vco) = v(:,vco)>=vpeak;
              if excitationClass>0
                u(v(:,vco)>=vpeak,vco) = u(v(:,vco)>=vpeak,vco)+d;
              else
                u(v(:,vco)>=vpeak,vco) = d;
              end
              oldv(v(:,vco)>=vpeak) = vpeak;
              v(v(:,vco)>=vpeak,vco) = c;
              if any(spikes(:,vco))
                if mod(VCOinds(vco),1000)==0
                  VCOSpikeTimes{vco}(length(VCOSpikeTimes{vco})+1000) = 0;
                end
                VCOinds(vco) = VCOinds(vco)+1;
                VCOSpikeTimes{vco}(VCOinds(vco)) = t;
              end
            end
          end
          spikesum = sum(spikes(1:nPostIn,:),1);

          if runNetworkBaseline
            %% baseline VCO network
            oldvb = vb;
            oldub = ub;
            if useGapJunctions
              if pcon==1
                % if all-to-all connected, the gap-junctions are equivalent
                % to sum(vb)-ncells*vb
                vb = oldvb + dt*(k*(oldvb-vr).*(oldvb-vt) - oldub + baseI + commonNoise + uniqueNoiseSTD*randn(ncells,1) + g*(sum(vb)-ncells*vb))/Cf;
              else
                % make matrix of differences in membrane potential
                vbdiffs = C.*(repmat(oldvb',ncells,1) - repmat(oldvb,1,ncells));
                vb = oldvb + dt*(k*(oldvb-vr).*(oldvb-vt) - oldub + baseI + commonNoise + uniqueNoiseSTD*randn(ncells,1) + g*sum(vbdiffs,2))/Cf;
                %% slowest way:
                % vb = oldvb + dt*(k*(oldvb-vr).*(oldvb-vt) - oldub + baseI + commonNoise + uniqueNoiseSTD*randn(ncells,1) + g*sum(C.*vbdiffs,2))/Cf;
              end
            else
              if pcon==1
                vb = oldvb + dt*(k*(oldvb-vr).*(oldvb-vt) - oldub + baseI + commonNoise + uniqueNoiseSTD*randn(ncells,1) + g*(sum(spikesb,1)-spikesb))/Cf;
              else
                vb = oldvb + dt*(k*(oldvb-vr).*(oldvb-vt) - oldub + baseI + commonNoise + uniqueNoiseSTD*randn(ncells,1) + g*(C*double(spikesb)))/Cf;
              end
            end
            ub = oldub + dt*a*(b*(oldvb-vr)-oldub);

            % save and reset spikes when v>=vpeak
            spikesb = vb>=vpeak;
            %       spikesb(:,tind) = vb>=vpeak;
            spikesumb = sum(spikesb(1:nPostIn));
            if excitationClass>0
              ub(vb>=vpeak) = ub(vb>=vpeak)+d;
            else
              ub(vb>=vpeak) = d;
            end
            oldvb(vb>=vpeak) = vpeak;
            vb(vb>=vpeak) = c;
%             if any(spikesb)
%               if mod(Baseind,1000)==0
%                 BaseSpikeTimes(length(BaseSpikeTimes)+1000) = 0;
%               end
%               Baseind = Baseind+1;
%               BaseSpikeTimes(Baseind) = t;
% 
%               if baselineLIFGating
%                 basegateCount = basegateDur/dt;
%               end
%             else
%               basegateCount = basegateCount-1;
%             end
% 
%             if basegateCount>0
%               basegate = 1;
%             else
%               basegate = 0;
%             end
          end
        end
      else
        spikesum = loadedSpikes(tind,1:end-1);
        spikesumb = loadedSpikes(tind,end);
      end

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

        if baselineLIFGating
          basegateCount = basegateDur/dt;
        end
      else
        basegateCount = basegateCount-1;
      end

      if basegateCount>0
        basegate = 1;
      else
        basegate = 0;
      end

      
      if ~isempty(saveVCOSpikesFileName)
        % each line: number of spikes from 1:nPostIn on this
        % time step for each VCO and the baseline
        % we'll buffer the output to speed things up
        buflin = [];
        if runNetworkVCOs
          for vcoi=1:nVCOs
            buflin = [buflin sum(spikes(1:nPostIn,vcoi))];
          end
        end
        if runNetworkBaseline
          buflin = [buflin sum(spikesb(1:nPostIn))];
        end
        buffer = [buffer; buflin];
        if size(buffer,1)==50
          fid = fopen(saveVCOSpikesFileName,'a');
          for br=1:size(buffer,1)
            for bc=1:size(buffer,2)
              fprintf(fid,'%d ',buffer(br,bc));
            end
            fprintf(fid,'\n');
          end
          fclose(fid);
          buffer = [];
        end
      end

      state(:,tind) = [v(1,:)'; vb(1)];

      if cheatHDRectification
        spikesum = spikesum.*(cos(prefHDs-curHD)>0);
      end
      
      % post cell 1 = non-gated lif
      if stablePost
        post(1,tind) = post(1,tind-1)*membraneDecay + postWeight*sum(spikesum)+baselineMult*baseWeight*spikesumb;
        if post(1,tind)>1
          post(1,tind) = -1e-10;
        end
      else
        GABAtau = 15; % ms
        % post cell 1 = non-gated simple model
        postinhib(1,tind) = postinhib(1,tind-1)*exp(-dt/GABAtau) + postWeight*sum(spikesum)+baselineMult*baseWeight*spikesumb;
        oldv = post(1,tind-1);
        oldu = postu(1,tind-1);
        post(1,tind) = oldv + dt*(k*(oldv-vr).*(oldv-vt) - oldu + baseI + postinhib(1,tind))/Cf;
        postu(1,tind) = oldu + dt*a*(b*(oldv-vr)-oldu);
        if post(1,tind)>vpeak
          postu(1,tind) = postu(1,tind) + d;
          post(1,tind-1) = vpeak;
          post(1,tind) = c;
        end
      end
      % post cell 2 = gated lif
      post(2,tind) = post(2,tind-1)*gatedmembraneDecay + basegate*gatedpostWeight*sum(spikesum);
      if post(2,tind)>1
        post(2,tind) = -1e-10;
      end
      % post cell 3 = non-gated resonant
      post(3,tind) = post(3,tind-1) + dt*((postc+postw*i)*post(3,tind-1) + respostWeight*sum(spikesum)+baselineMult*resbaseWeight*spikesumb);
      postspikes(tind) = real(post(3,tind))>1;
      if real(post(3,tind))>1
        post(3,tind) = i;
      end

      vpostb(tind) =  vpostb(tind-1)*membraneDecay + baseWeight*spikesumb;
      for vind=1:nVCOs
        vpost(vind,tind) = vpost(vind,tind-1)*membraneDecay + basegate*postWeight*spikesum(vind);
      end

      if runningPlot
        figure(runh);
        plot(pos(1,:),pos(2,:),'.',pos(1,find(postspikes>0)),pos(2,find(postspikes>0)),'R.'), title('network VCO spatial firing'); xlabel('Position (m)','fontsize',12), ylabel('Position (m)','fontsize',12);
        set(gca,'xlim',[-1 1],'ylim',[-1 1]);
        drawnow
      end
    end
    toc
    %% simulation is complete, now for some post-processing and plotting:

    if ~isempty(loadVCOSpikesFileName)
      for vcoi=1:nVCOs
        VCOSpikeTimes{vcoi} = dt*find(loadedSpikes(:,vcoi)>0);
      end
      BaseSpikeTimes = dt*find(loadedSpikes(:,end)>0);
    end
    
    % flush spike times buffer
    if ~isempty(saveVCOSpikesFileName)
      if size(buffer,1)>0
        fid = fopen(saveVCOSpikesFileName,'a');
        for br=1:size(buffer,1)
          for bc=1:size(buffer,2)
            fprintf(fid,'%d ',buffer(br,bc));
          end
        end
        fprintf(fid,'\n');
        fclose(fid);
        buffer = [];
      end
    end

    clear C

    if runNetwork
      for ce=1:npost
        spikeTimes = find(post(ce,:)==-1e-10);

        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

  % transform the simple model output so that our LIF spike-detection in
  % the figures works
  if stablePost==0
    post(post==vpeak) = -1e-10;
  end

  % tweak straight line trajectories to avoid error below
  if min(pos(1,:))==max(pos(1,:))
    pos(1,end) = pos(1,end)+eps;
  end
  if min(pos(2,:))==max(pos(2,:))
    pos(2,end) = pos(2,end)+eps;
  end
  

  %% FIGURES from here to the end of the file
  %% (and they are a bit of a mess!)
  % let's set up a few aspects of the appearance:
  set(0,'defaulttextfontsize',12);
  if journalChargesALotForColor
    % for images, white is largest value
    figure; set(0,'defaultfigurecolormap',flipud(gray)); close
    % we only ever plot <=3 colors at one time in the main figures
    set(0,'defaultaxescolororder',[0 0 1; 0 0.5 0; 1 0 0;]);
    %     set(0,'defaultaxescolororder',[0 0 0; .75 .75 .75]);
    %     set(0,'defaultaxeslinestyleorder',{'-','--'});
    %     % this means line 1 is solid/black, line 2 is solid/gray, line 3 dashed/black, line 4 dashed/gray
  end

  if plotVCOsAndPosts
    % Note: this won't plot if using loaded spikes, have to save this when
    % the VCO spikes are first generated!
    figure; plot(state(:,1:40000)'); set(gcf,'position',[82 404 3*257 125])
    xlabel('Time (ms)','fontsize',12); ylabel('Voltage (mV)','fontsize',12)
    set(gca,'ylim',[-60 50]);
    set(gca,'fontsize',12)
    set(gca,'box','off');
    if saveFigures && exist(['fig_RS' typePrefix '_' simprefix '_2D_vcos_trace_0a.eps'])~=2
      print_eps(['fig_RS' typePrefix '_' simprefix '_2D_vcos_trace_0a.eps'])
      saveas(gcf,['fig_RS' typePrefix '_' simprefix '_2D_vcos_trace_0a.fig'])
    end

    plotduration = 4000; % ms
    % plot non-gated lif post:
    figure; plot(dt:dt:plotduration,post(1,1:(plotduration/dt)));
    if any(post(1,1:(plotduration/dt))==-1e-10)
      hold on; plot(dt*find(post(1,1:(plotduration/dt))==-1e-10),1,'r.','MarkerSize',16);
    end
    set(gcf,'position',[82 404 3*257 125])
    xlabel('Time (ms)','fontsize',12); ylabel('Activity','fontsize',12)
    if stablePost
      set(gca,'ylim',[-0.1 1.1]);
    end
    set(gca,'fontsize',12)
    set(gca,'box','off');
    if saveFigures && exist(['fig_RS' typePrefix '_' simprefix '_2D_post_nogatelif_trace_0a.eps'])~=2
      print_eps(['fig_RS' typePrefix '_' simprefix '_2D_post_nogatelif_trace_0a.eps'])
      saveas(gcf,['fig_RS' typePrefix '_' simprefix '_2D_post_nogatelif_trace_0a.fig'])
    end

    % plot gated lif post:
    if nVCOs==1
      figure; plot(dt:dt:plotduration,post(2,1:(plotduration/dt))==-1e-10);
    else
%       figure; plot(dt:dt:plotduration,post(2,1:(plotduration/dt))==-1e-10);
      figure; plot(dt:dt:plotduration,post(2,1:(plotduration/dt)));      
    end
    if any(post(2,1:(plotduration/dt))==-1e-10)
      hold on; plot(dt*find(post(2,1:(plotduration/dt))==-1e-10),1,'r.','MarkerSize',16);
    end
    set(gcf,'position',[82 404 3*257 125])
    xlabel('Time (ms)','fontsize',12); ylabel('Activity','fontsize',12)
    set(gca,'ylim',[-0.1 1.1]);
    set(gca,'box','off');
    set(gca,'fontsize',12)
    if saveFigures && exist(['fig_RS' typePrefix '_' simprefix '_2D_post_gating_trace_0a.eps'])~=2
      print_eps(['fig_RS' typePrefix '_' simprefix '_2D_post_gating_trace_0a.eps'])
      saveas(gcf,['fig_RS' typePrefix '_' simprefix '_2D_post_gating_trace_0a.fig'])
    end

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

  if plotType==0
    % these in normalized coords:
    trajH = 3/4*1/2;
    trajW = 3/4*1/3;
    autoH = trajH;
    autoW = trajW;
    diffH = 2/3*1/3;
    diffW = trajW;

    trajL = 1/4*1/3;
    trajB = 1/6*1/2;
    autoL = 1/6*1/3;
    autoB = trajB;
    diffL = 1/5*1/3;
    diffB1 = .007+1/8*1/2;
    diffB2 = diffB1+0.005;
    diffB3 = trajB;

    figure('position',[82 404 647 484])
    if runNetwork
      posskip = 500;
      thr = 1.5;
      if stablePost
        subplot('position',[trajL+0 trajB+1/2 trajW trajH]); plot(pos(1,1:posskip:length(pos)),pos(2,1:posskip:length(pos)),'.','MarkerSize',8); hold on; plot(pos(1,find(post(2,:)==-1e-10)),pos(2,find(post(2,:)==-1e-10)),'R.','MarkerSize',15)
      else
        subplot('position',[trajL+0 trajB+1/2 trajW trajH]); plot(pos(1,1:posskip:length(pos)),pos(2,1:posskip:length(pos)),'.','MarkerSize',8); hold on; plot(pos(1,find(post(1,:)==-1e-10)),pos(2,find(post(1,:)==-1e-10)),'R.','MarkerSize',15)
      end
      if ~showAbstractFig
        title(['Spatial firing - 0 to ' num2str(simdur/1000) ' s '],'fontsize',12)
      else
        title('Network model spatial firing','fontsize',12)
      end
      xlabel('Position (m)','fontsize',12), ylabel('Position (m)','fontsize',12)
      axis equal
      set(gca,'xlim',[min(pos(1,:)) max(pos(1,:))],'ylim',[min(pos(2,:)) max(pos(2,:))])
      set(gca,'box','off');
      set(gca,'fontsize',12)

      edges{1} = linspace(min(pos(1,:)),max(pos(1,:)),50);
      edges{2} = linspace(min(pos(2,:)),max(pos(2,:)),50);
      if stablePost
        rate = hist3([pos(1,find(post(2,:)==-1e-10)); pos(2,find(post(2,:)==-1e-10))]','Edges',edges);
      else
        rate = hist3([pos(1,find(post(1,:)==-1e-10)); pos(2,find(post(1,:)==-1e-10))]','Edges',edges);
      end
      subplot('position',[autoL+1/3 autoB+1/2 autoW autoH]); imagesc(rot90(conv2(rate,rate,'same')))
      if showAbstractFig
        title('Network model autocorrelogram','fontsize',12)
      elseif stablePost
        title(['Spatial autocorrelogram - 0 to ' num2str(simdur/1000) ' s '],'fontsize',12)
      end
      set(gca,'xtick',[],'ytick',[]);
      set(gca,'fontsize',12)
      axis tight

      if ~showAbstractFig
        subplot('position',[trajL+0 trajB+0 trajW trajH]); plot(pos(1,1:posskip:length(pos)/4),pos(2,1:posskip:length(pos)/4),'.','MarkerSize',8);
        hold on;
        plot(pos(1,find(post(2,1:length(pos)/4)==-1e-10)),pos(2,find(post(2,1:length(pos)/4)==-1e-10)),'R.','MarkerSize',15)
        title('Spatial firing - 0 to 40 s ','fontsize',12)
        title(['Spatial firing - 0 to ' num2str(simdur/1000/4) ' s '],'fontsize',12)
        xlabel('Position (m)','fontsize',12), ylabel('Position (m)','fontsize',12)
        set(gca,'fontsize',12)
        axis equal
        set(gca,'xlim',[min(pos(1,:)) max(pos(1,:))],'ylim',[min(pos(2,:)) max(pos(2,:))])
        set(gca,'box','off');
        rate = hist3([pos(1,find(post(2,1:length(pos)/4)==-1e-10)); pos(2,find(post(2,1:length(pos)/4)==-1e-10))]','Edges',edges);
        subplot('position',[autoL+1/3 autoB+0 autoW autoH]); imagesc(rot90(conv2(rate,rate,'same')))
        title(['Spatial autocorrelogram - 0 to ' num2str(simdur/1000/4) ' s '],'fontsize',12)
        set(gca,'xtick',[],'ytick',[]);
        set(gca,'fontsize',12)
        axis tight
      end


      if ncells>1 && useGapJunctions==0 && useNoise>0
        % to keep filesizes lower when many noisy cells are synaptically
        % coupled:
        spikeskips = ncells/11;
      else
        spikeskips = 1;
      end
      if runNetworkVCOs
        if nVCOs>1
          subplot('position',[diffL+2/3 diffB1+2/3 diffW diffH]); plot(mod(vcohist(1,round(VCOSpikeTimes{1}(1:spikeskips:end)/dt)),2*pi),'.')
        else
          subplot('position',[diffL+2/3 trajB+1/2 diffW trajH]); plot(mod(vcohist(1,round(VCOSpikeTimes{1}(1:spikeskips:end)/dt)),2*pi),'.')
        end
        if plotVelocity
          hold on; plot(5-1e4*fdxs(round(VCOSpikeTimes{1}(1:spikeskips:end)/dt)),'k'); set(gca,'ylim',[0 6.2]);
        end
        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,'box','off');
        set(gca,'fontsize',12)
        if nVCOs>1
          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,'box','off');
          set(gca,'fontsize',12)
        end
      end

      if runNetworkBaseline
        if nVCOs>1
          subplot('position',[diffL+2/3 diffB3+0 diffW diffH]); plot(mod(basehist(round(BaseSpikeTimes(1:spikeskips:end)/dt)),2*pi),'.')
        else
          subplot('position',[diffL+2/3 diffB3+0 diffW trajH]); plot(mod(basehist(round(BaseSpikeTimes(1:spikeskips:end)/dt)),2*pi),'.')
        end
        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,'box','off');
        set(gca,'fontsize',12)
      end
    end
    if runAbstract && showAbstractFig
      spikeskip = 100;
      subplot('position',[trajL+0 trajB+0 trajW trajH]); plot(pos(1,1:posskip:length(pos)),pos(2,1:posskip:length(pos)),'.','MarkerSize',8);
      hold on; plot(pos(1,spikeskip*find(abstractGrid(1:spikeskip:end)>abstractThr)),pos(2,spikeskip*find(abstractGrid(1:spikeskip:end)>abstractThr)),'R.','MarkerSize',15)
      title('Abstract model spatial firing','fontsize',12)
      xlabel('Position (m)','fontsize',12), ylabel('Position (m)','fontsize',12)
      axis equal
      set(gca,'xlim',[min(pos(1,:)) max(pos(1,:))],'ylim',[min(pos(2,:)) max(pos(2,:))])
      set(gca,'box','off');
      set(gca,'fontsize',12)

      rate = hist3([pos(1,find(abstractGrid>abstractThr)); pos(2,find(abstractGrid>abstractThr))]','Edges',edges);
      %   rate = hist3([pos(1,find(abstractGrid>abstractThr)); pos(2,find(abstractGrid>abstractThr))]',[50 50]);
      subplot('position',[autoL+1/3 autoB+0 autoW autoH]); imagesc(rot90(conv2(rate,rate,'same')))
      title('Abstract model autocorrelogram','fontsize',12)
      set(gca,'xtick',[],'ytick',[]);
      set(gca,'fontsize',12)
      axis tight
    end
    if saveFigures && exist(['fig_RS' typePrefix '_' simprefix '_2D_all_0a.eps'])~=2
      save(['fig_RS' typePrefix '_' simprefix '_2D_variables_0a.mat'])
      % free some memory:
      %   clear pos abstractGrid vhist vcohist basehist post
      print_eps(['fig_RS' typePrefix '_' simprefix '_2D_all_0a.eps'])
      saveas(gcf,['fig_RS' typePrefix '_' simprefix '_2D_all_0a.fig'])
      close
    end
  elseif plotType==1
    % these in normalized coords:
    trajH = 3/4*1/2;
    trajW = 3/4*1/3;
    diffH = 2/3*1/3;
    diffW = trajW;

    col1L = 1/4*1/3;
    col2L = 1/3+1/5*1/3;
    col3L = 2/3+1/5*1/3;
    row2B = 1/6*1/2;
    diffL = 1/5*1/3;
    diffB1 = .007+1/8*1/2;
    diffB2 = diffB1+0.005;

    if ncells>1 && useNoise>0 % && useGapJunctions==0
      % to keep filesizes lower when many noisy cells are synaptically
      % coupled:
      spikeskips = ncells/23;
    else
      spikeskips = 1;
    end

    mainplot = 2;
    if mainplot==1
      if stablePost
        plottype = 'Integrator postsynaptic cell';
      else
        plottype = 'Inhibitory oscillators';
      end
    else
      plottype = 'Gated postsynaptic cell';
    end
    figure('position',[82 404 647 484])
    posskip = 500;
    if nVCOs==2
      % gated lif:
      subplot('position',[col1L+0 row2B+1/2 trajW trajH]);
      plot(pos(1,1:posskip:length(pos)),pos(2,1:posskip:length(pos)),'.','MarkerSize',8);
      hold on; plot(pos(1,find(post(mainplot,:)==-1e-10)),pos(2,find(post(mainplot,:)==-1e-10)),'R.','MarkerSize',15)
      title(plottype,'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,'box','off');
      set(gca,'fontsize',12)

      % 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,find(post(mainplot,:)==-1e-10)); pos(2,find(post(mainplot,:)==-1e-10))]','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([plottype ' autocorrelogram'],'fontsize',12)
      end
      set(gca,'xtick',[],'ytick',[]);
      set(gca,'fontsize',12)
      axis tight

      if showAbstractFig
        % abstract:
        abstractThr = 3;
        subplot('position',[col1L+0 row2B+0 trajW trajH]);
        plot(pos(1,1:posskip:length(pos)),pos(2,1:posskip: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,'box','off');
        set(gca,'fontsize',12)

        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
      else
        % non-gated res:
        subplot('position',[col1L+0 row2B+0 trajW trajH]);
        plot(pos(1,1:posskip:length(pos)),pos(2,1:posskip:length(pos)),'.','MarkerSize',8);
        hold on; plot(pos(1,find(postspikes>0)),pos(2,find(postspikes>0)),'R.','MarkerSize',15)
        title('Resonator postsynaptic 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,'box','off');
        set(gca,'fontsize',12)

        rate = hist3([pos(1,find(postspikes==1)); pos(2,find(postspikes==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:
      subplot('position',[col3L diffB2+2/3 diffW diffH]);
      plot(mod(vcohist(1,round(VCOSpikeTimes{1}(1:spikeskips:end)/dt)),2*pi),'.')
      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,'box','off');
      set(gca,'fontsize',12)

      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,'box','off');
      set(gca,'fontsize',12)

      % 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,'box','off');
      set(gca,'fontsize',12)
    end
    if saveFigures && exist(['fig_RS' typePrefix '_' simprefix '_2D_all_0a.eps'])~=2
      save(['fig_RS' typePrefix '_' simprefix '_2D_variables_0a.mat'])
      % free some memory:
      %   clear pos abstractGrid vhist vcohist basehist post
      print_eps(['fig_RS' typePrefix '_' simprefix '_2D_all_0a.eps'])
      saveas(gcf,['fig_RS' typePrefix '_' simprefix '_2D_all_0a.fig'])
      close
    end
  elseif plotType==2
    % these in normalized coords:
    trajH = 3/4*1/2; % <1/2
    trajW = 3/4*1/3; % <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

    figure('position',[82 404 647 484])
    posskip = 500;
    spikeskip = 1;
    % non-gated lif:
    subplot('position',[col1L+0 row2B+1/2 trajW trajH]);
    plot(pos(1,1:posskip:length(pos)),pos(2,1:posskip:length(pos)),'.','MarkerSize',8);
    hold on; plot(pos(1,spikeskip*find(post(1,1:spikeskip:end)==-1e-10)),pos(2,spikeskip*find(post(1,1:spikeskip:end)==-1e-10)),'R.','MarkerSize',15)
    if stablePost
      title('Integrator postsynaptic cell','fontsize',12)
    else
      title('Inhibitory oscillators','fontsize',12)
    end
    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,'box','off');
    set(gca,'fontsize',12)

    % gated lif:
    subplot('position',[col2L+0 row2B+1/2 trajW trajH]);
    plot(pos(1,1:posskip:length(pos)),pos(2,1:posskip:length(pos)),'.','MarkerSize',8);
    hold on; plot(pos(1,spikeskip*find(post(2,1:spikeskip:end)==-1e-10)),pos(2,spikeskip*find(post(2,1:spikeskip:end)==-1e-10)),'R.','MarkerSize',15)
    title('Gated postsynaptic 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,'box','off');
    set(gca,'fontsize',12)

    % non-gated res:
    subplot('position',[col3L+0 row2B+1/2 trajW trajH]);
    plot(pos(1,1:posskip:length(pos)),pos(2,1:posskip:length(pos)),'.','MarkerSize',8);
    hold on; plot(pos(1,spikeskip*find(postspikes(1:spikeskip:end)>0)),pos(2,spikeskip*find(postspikes(1:spikeskip:end)>0)),'R.','MarkerSize',15)
    title('Resonator postsynaptic 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,'box','off');
    set(gca,'fontsize',12)

    % abstract:
    if nVCOs==1
      abstractThr = 1.8;
    else
      abstractThr = 3;
    end
    spikeskip = 100;
    subplot('position',[col1L+0 row2B+0 trajW trajH]);
    plot(pos(1,1:posskip:length(pos)),pos(2,1:posskip:length(pos)),'.','MarkerSize',8);
    hold on; plot(pos(1,spikeskip*find(abstractGrid(1:spikeskip:end)>abstractThr)),pos(2,spikeskip*find(abstractGrid(1:spikeskip:end)>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,'box','off');
    set(gca,'fontsize',12)

    % vco 1 phase diff:
    spikeskips = 1;
    subplot('position',[col2L+0 row2B+0 trajW trajH]);
    plot(1:spikeskips:length(VCOSpikeTimes{1}),mod(vcohist(1,round(VCOSpikeTimes{1}(1:spikeskips:end)/dt)),2*pi),'.')
    title('VCO 1 phase drift','fontsize',12)
    xlabel('Time','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,'box','off');
    set(gca,'fontsize',12)

    % base phase diff:
    subplot('position',[col3L row2B+0 trajW trajH]);
    plot(1:spikeskips:length(BaseSpikeTimes),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,'box','off');
    set(gca,'fontsize',12)

    if saveFigures && exist(['fig_RS' typePrefix '_' simprefix '_2D_all_0a.eps'])~=2
      save(['fig_RS' typePrefix '_' simprefix '_2D_variables_0a.mat'])
      % if print_eps complains, you can free some memory: (or even clear all)
      %   clear pos abstractGrid vhist vcohist basehist post
      print_eps(['fig_RS' typePrefix '_' simprefix '_2D_all_0a.eps'])
      saveas(gcf,['fig_RS' typePrefix '_' simprefix '_2D_all_0a.fig'])
      close
    end
  elseif plotType==3
    % these in normalized coords:
    trajH = 3/4*1/3; % <1/2
    trajW = 3/4*1/3; % <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/3; % <1/3-trajH

    figure('position',[82 404 647 484])
    posskips = 500;
    %     if nVCOs==1
    % non-gated lif:
    subplot('position',[col1L+0 row2B+2/3 trajW trajH]);
    plot(pos(1,1:posskips:length(pos)),pos(2,1:posskips:length(pos)),'.','MarkerSize',8);
    hold on; plot(pos(1,find(post(1,:)==-1e-10)),pos(2,find(post(1,:)==-1e-10)),'R.','MarkerSize',15)
    if stablePost
      title('Integrator postsynaptic cell','fontsize',12)
    else
      title('Inhibitory oscillators','fontsize',12)
    end
    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,'box','off');
    set(gca,'fontsize',12)

    edges{1} = linspace(min(pos(1,:)),max(pos(1,:)),50);
    edges{2} = linspace(min(pos(2,:)),max(pos(2,:)),50);
    rate = hist3([pos(1,find(post(1,:)==-1e-10)); pos(2,find(post(1,:)==-1e-10))]','Edges',edges);
    subplot('position',[col1L+0 row2B+1/3 trajW trajH]); imagesc(rot90(conv2(rate,rate,'same')))
    title('Integrator autocorrelogram','fontsize',12)
    set(gca,'xtick',[],'ytick',[]);
    axis tight
    set(gca,'fontsize',12)

    % gated lif:
    subplot('position',[col2L+0 row2B+2/3 trajW trajH]);
    plot(pos(1,1:posskips:length(pos)),pos(2,1:posskips:length(pos)),'.','MarkerSize',8);
    hold on; plot(pos(1,find(post(2,:)==-1e-10)),pos(2,find(post(2,:)==-1e-10)),'R.','MarkerSize',15)
    title('Gated, integrator postsynaptic 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,'box','off');
    set(gca,'fontsize',12)

    rate = hist3([pos(1,find(post(2,:)==-1e-10)); pos(2,find(post(2,:)==-1e-10))]','Edges',edges);
    subplot('position',[col2L+0 row2B+1/3 trajW trajH]); imagesc(rot90(conv2(rate,rate,'same')))
    title('Gated, integrator autocorrelogram','fontsize',12)
    set(gca,'xtick',[],'ytick',[]);
    axis tight
    set(gca,'fontsize',12)

    % non-gated res:
    subplot('position',[col3L+0 row2B+2/3 trajW trajH]);
    plot(pos(1,1:posskips:length(pos)),pos(2,1:posskips:length(pos)),'.','MarkerSize',8);
    hold on; plot(pos(1,find(postspikes>0)),pos(2,find(postspikes>0)),'R.','MarkerSize',15)
    title('Resonator postsynaptic 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,'box','off');
    set(gca,'fontsize',12)

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

    % abstract:
    if nVCOs==1
      abstractThr = 1.8;
    elseif nVCOs==2
      abstractThr = 3;
    else
      abstractThr = nVCOs*1.5;
    end
    subplot('position',[col1L+0 row2B+0 trajW trajH]);
    plot(pos(1,1:posskips:length(pos)),pos(2,1:posskips: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,'box','off');
    set(gca,'fontsize',12)

    % vco 1 phase diff:
    if ncells<5000
      spikeskips = 1;
    else
      spikeskips = 25;
    end
    subplot('position',[col2L+0 row2B+0 trajW trajH]);
    plot(mod(vcohist(1,round(VCOSpikeTimes{1}(1:spikeskips:end)/dt)),2*pi),'.')
    title('VCO 1 phase drift','fontsize',12)
    xlabel('Time','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,'box','off');
    set(gca,'fontsize',12)

    % base phase diff:
    subplot('position',[col3L row2B+0 trajW trajH]);
    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,'box','off');
    set(gca,'fontsize',12)

    if saveFigures && exist(['fig_RS' typePrefix '_' simprefix '_2D_all_0a.eps'])~=2
      save(['fig_RS' typePrefix '_' simprefix '_2D_variables_0a.mat'])
      % if print_eps complains, you can free some memory: (or even clear all)
      %   clear pos abstractGrid vhist vcohist basehist post
      print_eps(['fig_RS' typePrefix '_' simprefix '_2D_all_0a.eps'])
      saveas(gcf,['fig_RS' typePrefix '_' simprefix '_2D_all_0a.fig'])
      close
    end
  end
end

% return

% plot spike phase with respect to baseline vs time
% spikeindsa = find(post(1,:)==-1e-10); name = 'LIF grid cell';
% spikeindsa = find(post(2,:)==-1e-10); name = 'Gated LIF grid cell';
[pks, spikeindsa] = find(postspikes>0); name = 'RAF grid cell';
spiketimesa = dt*spikeindsa; % ms 
spiketimesb = BaseSpikeTimes(diff(BaseSpikeTimes)>20); % ms
xs = cumsum(fdxs); % m
ys = cumsum(fdys); % m
spikeCoords = [xs(spikeindsa)' ys(spikeindsa)']; % m
spikePhases = zeros(1,length(spikeCoords));
% for each spike of vco 1
for spi=1:length(spiketimesa)
lowerbound = find(spiketimesb<spiketimesa(spi),1,'last');
upperbound = lowerbound+1;
if isempty(lowerbound) || lowerbound==length(spiketimesb)
spikePhases(spi) = NaN;
else
spikePhases(spi) = 2*pi*(spiketimesa(spi)-spiketimesb(lowerbound))/(spiketimesb(upperbound)-spiketimesb(lowerbound)); % rad
end
end
spikePhases = spikePhases(:);
figure; plot(spiketimesa,spikePhases,'.'); title(name); set(gca,'ylim',[0 2*pi]);
figure; plot(spikeCoords(:,1),spikePhases,'.'); title(name); set(gca,'ylim',[0 2*pi]);

spikeindsa = find(post(1,:)==-1e-10); name = 'LIF grid cell';
spiketimesa = dt*spikeindsa; % ms 
spiketimesb = BaseSpikeTimes(diff(BaseSpikeTimes)>20); % ms
xs = cumsum(fdxs); % m
ys = cumsum(fdys); % m
spikeCoords = [xs(spikeindsa)' ys(spikeindsa)']; % m
spikePhases = zeros(1,length(spikeCoords));
% for each spike of vco 1
for spi=1:length(spiketimesa)
lowerbound = find(spiketimesb<spiketimesa(spi),1,'last');
upperbound = lowerbound+1;
if isempty(lowerbound) || lowerbound==length(spiketimesb)
spikePhases(spi) = NaN;
else
spikePhases(spi) = 2*pi*(spiketimesa(spi)-spiketimesb(lowerbound))/(spiketimesb(upperbound)-spiketimesb(lowerbound)); % rad
end
end
spikePhases = spikePhases(:);
figure; plot(spiketimesa,spikePhases,'.'); title(name); set(gca,'ylim',[0 2*pi]);
figure; plot(spikeCoords(:,1),spikePhases,'.'); title(name); set(gca,'ylim',[0 2*pi]);

Loading data, please wait...