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