Added empty data folder

.gitignore

@@ -0,0 +1,6 @@
+# Folders
+data/*
+
+# Files
+matlab/pp_dir.txt
+!.gitkeep

matlab/@pricklypear/anadata.m

@@ -0,0 +1,59 @@
+function output = anadata(obj, varargin)
+%ANADATA Get the analog data from a PricklyPear run
+%   anadata(PP) returns the analog data (membrane potential) of the soma
+%   from a PricklyPear simulation run (therefore, PP needs an experiment
+%   and irec).
+%
+%   anadata(PP, channel) specifies the recording channel (c, d0d, ad...).
+%
+%   anadata(PP, channel, metric) specifies the metric (vm).
+%
+%   anadata(..., 'spt', spt) removes the spikes at spt (ms) by using
+%   linear interpolation over a [-1 1.2] ms window.
+%   
+%   anadata(..., 'spwin', [-1 1.2]) specifies the time window (ms) to be used
+%   for linear interpolation when suppressing spikes.
+
+    % Handle input
+    p = inputParser();
+    p.addOptional('channel', 'c', @ischar);
+    p.addOptional('metric', 'vm', @ischar);
+    p.addParamValue('spt', []);
+    p.addParamValue('spwin', [-1 1.2]);
+    p.parse(varargin{:});
+    p = p.Results;
+    
+    % The object must be a specific run
+    assert(~isempty(experiment(obj)), 'Please specify an experiment');
+    assert(~isempty(irec(obj)), 'Please specify a specific run');
+    
+    % Requested variables must be correct
+    assert(any(strcmp([obj.sim_channels {'all'}], p.channel)), sprintf('"%s" is not a valid channel', p.channel));
+    assert(any(strcmp([obj.sim_metrics {'all'}], p.metric)), sprintf('"%s" is not a valid metric', p.metric));
+    
+    % Handle output
+    output = {};
+    
+    % Generate the filename
+    if obj.irec == 0
+        fname = fullfile(obj.datadir(), sprintf('run_unsaved__vm_%s.txt', p.channel));
+    else
+        fname = fullfile(obj.datadir(), sprintf('run_%05i__vm_%s.txt', obj.irec, p.channel));
+    end
+    
+    % Read it if it exists
+    if exist(fname, 'file')
+        txt = fileread(fname);
+        data = textscan(txt, '%f', 'Delimiter', '\n');
+        output = data{1};
+        output(end) = [];
+    else
+        error('Data was not found');
+    end
+    
+    % Remove spikes
+    if ~isempty(p.spt)
+        output = cutFromWaveform(obj.dt, output, p.spt, p.spwin);
+    end
+    
+end

matlab/@pricklypear/basedir.m

@@ -0,0 +1,18 @@
+function output = basedir(obj)
+%DATADIR Get the pricklypear base directory
+%   directdir(PP)
+%
+%   See also pricklypear/datadir, pricklypear/modeldir
+    
+    fn = mfilename('fullpath');
+    fn = fileparts(fn);
+    fn = fileparts(fn);
+    fn = fullfile(fn, 'pp_dir.txt');
+    
+    fid = fopen(fn,'r');
+    txt = fscanf(fid, '%s');
+    fclose(fid);
+    
+    output = txt;
+
+end

matlab/@pricklypear/child.m

@@ -0,0 +1,11 @@
+function output = child(obj, irec)
+%CHILD Spawn a child of a pricklypear object
+%   PP.irec -> []
+%   P = child(PP, 2)
+%   P.irec -> 2
+%
+%   See also pricklypear/parent
+
+    output = pricklypear(experiment(obj), irec);
+
+end

matlab/@pricklypear/datadir.m

@@ -0,0 +1,15 @@
+function output = datadir(obj)
+%DATADIR Get the directory where data is stored
+%   datadir(PP) returns the directory where data is stored if the object
+%   has a specified experiment, otherwise it returns the direct directory
+%   (where the model directly outputs its data).
+%
+%   See also pricklypear/directdir, pricklypear/modeldir
+
+    if isempty(obj.sim_experiment)
+        output = obj.directdir();
+    else
+        output = fullfile(obj.basedir, 'data', obj.sim_experiment);
+    end
+
+end

matlab/@pricklypear/delete.m

@@ -0,0 +1,91 @@
+function output = delete(obj, varargin)
+%DELETE Delete the data of a run
+%   delete(PP) deletes the data of run PP if PP has a specified experiment
+%   and irec.
+
+    % Handle input
+    p = inputParser();
+    p.addParamValue('directMode', false);
+    p.parse(varargin{:});
+    p = p.Results;
+    
+    % Direct mode
+    fpath = obj.datadir();
+    irec = obj.irec;
+    if p.directMode
+        fpath = obj.directdir();
+        irec = obj.iter;
+    end
+
+    % The object must be a specific run
+    assert(~isempty(experiment(obj)), 'Please specify an experiment');
+    assert(~isempty(irec), 'Please specify a specific run');
+    
+    % Handle output
+    isDeleted = 0;
+    
+    % Variables
+    flExt = obj.sim_channels;
+    
+    % Delete every single channel known
+    for ii = 1:numel(flExt)
+        % Generate a filename
+        if irec == 0
+            fname = fullfile(fpath, sprintf('run_unsaved__vm_%s.txt', flExt{ii}));
+        else
+            fname = fullfile(fpath, sprintf('run_%05i__vm_%s.txt', irec, flExt{ii}));
+        end
+        
+        % Check if file exists
+        if exist(fname, 'file')
+            delete(fname);
+            isDeleted = 1;
+        end
+    end
+    % Delete other files
+    if irec == 0
+        fname = fullfile(fpath, sprintf('run_unsaved__evt_contra.txt'));
+        if exist(fname, 'file')
+            delete(fname);
+            isDeleted = 1;
+        end
+        fname = fullfile(fpath, sprintf('run_unsaved__evt_ipsi.txt'));
+        if exist(fname, 'file')
+            delete(fname);
+            isDeleted = 1;
+        end
+        fname = fullfile(fpath, sprintf('run_unsaved__properties.ppbin'));
+        if exist(fname, 'file')
+            delete(fname);
+            isDeleted = 1;
+        end
+    else
+        fname = fullfile(fpath, sprintf('run_%i__evt_contra.txt', irec));
+        if exist(fname, 'file')
+            delete(fname);
+            isDeleted = 1;
+        end
+        fname = fullfile(fpath, sprintf('run_%i__evt_ipsi.txt', irec));
+        if exist(fname, 'file')
+            delete(fname);
+            isDeleted = 1;
+        end
+        fname = fullfile(fpath, sprintf('run_%i__properties.ppbin', irec));
+        if exist(fname, 'file')
+            delete(fname);
+            isDeleted = 1;
+        end
+    end
+    
+    % Output
+    output = isDeleted;
+    if nargout == 0
+        switch isDeleted
+            case 0
+                fprintf('Run %i is not deleted\n', irec);
+            case 1
+                fprintf('Run %i is deleted\n', irec);
+        end
+    end
+
+end

matlab/@pricklypear/directdir.m

@@ -0,0 +1,9 @@
+function output = directdir(obj)
+%DATADIR Get the directory where the model directly stored its output data
+%   directdir(PP)
+%
+%   See also pricklypear/datadir, pricklypear/modeldir
+
+    output = fullfile(obj.basedir, 'data', '_direct');
+
+end

matlab/@pricklypear/disp.m

@@ -0,0 +1,32 @@
+function disp(obj)
+%DISP Display information about a PricklyPear object
+%   disp(PP)
+
+    disp(' ');
+    disp('==================');
+    disp('PricklyPear object');
+    disp('==================');
+    disp(' ');
+    
+    fprintf('Experiment: %s\n', obj.sim_experiment);
+    fprintf('      Irec: %i\n', obj.sim_irec);
+    disp(' ');
+    
+    p = properties(obj);
+    pt = obj.strpad(p);
+    
+    disp('Properties');
+    disp('==========');
+    for ii = 1:numel(pt)
+        if any(strcmp(obj.template_exclude, p{ii}))
+            continue;
+        end
+        fprintf('    %s -> %g\n', pt{ii}, obj.(p{ii}));
+    end
+    disp(' ');
+    
+%     disp('Methods');
+%     disp('=======');
+%     fprintf('    run\n');
+
+end

matlab/@pricklypear/exist.m

@@ -0,0 +1,73 @@
+function output = exist(obj, varargin)
+%EXIST Check if a run exists
+%   exist(PP) checks the existence of run PP if PP has a specified
+%   experiment and irec.
+%
+%   exist(PP, 'full', true) checks if all channels are actually present. If
+%   full is false, a single channel is sufficient to return true on whether
+%   the run exists. Eventtimes and other files are never checked. Default:
+%   false
+%
+%   exist(PP, 'directMode', true) checks if there is a run in the _direct
+%   folder, meaning it has just completed by the neuron program. Default:
+%   false
+
+    % Handle input
+    p = inputParser();
+    p.addParamValue('full', false);
+    p.addParamValue('directMode', false);
+    p.parse(varargin{:});
+    p = p.Results;
+    
+    % Direct mode
+    fpath = obj.datadir();
+    irec = obj.irec;
+    if p.directMode
+        fpath = obj.directdir();
+        irec = obj.iter;
+    end
+
+    % The object must be a specific run
+    assert(~isempty(experiment(obj)), 'Please specify an experiment');
+    assert(~isempty(irec), 'Please specify a specific run');
+    
+    % Handle output
+    doesExist = 0;
+    
+    % Variables
+    flExt = obj.sim_channels;
+    
+    % If a single channel exists, consider the run as existing
+    for ii = 1:numel(flExt)
+        % Generate a filename
+        if irec == 0
+            fname = fullfile(fpath, sprintf('run_unsaved__vm_%s.txt', flExt{ii}));
+        else
+            fname = fullfile(fpath, sprintf('run_%05i__vm_%s.txt', irec, flExt{ii}));
+        end
+        
+        % Check if file exists
+        if exist(fname, 'file')
+            doesExist = doesExist + 1;
+        end
+    end
+    
+    % Handle doesExist
+    if p.full
+        doesExist = doesExist >= numel(flExt);
+    else
+        doesExist = doesExist > 1;
+    end
+    
+    % Output
+    output = doesExist;
+    if nargout == 0
+        switch doesExist
+            case 0
+                fprintf('Run %i does not exist\n', irec);
+            case 1
+                fprintf('Run %i exists\n', irec);
+        end
+    end
+
+end

matlab/@pricklypear/experiment.m

@@ -0,0 +1,7 @@
+function output = experiment(obj)
+%EXPERIMENT Get experiment of pricklypear object
+%   experiment(PP) returns the experiment of PP
+    
+    output = obj.sim_experiment;
+    
+end

matlab/@pricklypear/irec.m

@@ -0,0 +1,7 @@
+function output = irec(obj)
+%IREC Get irec of pricklypear object
+%   irec(PP) returns the irec of PP
+    
+    output = obj.sim_irec;
+    
+end

matlab/@pricklypear/modeldir.m

@@ -0,0 +1,9 @@
+function output = modeldir(obj)
+%MODELDIR Get the directory where the model is stored
+%   modeldir(PP)
+%
+%   See also pricklypear/datadir, pricklypear/directdir
+
+    output = fullfile(obj.basedir, 'neuron');
+
+end

matlab/@pricklypear/notes.m

@@ -0,0 +1,54 @@
+function output = notes(obj, varargin)
+%NOTES See/edit notes for a pricklypear experiment or run
+%   notes(PP) for notes of an entire experiment (PP without irec)
+%   notes(PP) for notes of a single run (PP with irec)
+%   notes(..., 'edit') to edit the notes for PP
+
+    % Handle input
+    p = inputParser();
+    p.addOptional('mode', 'edit', @ischar);
+    p.parse(varargin{:});
+    p = p.Results;
+    
+    % The object must be an experiment or a saved run
+    assert(~isempty(experiment(obj)), 'Please specify an experiment');
+    
+    if ~isempty(obj.irec)
+        % The run must be saved
+        assert(obj.irec > 0, 'Only saved runs can have notes');
+        
+        % The run must exist
+        assert(exist(obj) == true, 'The run was not found');
+    else
+        % The experiment must exist
+        %assert(exists(obj) == true, 'The experiment was not found');
+    end
+    
+    % File name
+    if isempty(obj.irec) % Experiment
+        fname = fullfile(datadir(obj), 'experiment_notes.txt');
+    else % Run
+        fname = fullfile(datadir(obj), sprintf('run_%05i__notes.txt', obj.irec));
+    end
+    
+    % Mode
+    switch p.mode
+        case {'view' 'read'}
+            if ~exist(fname, 'file')
+                warning('pricklypear:noNotes', 'No notes found for this experiment/run');
+                return;
+            end
+            
+            txt = fileread(fname);
+            if nargout == 0
+                fprintf('\n%s\n\n', txt);
+            else
+                output = txt;
+            end
+        case {'edit' 'write'}
+            edit(fname);
+        otherwise
+            error('The provided action was not recognized (must be view/edit)');
+    end
+
+end

matlab/@pricklypear/parent.m

@@ -0,0 +1,11 @@
+function output = parent(obj)
+%PARENT Spawn the parent of a pricklypear object
+%   P.irec -> 2
+%   PP = parent(P)
+%   PP.irec -> []
+%
+%   See also pricklypear/child
+
+    output = pricklypear(experiment(obj));
+
+end

matlab/@pricklypear/plot.m

@@ -0,0 +1,83 @@
+function plot(obj, varargin)
+%PLOT Plot the data of a PricklyPear object
+%   plot(PP)
+%
+%   See also pricklypear/anadata
+
+    % Handle input
+    p = inputParser();
+    p.addOptional('channel', '', @(x)ischar(x)||iscell(x));
+    p.addOptional('metric', '', @(x)ischar(x)||iscell(x));
+    p.parse(varargin{:});
+    p = p.Results;
+    
+    if ~isempty(p.channel) && ischar(p.channel)
+        p.channel = {p.channel};
+    end
+    if ~isempty(p.metric) && ischar(p.metric)
+        p.metric = {p.metric};
+    end
+    
+    % Plot
+    figureFullScreen(1000, 800);
+    hold on;
+
+    clrs = {'green_500' 'green_700' 'green_900'...
+            'blue_700'...
+            'green_900' 'green_700' 'green_500'...
+            'red_700' 'red_500'};
+    
+    % Select channels
+    channels = {'d0d' 'd0m' 'd0p' 'c' 'd1p' 'd1m' 'd1d' 'ap' 'ad'};
+    if ~isempty(p.channel)
+        [channels iChan] = intersect(channels, p.channel);
+        clrs = clrs(iChan);
+    end
+    
+    % Spiketimes
+    spkt = spiketimes(obj);
+    
+    % SubplotXY options
+    opts = struct();
+    opts.margin = [0.1 0.08 0.08 0.08];
+    opts.spaceBetweenPlots = 0;
+    opts.nX = 1;
+    opts.nY = numel(channels);
+    opts.iX = 1;
+    opts.iY = 1;
+    
+    % Other variables
+    h = zeros(numel(channels), 1);
+    yl = [-80 60];
+    
+    % For each channel
+    for ii = 1:numel(channels)
+        opts.iY = ii;
+        h(ii) = subplotXY(opts);
+        hold on;
+        
+        % Plot spiketimes
+        if numel(spkt) > 0
+            plot([spkt spkt zeros(size(spkt))*nan]',...
+                 [ones(size(spkt))*yl(1) ones(size(spkt))*yl(2) ones(size(spkt))*yl(2)]',...
+                 'Color', [0.75 0.75 0.75]);
+        end
+        
+        % Plot data
+        plot(obj.x, anadata(obj, channels{ii})+(ii-1), 'Color', mdc(clrs{ii}));
+        box off;
+        ylim(yl);
+        ylabel(channels{ii});
+        
+        if mod(ii,2) == 0
+            set(gca, 'Color', [1 1 1]*0.95);
+        end
+        
+        if ii < numel(channels)
+            set(gca, 'xtick', []);
+        end
+    end
+    
+    linkaxes(h, 'x');
+
+end

matlab/@pricklypear/pricklypear.m

@@ -0,0 +1,146 @@
+classdef pricklypear < handle
+    %PRICKYLPEAR Interface to the MSO neuron model from matlab
+    %   PP = pricklypear() creates a new MSO neuron model interface.
+    %   PP = pricklypear(expName) specifies the experiment name.
+    %   PP = pricklypear(expName, irec) specifies an irec.
+    %
+    %   See also pricklypear/run
+
+    properties
+        seed = 1                % RNG seed
+        iter = 99               % Used by model (do not touch!)
+
+        tstop = 2000            % Length of run
+        dt = 0.025              % Temporal resolution
+
+        v_init = -65            % Initial membrane potential
+        celsius = 38            % Temperature (celsius)
+
+        axon_soma_distance = 45 % Axon soma distance (on ipsi dendrite)
+        dend_n = 2              % Number of dendrites
+        dend_n_syn = 10         % Number of synapses per dendrite
+        dend_n_seg = 20         % Number of segments per dendrite
+        dend_length = 200       % Length of dendrites (um)
+        dend_syn_spread = 0.5   % Spread of synapses (index)
+        dend_syn_offset = 0.45  % Offset of synapses (index)
+
+        dend_0_exc_G = 11       % Excitatory conductance contra dendrite
+        dend_1_exc_G = 11       % Excitatory conductance ipsi dendrite
+
+        dend_0_inh_G = 0        % Inhibitory conductance contra dendrite
+        dend_1_inh_G = 0        % Inhibitory conductance ipsi dendrite
+
+        dend_0_exc_gain = 0.0008    % Excitatory gain contra dendrite
+        dend_1_exc_gain = 0.0008    % Excitatory gain ipsi dendrite
+
+        dend_0_inh_gain = 0.001     % Inhibitory gain contra dendrite
+        dend_1_inh_gain = 0.001     % Inhibitory gain ipsi dendrite
+
+        spk_thres_exc_contra    = 0 % Excitatory spike threshold contra
+        spk_thres_exc_ipsi      = 0	% Excitatory spike threshold ipsi
+        spk_thres_inh_contra    = 1	% Inhibitory spike threshold contra
+        spk_thres_inh_ipsi      = 1	% Inhibitory spike threshold ipsi
+
+        eventtimes_contra   = 10+(0:75:1500)	% Input spiketimes contra (ms)
+        eventtimes_ipsi     = 10+(0:100:1500)	% Input spiketimes ipsi (ms)
+        eventtimes_inh_contra   = []            % Inh input spiketimes contra (ms)
+        eventtimes_inh_ipsi     = []            % Inh input spiketimes ipsi (ms)
+        model_name          = 'default'         % Model name
+        user_data           = []                % Arbitrary user data
+        
+        % Private
+        sim_channels = {'c' 'd0p' 'd0m' 'd0d' 'd1p' 'd1m' 'd1d' 'ap' 'ad'}
+        sim_metrics = {'vm'}
+        sim_experiment = ''
+        sim_irec = []
+        template_exclude = {'eventtimes_contra' 'eventtimes_ipsi'...
+                            'eventtimes_inh_contra' 'eventtimes_inh_ipsi'...
+                            'template_exclude' 'model_name', 'user_data'...
+                            'sim_irec' 'sim_experiment' 'sim_channels' 'sim_metrics'}
+    end
+    
+    methods
+        % Class constructor
+        function obj = pricklypear(varargin)
+            % Handle input
+            p = inputParser();
+            p.addOptional('expname', '', @ischar);
+            p.addOptional('irec', []);
+            p.parse(varargin{:});
+            p = p.Results;
+            
+            % Check expname
+            assert(~strcmp(p.expname, '_direct'), '_direct is a protected experiment name, pick any other.');
+            
+            % Handle expname
+            if ~isempty(p.expname)
+                obj.sim_experiment = p.expname;
+            end
+            
+            % Handle irec
+            if ~isempty(p.expname) && ~isempty(p.irec)
+                obj.sim_irec = p.irec;
+                
+                % Check if this run exists
+                if exist(obj)
+                    % Get properties from old runs
+                    % Compatibility issues will likely arise here...
+                    props = properties(obj);
+                    if p.irec == 0
+                        storedProps = load(fullfile(obj.datadir(), sprintf('run_unsaved__properties.ppbin')), '-mat');
+                    else
+                        storedProps = load(fullfile(obj.datadir(), sprintf('run_%05i__properties.ppbin', p.irec)), '-mat');
+                    end
+                    
+                    for ii = 1:numel(props)
+                        prop = props{ii};
+                        prop_ext = prop;
+                        % Skip irec as it has already been dealt with
+                        if isequal(prop_ext, 'sim_irec')
+                            continue;
+                        end
+                        %%% Backwards compatibility
+                        % Skip a property if not stored
+                        if ~isfield(storedProps, prop)
+                            continue;
+                        end
+                        % Renamed spiketimes to eventtimes
+                        if ~isempty(strfind(prop_ext, 'spiketimes'))
+                            prop_ext = strrep(prop_ext, 'spiketimes', 'eventtimes');
+                        end
+                        % Overwrite the field
+                        obj.(prop) = storedProps.(prop_ext);
+                    end
+                else
+                    %error('Run %i does not exist', p.irec);
+                end
+                
+            end
+            
+        end
+        
+        disp(obj)
+        plot(obj, varargin)
+        obj = save(obj)
+        obj = run(obj)
+        obj = irec(obj)
+        obj = experiment(obj)
+        output = struct(obj)
+        output = exist(obj, varargin)
+        output = delete(obj, varargin)
+        output = anadata(obj, varargin)
+        output = spiketimes(obj, varargin)
+        output = notes(obj, varargin)
+        output = strpad(obj, varargin)
+        output = basedir(obj)
+        output = datadir(obj)
+        output = directdir(obj)
+        output = x(obj)
+        output = subsref(obj, irec)
+        output = parent(obj)
+        output = child(obj, irec)
+        output = runSubthresholdModel(obj, varargin)
+    end
+    
+end
+

matlab/@pricklypear/run.m

@@ -0,0 +1,91 @@
+function obj = run(obj)
+%RUN Run PricklyPear model simulation
+%   run(PP)
+
+    % Variables
+    cdir = pwd();
+    
+    % Read the base hoc file
+    fname = fullfile(obj.modeldir(), sprintf('msomodel_%s.hoc', obj.model_name));
+    fid = fopen(fname);
+    txt = textscan(fid, '%s', 'delimiter', '\n');
+    txt = txt{1};
+    fclose(fid);
+    
+    % Replace the main placeholder
+    phLine = find(strcmp('//%%%TEMPLATE%%%', txt));
+    
+    props = properties(obj);
+    phNew = cell(numel(props), 1);
+    for ii = 1:numel(props)
+        if any(strcmp(obj.template_exclude, props{ii}))
+            continue;
+        end
+        phNew{ii} = sprintf('%s = %g', props{ii}, obj.(props{ii}));
+    end
+    phNew = [{'// PricklyPear parameters'};...
+             phNew;...
+             {'// If you see this file, something went wrong during pricklypear.run(). Debug and remove this file'}];
+    
+    txt = [txt(1:phLine-1); phNew; txt(phLine+1:end)];
+    
+    % Replace the ending placeholder
+    phLine = strcmp('//%%%TEMPLATE_ENDING%%%', txt);
+    txt{phLine} = 'quit()';
+    
+    % Write the temporary hoc file
+    fname = fullfile(obj.modeldir(), 'msomodel__temp.hoc');
+    fid = fopen(fname, 'w');
+    fprintf(fid, '%s\n', txt{:});
+    fclose(fid);
+    
+    % Write the input event times
+    fname = fullfile(obj.modeldir(), 'eventtimes_contra.txt');
+    fid = fopen(fname, 'w');
+    fprintf(fid, '%i\r\n', obj.eventtimes_contra);
+    fclose(fid);
+    fname = fullfile(obj.modeldir(), 'eventtimes_ipsi.txt');
+    fid = fopen(fname, 'w');
+    fprintf(fid, '%i\n', obj.eventtimes_ipsi);
+    fclose(fid);
+    
+    % Remove existing data files
+    if exist(obj, 'directMode', true)
+        isDeleted = delete(obj, 'directMode', true);
+    end
+    parobj = parent(obj);
+    if exist(child(parobj, 0))
+        isDeleted = delete(child(parobj, 0));
+    end
+    
+    % Execute
+    cd(obj.modeldir());
+    system('neuron -nogui msomodel__temp.hoc');
+    cd(cdir);
+    
+    % Wait for simulation to finish
+    while ~exist(obj, 'directMode', true, 'full', true)
+        pause(4);
+    end
+    
+    % Move data
+    if ~isempty(obj.sim_experiment)
+        [unused unused] = mkdir(obj.datadir());
+        for ii = 1:numel(obj.sim_channels)
+            movefile(fullfile(obj.directdir(), sprintf('run_%05i__vm_%s.txt', obj.iter, obj.sim_channels{ii})),...
+                     fullfile(obj.datadir(), sprintf('run_unsaved__vm_%s.txt', obj.sim_channels{ii})));
+        end
+    end
+    copyfile(fullfile(obj.modeldir(), sprintf('eventtimes_contra.txt')),...
+             fullfile(obj.datadir(),  sprintf('run_unsaved__evt_contra.txt')));
+    copyfile(fullfile(obj.modeldir(), sprintf('eventtimes_ipsi.txt')),...
+             fullfile(obj.datadir(),  sprintf('run_unsaved__evt_ipsi.txt')));
+	props = struct(obj);
+    save(fullfile(obj.datadir(), 'run_unsaved__properties.ppbin'), '-struct', 'props');
+    
+    % Clean up
+    delete(fullfile(obj.modeldir(), 'msomodel__temp.hoc'));
+    delete(fullfile(obj.modeldir(), 'eventtimes_contra.txt'));
+    delete(fullfile(obj.modeldir(), 'eventtimes_ipsi.txt'));
+
+end

matlab/@pricklypear/runSubthresholdModel.m

@@ -0,0 +1,83 @@
+function output = runSubthresholdModel(obj, S, varargin)
+%RUNSUBTHRESHOLDMODEL Run a simulation as defined by SubthresholdModel
+%   O = runSubthresholdModel(PP, S) uses the output of subthrmodel to run 
+%   new simulations mimicking that model as best as possible. The output O
+%   contains six PricklyPear instances corresponding to the 6 DZW
+%   conditions as well as all the spiketimes for each conditions.
+%
+%   runSubthresholdModel(..., 'axon_soma_distance', 45) sets the distance
+%   between axon and soma (on ipsi dendrite) in microns. Default: 45
+
+    % TODO Nicely implement this quickfix
+    add_delay = 200; % ms
+
+    % Handle input
+    p = inputParser();
+    p.addParamValue('axon_soma_distance', 45);
+    p.parse(varargin{:});
+    p = p.Results;
+
+    % Handle output
+    output = struct;
+    output.irec = [];
+    output.spt = struct;
+    
+    % The object must be an experiment (not a run)
+    assert(~isempty(experiment(obj)), 'Please specify an experiment');
+    assert(isempty(irec(obj)), 'Do not specify a run (irec)');
+    
+    % Log
+    fprintf('Starting runSubthresholdModel...\n');
+    
+    % For each of the conditions
+    for iCond = 1:numel(S.AllCond)
+        % Get the cond
+        cond = S.AllCond{iCond};
+        
+        % Log
+        fprintf('Condition %i: %s\n', iCond, cond);
+        fprintf('  - Setting parameters...\n');
+        
+        % Set PP parameters to S parameters
+        obj.axon_soma_distance = p.axon_soma_distance;
+        obj.eventtimes_contra = S.ev.(cond).C.evt + add_delay;
+        obj.eventtimes_ipsi   = S.ev.(cond).I.evt + add_delay;
+        obj.tstop = S.Param.Dur;
+        obj.user_data = struct;
+        obj.user_data.subthr = S;
+        obj.user_data.cond = cond;
+        
+        % Log
+        fprintf('  - Parameters set!\n');
+        fprintf('  - Simulation duration: %i seconds.\n', ceil(obj.tstop/1e3));
+        fprintf('  - Expected processing time: %i seconds.\n', ceil(obj.tstop/1e3)*5);
+        fprintf('  - Running simulation...\n');
+        
+        % Run simulation
+        run(obj);
+        
+        % Log
+        fprintf('  - Run finished!\n');
+        fprintf('  - Saving run...\n');
+        
+        % Save run
+        nextIrec = save(child(obj, 0));
+        output.irec(iCond) = nextIrec;
+        
+        % Log
+        fprintf('  - Run saved as run #%i!\n', nextIrec);
+        fprintf('  - Putting spiketimes in output...\n');
+        
+        % Put spiketimes in output struct
+        output.spt.(cond) = struct;
+        output.spt.(cond).spt = spiketimes(child(obj, nextIrec));
+        
+        % Log
+        fprintf('  - Done!\n');
+        
+    end
+    
+    % Log
+    fprintf('runSubthresholdModel is done!\n');
+
+end

matlab/@pricklypear/save.m

@@ -0,0 +1,44 @@
+function output = save(obj)
+%SAVE Save a PricklyPear run
+%   save(PP(0)) stores the unsaved run from experiment(PP) to the next
+%   available irec. Storing already saved runs is not allowed.
+
+    % The object must be the unsaved simulation (irec = 0)
+    assert(~isempty(experiment(obj)), 'Please specify an experiment');
+    assert(~isempty(obj.irec) && obj.irec == 0, 'Only the unsaved simulation (irec = 0) can be saved');
+    
+    % The unsaved experiment must exist
+    assert(exist(obj) == true, 'No unsaved simulation found');
+    
+    % What is the next available irec?
+    parobj = parent(obj);
+    nextIrec = [];
+    ii = 1;
+    while isempty(nextIrec)
+        if ~exist(child(parobj, ii), 'full', true)
+            nextIrec = ii;
+        else
+            ii = ii + 1;
+        end
+    end
+    
+    % Move data
+    [unused unused] = mkdir(obj.datadir());
+    for ii = 1:numel(obj.sim_channels)
+        movefile(fullfile(obj.datadir(), sprintf('run_unsaved__vm_%s.txt', obj.sim_channels{ii})),...
+                 fullfile(obj.datadir(), sprintf('run_%05i__vm_%s.txt', nextIrec, obj.sim_channels{ii})));
+    end
+    movefile(fullfile(obj.datadir(), sprintf('run_unsaved__evt_contra.txt')),...
+             fullfile(obj.datadir(), sprintf('run_%05i__evt_contra.txt', nextIrec)));
+    movefile(fullfile(obj.datadir(), sprintf('run_unsaved__evt_ipsi.txt')),...
+             fullfile(obj.datadir(), sprintf('run_%05i__evt_ipsi.txt', nextIrec)));
+    movefile(fullfile(obj.datadir(), sprintf('run_unsaved__properties.ppbin')),...
+             fullfile(obj.datadir(), sprintf('run_%05i__properties.ppbin', nextIrec)));
+    
+    % Notify the user or output the irec
+    output = nextIrec;
+    if nargout == 0
+        fprintf('Run saved as run %i\n', nextIrec);
+    end
+    
+end

matlab/@pricklypear/spiketimes.m

@@ -0,0 +1,25 @@
+function output = spiketimes(obj, varargin)
+%SPIKETIMES Spiketimes from pricklypear simulation
+%   spiketimes(PP)
+%   spiketimes(PP, thr)
+%   spiketimes(PP, thr, chan)
+
+    % Handle input
+    p = inputParser();
+    p.addOptional('threshold', 0);
+    p.addOptional('channel', 'ap', @ischar);
+    p.parse(varargin{:});
+    p = p.Results;
+    
+    % The object must be a specific run
+    assert(~isempty(experiment(obj)), 'Please specify an experiment');
+    assert(~isempty(irec(obj)), 'Please specify a specific run');
+    
+    % Get the raw recording
+    ad = anadata(obj, p.channel, 'vm');
+    
+    % Detect peaks
+    output = peakpicker(obj.dt, ad, p.threshold);
+    output = output{1};
+
+end

matlab/@pricklypear/strpad.m

@@ -0,0 +1,60 @@
+function S = strpad(obj, S, varargin)
+%STRPAD String rightpadding
+%   S = strpad(CA) applies rightpadding of spaces to all elements of the
+%   cell array CA so that all have the same length and returns the cell
+%   array S.
+%
+%   S = strpad(CA, n) sets extra rules for the padding. If n = 0 (default),
+%   the longest string determines the length. If n > 0, n will be the
+%   length until which padding is added. Note that if a string is longer
+%   than n, it will not be truncated. If n < 0, -n is a margin, a number of
+%   spaces added to the longest string and that then decides the final
+%   length for all elements.
+%
+%   S = strpad(CA, n, '-truncate') truncates string if their length
+%   exceeds the value of n > 0 (see above).
+
+    % Handle flags
+    [flg varargin] = flagParser(varargin, 'truncate');
+
+    % Handle input
+    p = inputParser();
+    p.addRequired('S');
+    p.addOptional('n', 0);
+    p.parse(S, varargin{:});
+    p = p.Results;
+    
+    S = p.S;
+    n = p.n;
+    
+    % If S is an array, make it a cell array
+    if isnumeric(S)
+        S = cellfun(@num2str, num2cell(S), 'UniformOutput', false);
+    end
+    
+    % If S is a string, make it a cell array
+    if ischar(S)
+        S = {S};
+    end
+    
+    % Padding method depends on n
+    if n > 0
+        maxLength = p.n;
+    else
+        maxLength = max(max(cellfun(@length, S))) - n;
+    end
+    
+    % For each element, add padding
+    S = cellfun(@(s)([s char(zeros(1, maxLength-numel(s))+32)]), S, 'UniformOutput', false);
+    
+    % Truncation
+    if flg.truncate
+        S = cellfun(@(s)(s(1:maxLength)), S, 'UniformOutput', false);
+    end
+    
+    % Display mode TODO
+    if nargout > 0
+        return;
+    end
+
+end

matlab/@pricklypear/struct.m

@@ -0,0 +1,24 @@
+function output = struct(obj)
+%STRUCT Generate a struct from a PricklyPear object
+%   struct(PP)
+
+%     % Handle input
+%     p = inputParser();
+%     p.addOptional('mode', 'public', @ischar);
+%     p.parse(varargin{:});
+%     p = p.Results;
+
+    % Handle output
+    output = struct;
+    
+    % Handle the properties
+    P = properties(obj);
+    
+    % Transfer properties to output
+    for ii = 1:numel(P)
+        p = P{ii};
+        output.(p) = obj.(p);
+    end
+
+end
+

matlab/@pricklypear/subsref.m

@@ -0,0 +1,20 @@
+function output = subsref(obj, S)
+%SUBSREF Handle the subreffing of PricklyPear objects
+%   Not to be called directly
+
+    output = obj;
+
+    for ii = 1:numel(S)
+        s = S(ii);
+        
+        switch s.type
+            case '()'
+                if numel(s.subs) > 0
+                    output = obj.child(s.subs{1});
+                end
+            case '.'
+                output = builtin('subsref', output, s);
+        end
+    end
+
+end

matlab/@pricklypear/x.m

@@ -0,0 +1,7 @@
+function output = x(obj)
+%X Generate a simulation timeline (used as x axis for plotting)
+%   x(PP)
+
+    output = (0:obj.dt:obj.tstop)';
+
+end

matlab/flagParser.m

@@ -0,0 +1,31 @@
+function [S args] = flagParser(args, flags)
+%FLAGPARSER Parse varargin for flags (BEFORE INPUTPARSER)
+%   [S args] = flagParser(args, flags) parses the cell array args,
+%   checking for any flags provided in the cell array flags, then returns
+%   the args without the flags so they can be processed in inputParser and
+%   the struct S, containing all the booleans specifying the absence or
+%   presence of the flags in the args.
+%
+%   Please provide the flags without the leading dash.
+    
+    % Convert input to cell array if needed
+    if ischar(flags)
+        flags = {flags};
+    end
+    
+    % tempFlags has a dash prepended to all flags
+    tempFlags = cellfun(@(x)(['-' x]), flags, 'UniformOutput', false);
+    
+    % tempArgs contains only strings
+    tempArgs = cellfun(@num2str, args, 'UniformOutput', false);
+
+    % Check for every potential flag
+    [hasFlags idxFlag] = ismember(tempFlags, tempArgs);
+    
+    % Make the output struct S
+    S = cell2struct(num2cell(hasFlags), flags, 2);
+    
+    % Remove the flags from the args (compatible with inputParser)
+    args(idxFlag(hasFlags)) = [];
+
+end

matlab/pp_dir.template.txt

@@ -0,0 +1 @@
+C:\matlab\pricklypear\

neuron/.gitignore

@@ -0,0 +1,5 @@
+*.o
+*.c
+*.lnk
+*.tmp
+*.dll

neuron/Alpha.mod

@@ -0,0 +1,68 @@
+COMMENT
+** This mod file is based on the alpha synapse in NEURON, except for a ratio factor con in the NET_RECEIVE block.**
+ENDCOMMENT
+
+NEURON {
+	POINT_PROCESS Alpha
+	RANGE tau1, tau2, e, i,con
+	NONSPECIFIC_CURRENT i
+
+	RANGE g
+	GLOBAL total
+}
+
+UNITS {
+	(nA) = (nanoamp)
+	(mV) = (millivolt)
+	(umho) = (micromho)
+}
+
+PARAMETER {
+	tau1=.1 (ms) <1e-9,1e9>
+	tau2 = 10 (ms) <1e-9,1e9>
+	e=0	(mV)
+        con=10
+}
+
+ASSIGNED {
+	v (mV)
+	i (nA)
+	g (umho)
+	factor
+	total (umho)
+}
+
+STATE {
+	A (umho)
+	B (umho)
+}
+
+INITIAL {
+	LOCAL tp
+	total = 0
+	if (tau1/tau2 > .9999) {
+		tau1 = .9999*tau2
+	}
+	A = 0
+	B = 0
+	tp = (tau1*tau2)/(tau2 - tau1) * log(tau2/tau1)
+	factor = -exp(-tp/tau1) + exp(-tp/tau2)
+	factor = 1/factor
+}
+
+BREAKPOINT {
+	SOLVE state METHOD cnexp
+	g = B - A
+	i = g*(v - e)
+}
+
+DERIVATIVE state {
+	A' = -A/tau1
+	B' = -B/tau2
+}
+
+NET_RECEIVE(weight (umho)) {
+	state_discontinuity(A, A + weight*factor*con)
+	state_discontinuity(B, B + weight*factor*con)
+	total = total+weight*con
+}

neuron/gaussstim.mod

@@ -0,0 +1,223 @@
+COMMENT
+This mod file is based on netstim.mod in NEURON
+modified by Yi Zhou to generate Gaussian interval distributions
+ENDCOMMENT
+
+NEURON {
+	POINT_PROCESS GaussStim
+	RANGE x
+	RANGE interval, number, start, factor, refrac
+	RANGE N_backward, N_forward, N_normal, N_total
+	RANGE rnd
+	RANGE rand
+}
+
+PARAMETER {
+	interval	= 10 (ms) <1e-9,1e9> 	: ideal time between spikes (ms)
+	number	 	= 10000 <0,1e9>			: maximum number of spikes
+	start		= 10 (ms)				: delay before first spiketime
+	factor  	= 4  					: fraction of mean to compute std (= mean / factor)
+	refrac		= 0.5 <0,3>    			: absolute refractory period (upper limit of firing = 1 kHz)
+										: refrac >= dt, otherwise x can't be reset to 0 
+}
+
+ASSIGNED {
+	x			: if x > 0, cell is in refractory period
+	event (ms)	: spiketime of next event in absolute time
+	on
+	end (ms)
+	m (ms)		: mean of Gaussian
+	diff (ms)
+	N_forward	: swap spike whose value exceed forwardly one interval 
+	N_backward	: swap spike whose value exceed backwardly one interval 
+	N_normal
+	N_total
+	rnd
+	rand
+}
+
+PROCEDURE seed(x) {
+	set_seed(x)
+}
+
+INITIAL {
+	on			= 0
+	x			= 0
+	diff		= 0
+	m 			= interval / 2	: each T has normal distribution N(interval/2, interval/2/factor)
+	N_forward	= 0
+	N_backward	= 0
+	N_normal	= 0
+	N_total		= 0
+
+	if (start >= 0 && number > 0) {
+		:: generate the first spiketime
+		event = start + invl(m) 
+		
+		:: make sure first spiketime is not before 0
+		if (event < 0) {
+			event = 0
+		}
+
+		:: skip directly to next event
+		net_send(event, 3)
+	}
+}
+
+PROCEDURE init_sequence(t(ms)) {
+:: Turn on spike generation if n_spikes_max > 0
+::	- "on" is set to 1
+::	- "event" is set to "t" (input of this call)
+::	- "end" is set to "interval * n_spikes_max"
+
+	if (number > 0) {
+		on = 1
+		event = t
+		end = t + 1e-6 + interval * (number-1)
+	}
+}
+
+PROCEDURE event_time() {
+:: Generate the next spiketime
+::	- call invl which calls normrand()
+
+	:: define local variables
+	LOCAL diff2, T
+
+	:: BAD diff gets overwritten in invl :facepalm:
+	diff2 = diff
+
+	:: if n_spikes_max is larger than 0
+	if (number > 0) {
+		:: generate a new interval until next spiketime
+		T = invl(m)
+
+		:: BAD useless
+		rnd = T
+
+		:: if new spiketime T is equal to an existing spiketime
+		if (T == 0 && diff2 == 0) {
+			T = T + dt
+		}
+
+		:: compute absolute time of next event (relative to 0 ms)
+		event = event + diff2 + T
+
+		:: augment total counter
+		N_total = N_total + 1
+ 	}
+
+	:: if next event is after the end, turn off spike generation
+	if (event > end) {
+		on = 0 :: BAD not the role of this function
+	}
+}
+
+FUNCTION invl(mean (ms)) (ms) {
+:: Returns a new interval invl
+::	Warning: also sets the global variable diff
+
+	:: define local variable
+	LOCAL std
+
+	:: compute standard deviation
+	if (mean <= 0.) {
+		mean = .01 (ms) 
+	}
+	std = mean / factor  
+
+	:: generate new interval invl
+	invl = normrand(mean, std)
+
+ 	:: make sure invl is within 0 and max interval
+	if (invl >= interval) {
+		:: decrease invl by n*interval
+		invl = fmod(invl, interval)
+		:: augment forward counter
+		N_forward = N_forward + 1
+	} else if (invl < 0) {
+		:: increase invl by n*interval
+		invl = fmod(invl, interval) + interval
+		:: augment backward counter
+		N_backward = N_backward + 1
+	} else {
+		:: augment normal counter
+		N_normal = N_normal + 1
+	}
+
+	:: overwrite global variable diff
+	diff = interval - invl
+}
+
+NET_RECEIVE(w) {
+	:: EXTERNAL CALL
+	if (flag == 0) {
+		if (w > 0 && on == 0) {
+			:: turn on spike generation
+			init_sequence(t)
+			net_send(0, 1)
+		} else if (w < 0 && on == 1) {
+			:: turn off spike generation
+			on = 0
+		}
+	}
+
+	:: MAIN CALL
+	if (flag == 1 && on == 1) {
+		if (x == 0) { :: time to spike (refractory period has passed)
+			:: set x to random uniform number (useless, just needs to > 0)
+			rand = scop_random() : [0, 1] uniform
+			x = rand
+
+			:: output event to synapse
+			net_event(t)
+
+			:: generate next spiketime
+			event_time()
+
+			:: if next spiketime is between refrac and refrac+dt
+			if (event-t >= refrac && event-t <= refrac+dt) {
+				:: move next spiketime up by one dt (just after refrac)
+				event = event + dt
+			}
+
+			:: if spike generation is turned on
+			if (on == 1) {
+				:: skip directly to next event
+				net_send(event - t, 1)
+			}
+
+			:: refractory period call
+			net_send(refrac, 2)
+		} else if (x != 0) { :: time to spike (within refractory period)
+			:: output event to synapse
+			net_event(t) :: BAD even still emitted even though refractory period, should be deleted
+
+			:: generate next spiketime
+			event_time() : although this spike will be ignored, still call next spike (independent cycles)
+
+			:: if spike generation is turned on
+			if (on == 1) {
+				:: skip directly to next event
+				net_send(event - t, 1)
+			}
+		}
+	}
+
+	:: REFRACTORY PERIOD CALL (END OF REFRACTORY PERIOD)
+	if (flag == 2) {
+		x = 0
+	}
+
+	:: INITIAL CALL
+	if (flag == 3) {
+		:: if spike generation is not yet turned on
+		if (on == 0) {
+			:: turn on spike generation
+			init_sequence(t)
+
+			:: call net_receive directly
+			net_send(0, 1)
+		}
+	}
+}

neuron/ih_VCN2003.mod

@@ -0,0 +1,80 @@
+TITLE Ih channels in VCN auditory neurons
+: ih=ghmax*r*(v-eh) 
+: based on Rothman and Manis (2003c)
+: Modifications by Yi Zhou for an MSO model
+
+INDEPENDENT {t FROM 0 TO 1 WITH 1 (ms)}
+
+NEURON {
+	SUFFIX ih_VCN2003
+	USEION h READ eh WRITE ih VALENCE 1	
+	RANGE ghbar 
+	RANGE r_inf,tau_r,r_exp
+	RANGE ih,gh
+	
+}
+
+
+UNITS {
+	(mA) = (milliamp)
+	(mV) = (millivolt)
+}
+
+PARAMETER {
+	ghbar	= 0.002	(mho/cm2)  
+	eh=-43		(mV) 
+	celsius =22		(degC)
+	dt              (ms)
+	v               (mV)
+	
+}
+
+STATE {
+	r
+}
+
+ASSIGNED {
+	gh	(mho/cm2)
+	ih	(mA/cm2)
+	r_inf
+	tau_r
+	r_exp
+	tadj
+}
+
+
+BREAKPOINT {
+	SOLVE states
+	gh=ghbar *r
+	ih  = gh*(v-eh)
+}
+
+
+PROCEDURE states() {	: this discretized form is more stable
+	evaluate_fct(v)
+	r = r + r_exp * (r_inf - r)
+	VERBATIM
+	return 0;
+	ENDVERBATIM
+}
+
+UNITSOFF
+INITIAL {
+:
+:  Q10 was assumed to be 3 for both currents
+:
+	tadj = 3.0 ^ ((celsius-22)/ 10 )
+	evaluate_fct(v)
+	r= r_inf
+	}
+
+PROCEDURE evaluate_fct(v(mV)) {
+	
+	tau_r = (100000/(237*exp((v+60)/12)+17*exp(-(v+60)/14))+25)/ tadj
+	r_inf = 1/(1+exp((v+76)/7))
+   
+	r_exp = 1 - exp(-dt/tau_r)
+	
+}
+
+UNITSON

neuron/kHT_VCN2003.mod

@@ -0,0 +1,95 @@
+TITLE high threshold potassium channels in VCN auditory neurons
+  
+: k_HT=ght*(rr*n^2+(1-rr)*p)*(v-Ek)
+: based on Rothman and Manis (2003c)
+:
+: Modifications by Yi Zhou for an MSO model
+
+INDEPENDENT {t FROM 0 TO 1 WITH 1 (ms)}
+
+NEURON {
+	SUFFIX kHT_VCN2003
+	USEION k READ ek WRITE ik
+	RANGE gkbar 
+	RANGE n_inf,p_inf
+	RANGE tau_n,tau_p
+	RANGE n_exp,p_exp
+	RANGE ik,gk
+	
+}
+
+
+UNITS {
+	(mA) = (milliamp)
+	(mV) = (millivolt)
+}
+
+PARAMETER {
+	gkbar	= 0.03	(mho/cm2)  
+	ek=-70		(mV) 
+	celsius =22		(degC)
+	dt              (ms)
+	v               (mV)
+	
+}
+
+STATE {
+	n p rr
+}
+
+ASSIGNED {
+	gk(mho/cm2)
+	ik	(mA/cm2)
+	n_inf
+	p_inf
+	tau_n
+	tau_p
+	n_exp
+	p_exp
+	tadj
+}
+
+
+BREAKPOINT {
+	SOLVE states
+	gk=gkbar *(rr*n^2+(1-rr)*p)
+	ik  = gk*(v-ek)
+}
+
+
+
+PROCEDURE states() {	: this discretized form is more stable
+	evaluate_fct(v)
+	n = n + n_exp * (n_inf - n)
+	p = p + p_exp * (p_inf - p)
+	VERBATIM
+	return 0;
+	ENDVERBATIM
+}
+
+UNITSOFF
+INITIAL {
+:
+:  Q10 was assumed to be 3 for both currents
+:
+	tadj = 3.0 ^ ((celsius-22)/ 10 )
+	evaluate_fct(v)
+	n= n_inf
+	p= p_inf
+        rr=0.85
+	}
+
+PROCEDURE evaluate_fct(v(mV)) {
+	
+	tau_n = (100/(11*exp((v+60)/24)+21*exp(-(v+60)/23))+0.7)/ tadj
+	n_inf = 1/(1+exp(-(v+15)/5))^2
+        
+	tau_p = (100/(4*exp((v+60)/32)+5*exp(-(v+60)/22))+5)/ tadj
+	p_inf = 1/(1+exp(-(v+23)/6))
+	
+	n_exp = 1 - exp(-dt/tau_n)
+	p_exp = 1 - exp(-dt/tau_p)
+	
+}
+
+UNITSON

neuron/kLT_VCN2003.mod

@@ -0,0 +1,93 @@
+TITLE low threshold potassium channels in VCN auditory neurons 
+: k_LT=glt*w^4*z*(v-Ek)
+: based on Rothman and Manis (2003c)
+: Modifications by Yi Zhou for an MSO model
+
+INDEPENDENT {t FROM 0 TO 1 WITH 1 (ms)}
+
+NEURON {
+	SUFFIX kLT_VCN2003
+	USEION k READ ek WRITE ik
+	RANGE gkbar 
+	RANGE w_inf,z_inf
+	RANGE tau_w,tau_z
+	RANGE w_exp,z_exp
+	RANGE ik,gk
+	
+}
+
+
+UNITS {
+	(mA) = (milliamp)
+	(mV) = (millivolt)
+}
+
+PARAMETER {
+	gkbar	= 0.04	(mho/cm2)  
+	ek=-70		(mV) 
+	celsius =22		(degC)
+	dt              (ms)
+	v               (mV)
+	
+}
+
+STATE {
+	w z
+}
+
+ASSIGNED {
+	gk(mho/cm2)
+	ik(mA/cm2)
+	w_inf
+	z_inf
+	tau_w
+	tau_z
+	w_exp
+	z_exp
+	tadj
+}
+
+
+BREAKPOINT {
+	SOLVE states
+	gk=gkbar *w^4*z
+	ik  = gk*(v-ek)
+}
+
+
+
+PROCEDURE states() {	: this discretized form is more stable
+	evaluate_fct(v)
+	w = w + w_exp * (w_inf - w)
+	z = z + z_exp * (z_inf - z)
+	VERBATIM
+	return 0;
+	ENDVERBATIM
+}
+
+UNITSOFF
+INITIAL {
+:
+:  Q10 was assumed to be 3 for both currents
+:
+	tadj = 3.0 ^ ((celsius-22)/ 10 )
+	evaluate_fct(v)
+	w= w_inf
+	z= z_inf
+	}
+
+PROCEDURE evaluate_fct(v(mV)) {LOCAL h
+	
+	tau_w = (100/(6*exp((v+60)/6)+16*exp(-(v+60)/45))+1.5)/ tadj
+	w_inf = 1/(1+exp(-(v+48)/6))^0.25
+        
+	tau_z = (1000/(exp((v+60)/20)+exp(-(v+60)/8))+50)/ tadj
+	h=0.5
+	z_inf = (1-h)/(1+exp((v+71)/10))+h
+	
+	w_exp = 1 - exp(-dt/tau_w)
+	z_exp = 1 - exp(-dt/tau_z)
+	
+}
+
+UNITSON

neuron/mosinit.hoc

@@ -0,0 +1,34 @@
+load_file("nrngui.hoc")
+
+strdef tstr
+
+xpanel("Figures for Zhou et al 2005")
+xlabel("Press a button to create the figure")
+xradiobutton("Main simulation", "restart(\"MSO_Zhouetal_2005\")")
+xradiobutton("Additional simulation", "restart(\"test_inputs\")")
+xlabel("After you click on one of the above buttons you will be asked")
+xlabel("to enter a number which will be used to seed the random number")
+xlabel("generator.  After you have entered the number please press the")
+xlabel("\"Single Run\" button which will appear in a Control window")
+xlabel("Note: the main simulation takes 3 1/4 minutes on a 2Ghz machine")
+xlabel("and the additional simulation takes less than 15 seconds.")
+xpanel(5,100)
+
+pwmcnt = PWManager[0].count // the initial gui should not be dismissed
+
+proc restart() {local i
+	objref graphItem, save_window_
+	for i=0, n_graph_lists-1 {
+		graphList[i].remove_all()
+	}
+	flush_list.remove_all()
+	fast_flush_list.remove_all()
+	doNotify()
+	for (i= PWManager[0].count-1; i >= pwmcnt; i -= 1) {
+		PWManager[0].close(i)
+		doNotify()
+	}
+	
+	sprint(tstr, "%s.hoc", $s1)
+	load_file(1, tstr)
+}

neuron/msomodel_base.hoc

@@ -0,0 +1,378 @@
+// Init measurements
+objref vm_s, vm_d0p, vm_d0m, vm_d0d, vm_d1p, vm_d1m, vm_d1d, vm_ap, vm_ad
+vm_s = new Vector(tstop/dt+2,0)
+vm_d0p = new Vector(tstop/dt+2,0)
+vm_d0m = new Vector(tstop/dt+2,0)
+vm_d0d = new Vector(tstop/dt+2,0)
+vm_d1p = new Vector(tstop/dt+2,0)
+vm_d1m = new Vector(tstop/dt+2,0)
+vm_d1d = new Vector(tstop/dt+2,0)
+vm_ap = new Vector(tstop/dt+2,0)
+vm_ad = new Vector(tstop/dt+2,0)
+vm_s.record(&maincell.soma.v(0.5), dt)
+vm_d0p.record(&maincell.dend[0].v(0.1), dt)
+vm_d0m.record(&maincell.dend[0].v(0.5), dt)
+vm_d0d.record(&maincell.dend[0].v(0.9), dt)
+vm_d1p.record(&maincell.dend[1].v(0.1), dt)
+vm_d1m.record(&maincell.dend[1].v(0.5), dt)
+vm_d1d.record(&maincell.dend[1].v(0.9), dt)
+vm_ap.record(&maincell.axon.v(0.1), dt)
+vm_ad.record(&maincell.axon.v(0.5), dt)
+
+
+// Variables left over from obsolete functions
+// TODO find obsolete variables
+bin = 0.1 	 // [ms]  binwidth of APs
+cvode.active(0) // constant integration time step
+
+access maincell.axon
+
+objref apc, v_ISI, v_voltage
+objref pre_E, pre_I_contra, pre_I_ipsi
+
+v_ISI = new Vector()
+v_voltage = new Vector(tstop/dt+2,0)
+
+pre_E = new Vector()
+pre_I_contra = new Vector()
+pre_I_ipsi = new Vector()
+
+//==== counting APs
+AP_threshold = -10 	   //[mV] AP detection threshold
+apc = new APCount(0.5)   // measure action potenials at the midpoint of the axon
+apc.thresh = AP_threshold
+
+//==== saving event times and membrane potentials
+apc.record(v_ISI) //record event time
+v_voltage.record(&maincell.axon.v(0.5), dt)
+
+//==== monitoring input events
+netcon0[0].record(pre_E)		//record input event time
+netcon_inhi0[0].record(pre_I_contra)  //record input event time
+netcon_inhi1[0].record(pre_I_ipsi)    //record input event time
+
+objref v_save,v_all,v_psth
+objref g_psth,x_psth
+objref d, xtime
+
+bin_psth=bin
+
+objref v_pth, g_pth, x_pth
+
+objref hist_ISI,xtime,g_ISI,hist
+objref mean_T, sq_T, num_T
+objref temp,t1, v1,v2   // operating vectors
+
+mean_ISI=0
+OVF=0    // overflow
+EI=0
+
+mean_ISI2=0
+std_ISI=0
+sq_sum=0
+CV_ISI=0
+
+objref v_AP, mean_T, sq_T, num_T     		// rates and spike intervals at one ITD over N_run
+objref r_mean,std_mean,vs_ITD, phase_ITD    	// rates at all ITDs
+objref T_ITD,std_ITD,CV_ITD   	     	// spike interval statistics at all ITDs
+objref ITD,g_Rate,g_VS,data
+
+N_run = n_runs_total
+ITD_max=1 //[ms] maximum ITD tested
+
+
+// Functions
+proc reset_membrane() {
+	finitialize(v_init)
+	fcurrent()
+	
+	maincell.soma {
+		gl_na_VCN2003
+		el_na_VCN2003=(ina+ik+ih+gl_na_VCN2003*v_init)/gl_na_VCN2003
+		g_rest1=gl_na_VCN2003+gk_kLT_VCN2003+gk_kHT_VCN2003+gna_na_VCN2003+gh_ih_VCN2003
+		cm/(g_rest1*1000)
+	}
+	maincell.axon {
+		el_na_VCN2003=(ina+ik+ih+gl_na_VCN2003*v_init)/gl_na_VCN2003
+		g_rest3=gl_na_VCN2003+gk_kLT_VCN2003+gk_kHT_VCN2003+gna_na_VCN2003+gh_ih_VCN2003
+		cm/(g_rest3*1000)
+	}
+	maincell.dend[0] {
+		e_pas=v_init
+		g_rest2=g_pas
+		cm/(g_rest2*1000)
+	}
+	maincell.dend[1] {
+		e_pas=v_init
+	}
+
+	fcurrent()
+}
+// proc reset_membrane() {
+// 	finitialize(v_init)
+// 	fcurrent()  // make sure that leak current balance the non-leak current
+// 	maincell.soma {
+// 		print "gl_na_VCN2003=[S/cm^2]=", gl_na_VCN2003
+// 		print "soma el=", el_na_VCN2003=(ina+ik+ih+gl_na_VCN2003*v_init)/gl_na_VCN2003
+// 		print "g_rest=S/cm2  ", g_rest1=gl_na_VCN2003+gk_kLT_VCN2003+gk_kHT_VCN2003+gna_na_VCN2003+gh_ih_VCN2003
+// 		print "tau_rest=", cm/(g_rest1*1000)
+// 	}
+// 	maincell.axon {
+// 		print "axon el=", el_na_VCN2003=(ina+ik+ih+gl_na_VCN2003*v_init)/gl_na_VCN2003
+// 		print "g_rest=S/cm2  ", g_rest3=gl_na_VCN2003+gk_kLT_VCN2003+gk_kHT_VCN2003+gna_na_VCN2003+gh_ih_VCN2003
+// 		print "tau_rest=", cm/(g_rest3*1000)
+// 	}
+// 	maincell.dend[0] {
+// 		print "dendrite el=", e_pas=v_init
+// 		print " dend g_rest=S/cm2  ", g_rest2=g_pas
+// 		print "tau_rest=", cm/(g_rest2*1000)
+// 	}
+// 	maincell.dend[1] {
+// 		e_pas=v_init
+// 	}
+
+// 	fcurrent()
+// }
+
+reset_membrane()
+
+proc reset() {
+	reset_membrane()
+
+	//==== Reset all vec and var
+	N_run = n_runs_total
+
+	v_voltage=new Vector(tstop/dt+2,0)
+	v_voltage.record(&maincell.axon.v(0.5),dt)   // measure the midpoint of the axon
+	
+	vm_s = new Vector(tstop/dt+2,0)
+	vm_d0p = new Vector(tstop/dt+2,0)
+	vm_d0m = new Vector(tstop/dt+2,0)
+	vm_d0d = new Vector(tstop/dt+2,0)
+	vm_d1p = new Vector(tstop/dt+2,0)
+	vm_d1m = new Vector(tstop/dt+2,0)
+	vm_d1d = new Vector(tstop/dt+2,0)
+	vm_ap = new Vector(tstop/dt+2,0)
+	vm_ad = new Vector(tstop/dt+2,0)
+	vm_s.record(&maincell.soma.v(0.5), dt)
+	vm_d0p.record(&maincell.dend[0].v(0.1), dt)
+	vm_d0m.record(&maincell.dend[0].v(0.5), dt)
+	vm_d0d.record(&maincell.dend[0].v(0.9), dt)
+	vm_d1p.record(&maincell.dend[1].v(0.1), dt)
+	vm_d1m.record(&maincell.dend[1].v(0.5), dt)
+	vm_d1d.record(&maincell.dend[1].v(0.9), dt)
+	vm_ap.record(&maincell.axon.v(0.1), dt)
+	vm_ad.record(&maincell.axon.v(0.5), dt)
+
+  	v_AP=new Vector(N_run,0)
+  	mean_T=new Vector(N_run,0)
+  	sq_T=new Vector(N_run,0)
+  	num_T=new Vector(N_run,0)
+}
+
+objref fl
+objref flPG
+objref spkt_contra
+objref spkt_ipsi
+strdef flname
+strdef cmd
+fl = new File()
+flPG = new File()
+
+ispkt_contra = 0
+ispkt_ipsi   = 0
+spkt_contra  = new Vector()
+spkt_ipsi    = new Vector()
+
+proc rerun() {
+	reset()
+
+	//*********************************
+	while (N_run>0) {
+		v_ISI=new Vector()
+		v_voltage=new Vector(tstop/dt+2,0)
+		apc.record(v_ISI) 				//record event time
+		v_voltage.record(&maincell.axon.v(0.5),dt)    // measure APs at the midpoint of the axon
+
+		vm_s = new Vector(tstop/dt+2,0)
+		vm_d0p = new Vector(tstop/dt+2,0)
+		vm_d0m = new Vector(tstop/dt+2,0)
+		vm_d0d = new Vector(tstop/dt+2,0)
+		vm_d1p = new Vector(tstop/dt+2,0)
+		vm_d1m = new Vector(tstop/dt+2,0)
+		vm_d1d = new Vector(tstop/dt+2,0)
+		vm_ap = new Vector(tstop/dt+2,0)
+		vm_ad = new Vector(tstop/dt+2,0)
+		vm_s.record(&maincell.soma.v(0.5), dt)
+		vm_d0p.record(&maincell.dend[0].v(0.1), dt)
+		vm_d0m.record(&maincell.dend[0].v(0.5), dt)
+		vm_d0d.record(&maincell.dend[0].v(0.9), dt)
+		vm_d1p.record(&maincell.dend[1].v(0.1), dt)
+		vm_d1m.record(&maincell.dend[1].v(0.5), dt)
+		vm_d1d.record(&maincell.dend[1].v(0.9), dt)
+		vm_ap.record(&maincell.axon.v(0.1), dt)
+		vm_ad.record(&maincell.axon.v(0.5), dt)
+		
+		flPG.ropen("eventtimes_contra.txt")
+		spkt_contra.scanf(flPG)
+		flPG.close()
+		flPG.ropen("eventtimes_ipsi.txt")
+		spkt_ipsi.scanf(flPG)
+		flPG.close()
+
+		for i = 0,9 {
+			gen0[i].timeToNextEvt = spkt_contra.x[0] - dt
+			gen1[i].timeToNextEvt = spkt_ipsi.x[0] - dt
+		}
+
+		finitialize(v_init)
+		fcurrent()
+
+		integrate()
+		// get_pth()
+		// get_ISI()
+		v_AP.x[N_run-1] = apc.n
+		N_run = N_run-1
+		
+		// Export to temporary files
+		sprint(flname, "../data/_direct/run_%05i__vm_c.txt.tmp", iter)
+		fl.wopen(flname)
+		vm_s.printf(fl)
+		fl.close()
+		
+		sprint(flname, "../data/_direct/run_%05i__vm_d0p.txt.tmp", iter)
+		fl.wopen(flname)
+		vm_d0p.printf(fl)
+		fl.close()
+		
+		sprint(flname, "../data/_direct/run_%05i__vm_d0m.txt.tmp", iter)
+		fl.wopen(flname)
+		vm_d0m.printf(fl)
+		fl.close()
+		
+		sprint(flname, "../data/_direct/run_%05i__vm_d0d.txt.tmp", iter)
+		fl.wopen(flname)
+		vm_d0d.printf(fl)
+		fl.close()
+		
+		sprint(flname, "../data/_direct/run_%05i__vm_d1p.txt.tmp", iter)
+		fl.wopen(flname)
+		vm_d1p.printf(fl)
+		fl.close()
+		
+		sprint(flname, "../data/_direct/run_%05i__vm_d1m.txt.tmp", iter)
+		fl.wopen(flname)
+		vm_d1m.printf(fl)
+		fl.close()
+		
+		sprint(flname, "../data/_direct/run_%05i__vm_d1d.txt.tmp", iter)
+		fl.wopen(flname)
+		vm_d1d.printf(fl)
+		fl.close()
+		
+		sprint(flname, "../data/_direct/run_%05i__vm_ap.txt.tmp", iter)
+		fl.wopen(flname)
+		vm_ap.printf(fl)
+		fl.close()
+		
+		sprint(flname, "../data/_direct/run_%05i__vm_ad.txt.tmp", iter)
+		fl.wopen(flname)
+		vm_ad.printf(fl)
+		fl.close()
+
+		// Rename files from temporary name to destined name
+		sprint(cmd, "rename ..\\data\\_direct\\run_%05i__vm_c.txt.tmp run_%05i__vm_c.txt", iter, iter)
+		print "Cmd: ",cmd
+		system(cmd)
+		sprint(cmd, "rename ..\\data\\_direct\\run_%05i__vm_d0p.txt.tmp run_%05i__vm_d0p.txt", iter, iter)
+		print "Cmd: ",cmd
+		system(cmd)
+		sprint(cmd, "rename ..\\data\\_direct\\run_%05i__vm_d0m.txt.tmp run_%05i__vm_d0m.txt", iter, iter)
+		system(cmd)
+		sprint(cmd, "rename ..\\data\\_direct\\run_%05i__vm_d0d.txt.tmp run_%05i__vm_d0d.txt", iter, iter)
+		system(cmd)
+		sprint(cmd, "rename ..\\data\\_direct\\run_%05i__vm_d1p.txt.tmp run_%05i__vm_d1p.txt", iter, iter)
+		system(cmd)
+		sprint(cmd, "rename ..\\data\\_direct\\run_%05i__vm_d1m.txt.tmp run_%05i__vm_d1m.txt", iter, iter)
+		system(cmd)
+		sprint(cmd, "rename ..\\data\\_direct\\run_%05i__vm_d1d.txt.tmp run_%05i__vm_d1d.txt", iter, iter)
+		system(cmd)
+		sprint(cmd, "rename ..\\data\\_direct\\run_%05i__vm_ap.txt.tmp run_%05i__vm_ap.txt", iter, iter)
+		system(cmd)
+		sprint(cmd, "rename ..\\data\\_direct\\run_%05i__vm_ad.txt.tmp run_%05i__vm_ad.txt", iter, iter)
+		system(cmd)
+	}
+
+	N_run = n_runs_total
+
+}
+
+proc integrate() {
+	// while the simulation has not yet ended
+	while (t < tstop) {
+		// CONTRA
+		// if t is earlier than the last spiketime
+		if (t < spkt_contra.x[spkt_contra.size()-1]) {
+			// if t is right after a spiketime has been issued, issue the next one
+			if (t >= spkt_contra.x[ispkt_contra-1]+dt && t < spkt_contra.x[ispkt_contra-1]+dt*2) {
+				if (ispkt_contra == spkt_contra.size()-1) {
+					for i = 0,9 {
+						gen0[i].timeToNextEvt = 0
+					}
+				} else {
+					for i = 0,9 {
+						gen0[i].timeToNextEvt = spkt_contra.x[ispkt_contra+1] - spkt_contra.x[ispkt_contra]
+					}
+				}
+
+				ispkt_contra = ispkt_contra + 1
+			} else if (t >= dt && t < dt*2) { // if t is right the beginning of the simulation
+				for i = 0,9 {
+					gen0[i].timeToNextEvt = spkt_contra.x[ispkt_contra+1] - spkt_contra.x[ispkt_contra] - dt
+				}
+				ispkt_contra = ispkt_contra + 1
+			}
+		}
+
+		// IPSI
+		// if t is earlier than the last spiketime
+		if (t < spkt_ipsi.x[spkt_ipsi.size()-1]) {
+			// if t is right after a spiketime has been issued, issue the next one
+			if (t >= spkt_ipsi.x[ispkt_ipsi-1]+dt && t < spkt_ipsi.x[ispkt_ipsi-1]+dt*2) {
+				if (ispkt_ipsi == spkt_ipsi.size()-1) {
+					for i = 0,9 {
+						gen1[i].timeToNextEvt = 0
+					}
+				} else {
+					for i = 0,9 {
+						gen1[i].timeToNextEvt = spkt_ipsi.x[ispkt_ipsi+1] - spkt_ipsi.x[ispkt_ipsi]
+					}
+				}
+
+				ispkt_ipsi = ispkt_ipsi + 1
+			} else if (t >= dt && t < dt*2) { // if t is right the beginning of the simulation
+				for i = 0,9 {
+					gen1[i].timeToNextEvt = spkt_ipsi.x[ispkt_ipsi+1] - spkt_ipsi.x[ispkt_ipsi] - dt
+				}
+				ispkt_ipsi = ispkt_ipsi + 1
+			}
+		}
+
+		fadvance() // advance solution
+	}
+	print "N_run=",N_run   // show run time
+	print "spike number=",apc.n
+	// print "-------------------------"
+	//print "Number_input_E=",pre_E.size
+	//print "Number_input_I_contra=",pre_I_contra.size
+	//print "Number_input_I_ipsi=",pre_I_ipsi.size
+	//print "*************"
+
+	if(v_ISI.size<=2) { v_ISI=new Vector(3,0) }
+	L_ISI=v_ISI.size
+	//print "L_ISI=",L_ISI
+	L_PTH=v_voltage.size
+	//print "L_PTH=",L_PTH
+}
+
+rerun()

neuron/msomodel_default.hoc

@@ -0,0 +1,249 @@
+//----------
+// VARIABLES
+//----------
+
+seed_x = 1
+
+iter = 1
+
+n_runs_total = 1
+tstop = 2000
+dt = 0.025
+
+v_init = -65
+celsius = 38
+
+axon_soma_distance = 45
+dend_n = 2
+dend_n_syn = 10
+dend_n_seg = 20
+dend_length = 200
+dend_syn_spread = 0.5
+dend_syn_offset = 0.45
+
+dend_0_exc_G = 11
+dend_1_exc_G = 11
+
+dend_0_inh_G = 0
+dend_1_inh_G = 0
+
+dend_0_exc_gain = 0.0008
+dend_1_exc_gain = 0.0008
+
+dend_0_inh_gain = 0.001
+dend_1_inh_gain = 0.001
+
+// Spike rate / threshold
+spk_thres_exc_contra = 0
+spk_thres_exc_ipsi   = 0
+spk_thres_inh_contra = 1
+spk_thres_inh_ipsi   = 1
+
+// Overrule variables using templating
+//%%%TEMPLATE%%%
+
+
+//----------------
+// CELL DEFINITION
+//----------------
+
+begintemplate  Cell
+ 	public soma, axon, dend
+	create soma, axon, dend[1]
+
+  	proc init() {
+        ndend = $1
+        dend_length = $2
+        axon_soma_distance = $3
+		create soma, axon, dend[ndend]
+
+		soma {
+			nseg = 1
+			diam = 20  // um
+			L    = 40  // um
+			Ra   = 200 // ohm.cm
+			cm   = 1   // uF/cm2
+
+			insert kHT_VCN2003
+			insert kLT_VCN2003
+			insert na_VCN2003
+			insert ih_VCN2003
+
+			ek  = -70 // mV
+			ena = 55  // mV
+			eh  = -43 // mV
+
+			gkbar_kHT_VCN2003 = 0     // S/cm2
+			gkbar_kLT_VCN2003 = 0     // S/cm2
+			gnabar_na_VCN2003 = 0.1   // S/cm2
+			ghbar_ih_VCN2003  = 0     // S/cm2
+			gl_na_VCN2003     = 0.002 // S/cm2
+		}
+
+		axon {
+			nseg = 51
+			diam = 2   // um
+			L    = 400 // um^2
+			Ra   = 200 // ohm.cm
+			cm   = 1   // uF/cm2
+
+			insert kHT_VCN2003
+			insert kLT_VCN2003
+			insert na_VCN2003
+			insert ih_VCN2003
+			
+			ek  = -70 // mV
+			ena = 55  // mV
+			eh  = -43 // mV
+
+			// total G[nS]=gkbar*area*10
+
+			gkbar_kHT_VCN2003 = 0.02   // S/cm2
+			gkbar_kLT_VCN2003 = 0.03   // S/cm2
+			gnabar_na_VCN2003 = 0.3    // S/cm2
+			ghbar_ih_VCN2003  = 0.0015 // S/cm2
+			gl_na_VCN2003     = 0.002  // S/cm2
+    	}
+
+		for i = 0, ndend-1 dend[i] {
+        	nseg = 20
+        	diam = 3   // um
+        	L    = 200 // um
+        	Ra   = 200 // ohm.cm
+			cm   = 1   // uF/cm2
+
+        	insert pas
+			e_pas = -65
+			g_pas = 0.002
+    	}
+
+     	connect dend[0](0), soma(0)
+     	connect dend[1](0), soma(1)
+		connect axon(0), dend[1](axon_soma_distance/dend_length)
+	}
+endtemplate Cell
+
+
+//--------------------
+// CELL INITIALISATION
+//--------------------
+
+objectvar maincell
+maincell = new Cell(dend_n, dend_length, axon_soma_distance)
+access maincell.soma
+
+
+//--------------------
+// EXCITATORY SYNAPSES
+//--------------------
+
+// Variables
+objectvar syn0[10], gen0[10], syn1[10], gen1[10]
+
+// Event generators
+for i = 0,9 {
+	// Contra
+	gen0[i] = new RelayStim(0.5)
+	gen0[i].seed(seed_x+i)
+
+	// Ipsi
+	gen1[i] = new RelayStim(0.5)
+	gen1[i].seed(seed_x+10+i)
+}
+
+// Contra synapses
+maincell.dend[0] {
+	for i = 0,9 {
+		syn0[i] = new Alpha(i/((dend_n_syn-1)/dend_syn_spread)+dend_syn_offset)
+
+		syn0[i].tau1 = 0.1	// ms
+		syn0[i].tau2 = 0.1	// ms
+		syn0[i].con  = dend_0_exc_G	// nS-mho
+		syn0[i].e    = 0		// mV
+	}
+}
+
+// Ipsi synapses
+maincell.dend[1] {
+	for i = 0,9 {
+		syn1[i] = new Alpha(i/((dend_n_syn-1)/dend_syn_spread)+dend_syn_offset)
+
+		syn1[i].tau1 = 0.1	// ms
+		syn1[i].tau2 = 0.1	// ms
+		syn1[i].con  = dend_1_exc_G	// nS-mho
+		syn1[i].e    = 0		// mV
+	}
+}
+
+// Connect generators and synapses
+objref netcon0[10], netcon1[10]
+e_delay = 0 // Synaptic delay (ms)
+
+for i = 0,9  {
+	netcon0[i] = new NetCon(gen0[i], syn0[i], spk_thres_exc_contra, e_delay, dend_0_exc_gain)
+	netcon1[i] = new NetCon(gen1[i], syn1[i], spk_thres_exc_ipsi,   e_delay, dend_1_exc_gain)
+}
+
+
+//--------------------
+// INHIBITORY SYNAPSES
+//--------------------
+
+// Variables
+objectvar syn_inhi0[10], gen_inhi0[10], syn_inhi1[10], gen_inhi1[10]
+
+// Event generators
+for i = 0,9 {
+	gen_inhi0[i]=new GaussStim(0.5)  //  contralateral inputs
+	gen_inhi0[i].interval=2
+	gen_inhi0[i].start=0
+	gen_inhi0[i].number=10000
+	gen_inhi0[i].factor=10
+	gen_inhi0[i].seed(seed_x+30+i)
+
+	gen_inhi1[i]=new GaussStim(0.5)  //  ipsilateral inputs
+	gen_inhi1[i].interval=2
+	gen_inhi1[i].start=0
+	gen_inhi1[i].number=10000
+	gen_inhi1[i].factor=10
+	gen_inhi1[i].seed(seed_x+50+i)
+}
+
+// Somatic synapses
+maincell.soma {
+	for i = 0,9 {
+		syn_inhi0[i]      = new Alpha(0.5)
+		syn_inhi0[i].tau1 = 0.1 //ms
+		syn_inhi0[i].tau2 = 2   //ms
+		syn_inhi0[i].con  = dend_0_inh_G //umho
+		syn_inhi0[i].e    = -70 // mV
+	}
+
+	for i = 0,9 {
+		syn_inhi1[i]      = new Alpha(0.5)
+		syn_inhi1[i].tau1 = 0.1 //ms
+		syn_inhi1[i].tau2 = 2   //ms
+		syn_inhi1[i].con  = dend_1_inh_G //umho
+		syn_inhi1[i].e    = -70 // mV
+	}
+}
+
+// Connect generators and synapses
+objref netcon_inhi0[10],netcon_inhi1[10]
+I_delay0 = 0 // Synaptic delay
+I_delay1 = 0 // Synaptic delay
+
+for i = 0,9{
+	netcon_inhi0[i] = new NetCon(gen_inhi0[i], syn_inhi0[i], spk_thres_inh_contra, I_delay0, dend_0_inh_gain)
+	netcon_inhi1[i] = new NetCon(gen_inhi1[i], syn_inhi1[i], spk_thres_inh_ipsi,   I_delay1, dend_1_inh_gain)
+}
+
+
+//---------------------
+// LOAD SIMULATION CODE
+//---------------------
+
+load_file("msomodel_base.hoc")
+
+// Ending template (useful to quit automated runs)
+//%%%TEMPLATE_ENDING%%%

neuron/msomodel_passive_axon.hoc

@@ -0,0 +1,236 @@
+//----------
+// VARIABLES
+//----------
+
+seed_x = 1
+
+iter = 1
+
+n_runs_total = 1
+tstop = 2000
+dt = 0.025
+
+v_init = -65
+celsius = 38
+
+axon_soma_distance = 45
+dend_n = 2
+dend_n_syn = 10
+dend_n_seg = 20
+dend_length = 200
+dend_syn_spread = 0.5
+dend_syn_offset = 0.45
+
+dend_0_exc_G = 11
+dend_1_exc_G = 11
+
+dend_0_inh_G = 0
+dend_1_inh_G = 0
+
+dend_0_exc_gain = 0.0008
+dend_1_exc_gain = 0.0008
+
+dend_0_inh_gain = 0.001
+dend_1_inh_gain = 0.001
+
+// Spike rate / threshold
+spk_thres_exc_contra = 0
+spk_thres_exc_ipsi   = 0
+spk_thres_inh_contra = 1
+spk_thres_inh_ipsi   = 1
+
+// Overrule variables using templating
+//%%%TEMPLATE%%%
+
+
+//----------------
+// CELL DEFINITION
+//----------------
+
+begintemplate  Cell
+ 	public soma, axon, dend
+	create soma, axon, dend[1]
+
+  	proc init() {
+        ndend = $1
+        dend_length = $2
+        axon_soma_distance = $3
+		create soma, axon, dend[ndend]
+
+		soma {
+			nseg = 1
+			diam = 20  // um
+			L    = 40  // um
+			Ra   = 200 // ohm.cm
+			cm   = 1   // uF/cm2
+
+			insert kHT_VCN2003
+			insert kLT_VCN2003
+			insert na_VCN2003
+			insert ih_VCN2003
+
+			ek  = -70 // mV
+			ena = 55  // mV
+			eh  = -43 // mV
+
+			gkbar_kHT_VCN2003 = 0     // S/cm2
+			gkbar_kLT_VCN2003 = 0     // S/cm2
+			gnabar_na_VCN2003 = 0.1   // S/cm2
+			ghbar_ih_VCN2003  = 0     // S/cm2
+			gl_na_VCN2003     = 0.002 // S/cm2
+		}
+
+		axon {
+			nseg = 51
+			diam = 2   // um
+			L    = 400 // um^2
+			Ra   = 200 // ohm.cm
+			cm   = 1   // uF/cm2
+
+        	insert pas
+			e_pas = -65
+			g_pas = 0.002
+    	}
+
+		for i = 0, ndend-1 dend[i] {
+        	nseg = 20
+        	diam = 3   // um
+        	L    = 200 // um
+        	Ra   = 200 // ohm.cm
+			cm   = 1   // uF/cm2
+
+        	insert pas
+			e_pas = -65
+			g_pas = 0.002
+    	}
+
+     	connect dend[0](0), soma(0)
+     	connect dend[1](0), soma(1)
+		connect axon(0), dend[1](axon_soma_distance/dend_length)
+	}
+endtemplate Cell
+
+
+//--------------------
+// CELL INITIALISATION
+//--------------------
+
+objectvar maincell
+maincell = new Cell(dend_n, dend_length, axon_soma_distance)
+access maincell.soma
+
+
+//--------------------
+// EXCITATORY SYNAPSES
+//--------------------
+
+// Variables
+objectvar syn0[10], gen0[10], syn1[10], gen1[10]
+
+// Event generators
+for i = 0,9 {
+	// Contra
+	gen0[i] = new RelayStim(0.5)
+	gen0[i].seed(seed_x+i)
+
+	// Ipsi
+	gen1[i] = new RelayStim(0.5)
+	gen1[i].seed(seed_x+10+i)
+}
+
+// Contra synapses
+maincell.dend[0] {
+	for i = 0,9 {
+		syn0[i] = new Alpha(i/((dend_n_syn-1)/dend_syn_spread)+dend_syn_offset)
+
+		syn0[i].tau1 = 0.1	// ms
+		syn0[i].tau2 = 0.1	// ms
+		syn0[i].con  = dend_0_exc_G	// nS-mho
+		syn0[i].e    = 0		// mV
+	}
+}
+
+// Ipsi synapses
+maincell.dend[1] {
+	for i = 0,9 {
+		syn1[i] = new Alpha(i/((dend_n_syn-1)/dend_syn_spread)+dend_syn_offset)
+
+		syn1[i].tau1 = 0.1	// ms
+		syn1[i].tau2 = 0.1	// ms
+		syn1[i].con  = dend_1_exc_G	// nS-mho
+		syn1[i].e    = 0		// mV
+	}
+}
+
+// Connect generators and synapses
+objref netcon0[10], netcon1[10]
+e_delay = 0 // Synaptic delay (ms)
+
+for i = 0,9  {
+	netcon0[i] = new NetCon(gen0[i], syn0[i], spk_thres_exc_contra, e_delay, dend_0_exc_gain)
+	netcon1[i] = new NetCon(gen1[i], syn1[i], spk_thres_exc_ipsi,   e_delay, dend_1_exc_gain)
+}
+
+
+//--------------------
+// INHIBITORY SYNAPSES
+//--------------------
+
+// Variables
+objectvar syn_inhi0[10], gen_inhi0[10], syn_inhi1[10], gen_inhi1[10]
+
+// Event generators
+for i = 0,9 {
+	gen_inhi0[i]=new GaussStim(0.5)  //  contralateral inputs
+	gen_inhi0[i].interval=2
+	gen_inhi0[i].start=0
+	gen_inhi0[i].number=10000
+	gen_inhi0[i].factor=10
+	gen_inhi0[i].seed(seed_x+30+i)
+
+	gen_inhi1[i]=new GaussStim(0.5)  //  ipsilateral inputs
+	gen_inhi1[i].interval=2
+	gen_inhi1[i].start=0
+	gen_inhi1[i].number=10000
+	gen_inhi1[i].factor=10
+	gen_inhi1[i].seed(seed_x+50+i)
+}
+
+// Somatic synapses
+maincell.soma {
+	for i = 0,9 {
+		syn_inhi0[i]      = new Alpha(0.5)
+		syn_inhi0[i].tau1 = 0.1 //ms
+		syn_inhi0[i].tau2 = 2   //ms
+		syn_inhi0[i].con  = dend_0_inh_G //umho
+		syn_inhi0[i].e    = -70 // mV
+	}
+
+	for i = 0,9 {
+		syn_inhi1[i]      = new Alpha(0.5)
+		syn_inhi1[i].tau1 = 0.1 //ms
+		syn_inhi1[i].tau2 = 2   //ms
+		syn_inhi1[i].con  = dend_1_inh_G //umho
+		syn_inhi1[i].e    = -70 // mV
+	}
+}
+
+// Connect generators and synapses
+objref netcon_inhi0[10],netcon_inhi1[10]
+I_delay0 = 0 // Synaptic delay
+I_delay1 = 0 // Synaptic delay
+
+for i = 0,9{
+	netcon_inhi0[i] = new NetCon(gen_inhi0[i], syn_inhi0[i], spk_thres_inh_contra, I_delay0, dend_0_inh_gain)
+	netcon_inhi1[i] = new NetCon(gen_inhi1[i], syn_inhi1[i], spk_thres_inh_ipsi,   I_delay1, dend_1_inh_gain)
+}
+
+
+//---------------------
+// LOAD SIMULATION CODE
+//---------------------
+
+load_file("msomodel_base.hoc")
+
+// Ending template (useful to quit automated runs)
+//%%%TEMPLATE_ENDING%%%

neuron/na_VCN2003.mod

@@ -0,0 +1,107 @@
+TITLE Na channels in VCN auditory neurons of guinea pig
+ 
+: na=gna*m^3*h   
+: based on Rothman and Manis 2003
+: Modifications by Yi Zhou for an MSO model
+
+
+INDEPENDENT {t FROM 0 TO 1 WITH 1 (ms)}
+
+NEURON {
+	SUFFIX na_VCN2003
+	USEION na READ ena WRITE ina
+	NONSPECIFIC_CURRENT il
+	RANGE gnabar 
+	RANGE m_inf,h_inf
+	RANGE tau_m,tau_h
+	RANGE m_exp,h_exp
+	RANGE ina,gna
+	RANGE gl,el
+	
+}
+
+
+UNITS {
+	(mA) = (milliamp)
+	(mV) = (millivolt)
+}
+
+PARAMETER {
+	gnabar	= 0.2	(mho/cm2) 
+	ena=55		(mV) 
+	gl= 4.0e-4	(mho/cm2)   
+        el=-57		(mV)
+	celsius=22 		(degC)
+	dt              (ms)
+	v               (mV)
+	
+}
+
+STATE {
+	m h
+}
+
+ASSIGNED {
+	gna (mho/cm2)
+	ina	(mA/cm2)
+	il	(mA/cm2)
+	m_inf
+	h_inf
+	tau_m
+	tau_h
+	m_exp
+	h_exp
+	tadj3
+	
+}
+
+
+BREAKPOINT {
+	SOLVE states
+	gna=gnabar * m*m*h
+	ina  = gna * (v - ena)
+	il=gl*(v-el)
+}
+
+
+
+PROCEDURE states() {	: this discretized form is more stable
+	evaluate_fct(v)
+	m = m + m_exp * (m_inf - m)
+	h = h + h_exp * (h_inf - h)
+	VERBATIM
+	return 0;
+	ENDVERBATIM
+}
+
+UNITSOFF
+INITIAL {
+:
+:  Q10 was assumed to be 3 for both currents
+:
+	tadj3 = 3.0 ^ ((celsius-22)/ 10 )
+	
+	evaluate_fct(v)
+	m= m_inf
+	h= h_inf
+	}
+
+PROCEDURE evaluate_fct(v(mV)) { 
+	
+	m_inf = 1 / (1+exp(-(v + 38) / 7))
+    	h_inf = 1 / (1+exp((v + 65) / 6))
+
+    	tau_m =  (10 / (5*exp((v+60) / 18) + 36*exp(-(v+60) / 25))) + 0.04
+    	tau_h =  (100 / (7*exp((v+60) / 11) + 10*exp(-(v+60) / 25))) + 0.6
+	
+	tau_m=tau_m/tadj3	
+	tau_h=tau_h/tadj3
+
+	m_exp = 1 - exp(-dt/tau_m)
+	h_exp = 1 - exp(-dt/tau_h)
+	
+}
+
+UNITSON
+
+

neuron/readme.txt

@@ -0,0 +1,30 @@
+This model package provides NEURON codes for a multi-compartment model described 
+in Zhou Y, Carney LH, Colburn HS (2005). The model simulates recent physiological 
+results from the gerbil medial superior olive (MSO) that reveal that blocking 
+glycinergic inhibition can shift the tuning for the interaural time difference (ITD)
+of the cell (Brand et al., 2002). 
+
+The model has a bipolar structure and incorporates two anatomic observations 
+in the MSO: (1) the axon arises from the dendrite that receives ipsilateral 
+inputs (Smith, 1995); and (2) inhibitory synapses are located primarily 
+on the soma in adult animals. Specific Hodgkin-Huxley type of ionic channels 
+are inserted into the axon to generate action potentials: fast sodium (Ina), 
+low threshold potassium (ILTK), high threshold potassium (IHTK), 
+hyperpolarization-active inward current (Ih), and the leakage channel (Il). 
+
+The kinetics of all channels are based on the model of Type II cells in the 
+anteroventral cochlear nuclei (Rothman and Manis, 2003c). The soma has only 
+fast sodium channels. Ionic channels and the driving function for synaptic 
+inputs are described in .mod files. The hoc files describe simulation procedures 
+for modeling neural responses to ITDs (MSO_Zhouetal_2005.hoc) and a test kit 
+for exploring input characteristics (test_inputs.hoc). Use NEURON command nrnivmodl 
+to compile *.mod files, and then execute a *.hoc file. Please refer to NEURON 
+documentations for further instructions on compiling mod files and 
+launching hoc files on different operation systems.  
+
+The file MSO_Zhouetal_2005.hoc produces the rate-ITD responses 
+in Figures. 5, 6, 8, 11 in Zhou Y, Carney LH, Colburn HS (2005). 
+
+Go to http://www.bu.edu/dbin/binaural/MSOdemo for more instructions on using the model. 
+
+Any questions for implementing the model should be directed to Yi Zhou (yizhou@bu.edu). 

neuron/relaystim.mod

@@ -0,0 +1,75 @@
+COMMENT
+This mod file is based on netstim.mod in NEURON
+modified by Yarmo Mackenbach to relay spike times provided by the hoc script
+ENDCOMMENT
+
+NEURON {
+	POINT_PROCESS RelayStim
+	RANGE x, timeToNextEvt, refrac
+}
+
+PARAMETER {
+	timeToNextEvt	= 10 (ms) <1e-9,1e9> 	: time until next spiketime (ms)
+	refrac			= 0.1 (ms) <0,3>    	: absolute refractory period (ms)
+}
+
+ASSIGNED {
+	x	: unknown significance
+	on	: is spike generation active?
+}
+
+PROCEDURE seed(x) {
+	set_seed(x)
+}
+
+INITIAL {
+	: init variables
+	x	= 0
+	on	= 0
+
+	: wait until first event
+	net_send(timeToNextEvt, 3)
+}
+
+NET_RECEIVE(w) {
+	: MAIN CALL
+	if (flag == 1 && on == 1) {
+		: because somehow this is a crucial step
+		x = 1
+
+		: output event to synapse
+		net_event(t)
+
+		: disable spike generation if timeToNextEvt is null
+		if (timeToNextEvt == 0) {
+			on = 0
+		}
+
+		: if spike generation is turned on
+		if (on == 1) {
+			: wait until next event
+			net_send(timeToNextEvt, 1)
+		}
+
+		: refractory period call (somehow required)
+		net_send(refrac, 2)
+	}
+
+	: REFRACTORY PERIOD CALL (END OF REFRACTORY PERIOD)
+	if (flag == 2) {
+		: unknown signifance (but a required step)
+		x = 0
+	}
+
+	: INITIAL CALL
+	if (flag == 3) {
+		: if spike generation is not yet turned on
+		if (on == 0) {
+			: turn on spike generation
+			on = 1
+
+			: call net_receive directly
+			net_send(0, 1)
+		}
+	}
+}

neuron/test_inputs.hoc

@@ -0,0 +1,601 @@
+// *****************  Test Inputs  **********************
+
+// ------------  Introduction   --------------------------------  
+// 01/01/2005
+// This program generats periodic input spikes that simulate afferent inputs to cells in the Medial Superio Olive (MSO) 
+// of the mammalian brainstem as described in 
+// Zhou, Carney, and Colburn (2005) J. Neuroscience, 25(12):3046-3058   
+
+
+
+
+
+// ******* Codes start here ****************************************
+
+seed_x=xred("Please choose a random seed number <1,1000>", 1, 1, 1000) //string,default,min,max
+
+strdef inputDir,outputDir
+outputDir="D:\msomodel\"
+inputDir="D:\msomodel\"
+
+  
+  load_file("nrngui.hoc")
+ 
+  create cell
+  
+  access cell
+
+
+
+
+
+  //**** Excitatory and inhibitory synapses **************** 
+  //---------------------------------------------------------
+  
+  //----     parameter description     ----------------------
+  //syn0[]: E alpha synapses on contralateral dendrite, dend0
+  //syn1[]: E alpha synapses on ipsilateral dendrite, dend1
+  //gen0[]: E stimuli to syn0
+  //gen1[]: E stimuli to syn1
+  //netcon0[]: input- E synapse connection on dend0
+  //netcon1[]: input- E synapse connection on dend1
+
+  //syn_inhi0[]:  I alpha synapse on the soma
+  //gen_inhi0[]:  I stimuli to syn_inhi0
+  //netcon_inhi0[]: input- I synapse connection on the soma
+  
+ 
+  // ************ Excitatory synapses on two dendrites  ********
+  objectvar syn0[10], gen0[10],syn1[10], gen1[10]
+  
+     //----- add spike, virtual spikes
+     for i=0, 9{
+                gen0[i] = new GaussStim(0.5)
+		gen0[i].interval=2   // [ms] input period
+  		gen0[i].start=0	  //  [ms] arrival time of input spikes
+		gen0[i].number=10000
+		gen0[i].factor=20	   // phase-locking factor F
+		
+		gen1[i] = new GaussStim(0.5) 
+		gen1[i].interval=2   
+  		gen1[i].start=0	 
+  		gen1[i].number=10000
+		gen1[i].factor=20	   
+		
+		gen0[i].seed(seed_x+i)
+                gen1[i].seed(seed_x+10+i)
+  	} 
+ 
+ 	//----- add EE expsyn at dend0
+  	cell {  //contra side
+		for i=0,9 { 
+		n=20  //n=nseg  doesn't change after nseg
+		syn0[i] = new Alpha(i*1/n+1/n/2)  // scatter along the firsthalf dend
+
+		syn0[i].tau1 =0.1  //ms    
+		syn0[i].tau2 =0.1   //ms
+		syn0[i].con=10    //nS- mho
+		syn0[i].e =0   // mV
+			}                            
+   		}
+
+  	cell { //ipsi side
+  		for i=0,9 {
+		n=20  //n=nseg  
+		syn1[i] = new Alpha(i*1/n+1/n/2)
+		
+		syn1[i].tau1 =0.1  //ms    may add some jitters
+		syn1[i].tau2 =0.1   //ms
+		syn1[i].con=0    //nS-mho
+		syn1[i].e =0   // mV
+    			}     
+		}
+  
+   
+   //---  Connect gen[] with syn[] 
+   	objref netcon0[10],netcon1[10]
+
+    	Rave_E_contra=240 // spikes/sec
+	Rave_E_ipsi=240 
+
+    	x_thres_contra=1-Rave_E_contra*gen0[0].interval/1000        //   prob(x>1-p)=p
+    	x_thres_ipsi=1-Rave_E_ipsi*gen1[0].interval/1000        
+      
+       if (x_thres_contra<=0) {x_thres_contra=1e-4} 
+       if (x_thres_ipsi<=0) {x_thres_ipsi=1e-4} 
+
+        
+    	e_delay=0  // no synaptic delay
+    
+    	for i=0,9  { 
+    		netcon0[i]=new NetCon(gen0[i], syn0[i],x_thres_contra,e_delay,0.001) 
+    										//thres delay gain[uS]
+		netcon1[i]=new NetCon(gen1[i], syn1[i],x_thres_ipsi,e_delay,0.001) 
+										//thres delay gain[uS]
+       	}
+  
+ 
+
+  
+  v_init=-65
+  
+
+    
+ //----  parameter and data vector descriptions     ------------ 
+
+  // v_ISI: [ms]  AP event time
+  // v_voltage:  occurrence of events at each dt
+ 
+  
+  tstop=1000   	// [ms]  stimulus duration
+  dt=0.025   	// [ms]  simulation time step
+  bin=0.1 	// [ms]  binwidth of APs  
+  
+  cvode.active(0)   // constant integration time step
+
+
+  
+
+  objref v_ISI,v_voltage
+  objref pre_E
+  
+  v_ISI=new Vector()
+  v_voltage=new Vector(tstop/dt+2,0)
+
+  
+  
+  //==== for input ISIs
+  netcon0[0].record(v_ISI)//record input event time
+  v_voltage.record(&gen0[0].x,dt) //record the event
+  AP_threshold=x_thres_contra  // reset if record inputs
+  
+ 
+  
+  //****  function  "pth": period time histogram  **** 
+  //-----------------------------------------------
+  //----  data vector description     ------------
+  // v_save: bin v_voltage for each run
+  // v_pth: wrap v_all with length=gen.interval/bin and normalize
+  // vs_real: vector strength
+  // phase_real: mean phase 
+
+  // get_pth(): get v_save from each run
+  
+  
+  objref v_save, v_all, v_pth, g_pth, x_pth, v_stand0
+  
+  proc get_pth() { 
+		v_save=new Vector(L_PTH)
+  		v_save.spikebin(v_voltage,AP_threshold) //(v_source, threshold)
+		v_save.rebin(v_save,bin/dt)  // rebin the time axis
+		v_all.add(v_save)   // add all trials
+  	}
+  	
+  proc pth() { 
+  		x=0
+  		y=0
+  		vs_real=0
+  		angel_real=0
+		length=gen0[0].interval/bin  //length of period points
+		
+		N_fold=v_all.size/length
+		v_pth= new Vector(length,0)
+		x_pth= new Vector(length,0)  // phase vector
+		x_pth.indgen(1/length)
+		
+		//==== wrapping over one period
+		for i=0,N_fold-1 { 
+			for j=0,length-1 {
+					v_pth.x[j]=v_pth.x[j]+v_all.x[length*i+j]
+							}
+			}
+		v_pth.div(v_pth.sum)
+
+		v_stand0=new Vector(x_pth.size,0)     // dend[0] input PSH
+		v_stand0.indgen(1/length)
+  		v_stand0.apply("Gauss0")
+		v_stand0.rotate(gen0[0].start/bin)
+		
+		//==== calculate the VS and mean phase of responses
+		for i=0,length-1 {
+			x=x+v_pth.x[i]*cos(2*PI*i/length)
+			y=y+v_pth.x[i]*sin(2*PI*i/length)
+		}
+		vs_real=sqrt(x^2+y^2)     //compute VS
+		print "VS=",vs_real
+		
+		angel_real=atan(y/x)
+			if(x<0 ) { angel_real=angel_real+PI
+					} else if (y<0 && x>0) {
+						angel_real=angel_real+2*PI
+						}	
+		print "Phase=",angel_real/(2*PI)
+		
+		
+  		g_pth.size(0,1,0,0.3)
+  		g_pth.align(0.5,0.2)
+  		//g_pth.label(0.9,0.1,"Phase")
+		g_pth.color(1)
+  		g_pth.label(0.3,0.9,"Period Histogram")
+  		g_pth.begin() 
+		v_pth.mark(g_pth,x_pth,"+",9,1)
+		v_pth.line(g_pth,x_pth,1,2)
+		g_pth.color(2)
+		v_stand0.line(g_pth,x_pth,2,2)
+		
+  				  		
+}
+
+  proc VS() {
+		vs=exp(-PI^2/(2*gen0[0].factor^2))
+	}
+  
+  func Gauss0() { //for PST
+        sigma2=0.5/gen0[0].factor
+        return exp(-($1-0.5)^2/(2*sigma2^2))/((sqrt(2*PI)*sigma2)*length)
+}	
+  
+  
+  //-------------- END pth () -------------------------
+  //---------------------------------------------------
+
+
+
+
+
+  //****  function  "ISI": interspike interval  **** 
+  //-----------------------------------------------
+  //----  data vector description     ------------
+  
+  // v_ISI: event times
+  // bin_ISI : bin width of ISI	
+  // hist_ISI: histogram of ISI for all run
+  // hist: histogram of ISI for each run
+  // mean_T, sq_T, num_T: ISI statistics for each run
+
+  // get_ISI(): get data for each run
+  
+  objref hist_ISI,xtime,g_ISI,hist, v_stand
+  objref mean_T, sq_T, num_T
+  objref temp,t1, v1,v2   // operating vectors 
+  
+  mean_ISI=0
+  OVF=0    // overflow
+  EI=0
+  
+  mean_ISI2=0
+  std_ISI=0
+  sq_sum=0
+  CV_ISI=0
+  
+  proc get_ISI() {
+  		
+  		v2=new Vector(L_ISI,0)
+  		v1=new Vector(L_ISI,0)
+  		
+  		
+  		v1.add(v_ISI)
+  		v2.add(v_ISI)
+  		v1.rotate(1) // right shift
+  		v2.sub(v1)
+  		v2.rotate(-1)
+  		v2.resize(v2.size-1)  // cut off the initial interval
+  		print "mean ISI interval   ", v2.mean
+  		print "ISI std    ",v2.stdev
+  		 
+  		temp = new Vector(v2.size,0)
+  		temp.add(v2)
+  		temp.sub(v2.mean)  // for calculate the square sum
+  		
+  		mean_T.x[N_run-1]=v2.mean
+  		num_T.x[N_run-1]=v2.size
+  		sq_T.x[N_run-1]=temp.sumsq()
+  		
+  		mean_ISI=v2.mean + mean_ISI  // mean ISI for each run
+  
+  		print "**********************************"
+  		
+  		hist=v2.histogram(0, topbin+OVFbin, bin_ISI) 
+		
+  		if(hist.size-hist_ISI.size==1) { hist.resize(hist_ISI.size)} 
+  		hist_ISI.add(hist)
+		
+  	}	
+  	
+  proc ISI() {
+  		hist_ISI.div(hist_ISI.sum)  //normalize
+  		OVF=hist_ISI.sum(topbin/bin_ISI, (topbin+OVFbin)/bin_ISI-1)
+  		EI=hist_ISI.sum(0,1.5*T/bin_ISI-1)   // the first whole interval
+  		
+		hist_ISI.resize(topbin/bin_ISI+1)
+		
+  		xtime=new Vector(hist_ISI.size,0)
+  		xtime.indgen(bin_ISI)
+
+		v_stand=new Vector(hist_ISI.size,0)
+  		v_stand.indgen(bin_ISI)
+  		v_stand.apply("Gauss")
+  		
+  		g_ISI.size(0,topbin,0,0.2)
+  		g_ISI.align(0.5,0.5)
+  		//g_ISI.label(0.9,0.2,"T")
+  		g_ISI.label(0.2,0.9,"ISI")
+  		g_ISI.beginline()  
+  		hist_ISI.mark(g_ISI,xtime,"o",6,1)
+  		hist_ISI.line(g_ISI,xtime,1,2)
+		v_stand.line(g_ISI,xtime,2,2)
+		
+  		
+  		//==== calculate the mean and the std of intervals  
+   		t1=new Vector(N_run,0)
+   		t1.add(mean_T)
+   		t1.mul(num_T)
+   		
+   
+   		t2=t1.sum()
+   		mean_ISI2=t2/num_T.sum()  //overall ISI mean 
+  		
+  		t1=new Vector(N_run,0)
+  		t1.add(mean_T)
+  		t1.sub(mean_ISI2)
+  		t1.mul(t1)
+  		t1.mul(num_T)
+  		sq_sum=sq_T.sum()+t1.sum()
+  		
+  		if(num_T.sum()>1) {
+  			std_ISI=sqrt(sq_sum/(num_T.sum-1))   //Overall ISI std
+  		} else { std_ISI=sqrt(sq_sum) } 
+  		
+  		if(mean_ISI2>0) {
+  				CV_ISI=std_ISI/mean_ISI2   //coefficient of variation
+  			} else{ CV_ISI=1
+  		}            
+  		  
+   }
+
+  func Gauss() { //for ISI
+        sigma=sqrt(2)/2*(gen0[0].interval/gen0[0].factor)  //std=sqrt(2)*T/2/factor
+        return exp(-($1-gen0[0].interval)^2/(2*sigma^2))*bin/(sqrt(2*PI)*sigma)
+}
+  
+
+  //----------- END ISI() -------------------------
+  //-----------------------------------------------
+
+
+  N_run=10   // ten trials
+ 
+
+   //----- reset function for all simulations
+   //------------------------------------------
+  
+  proc reset() {
+  	
+       
+        //==== Reset synapse, and inputs
+        
+	getcon()
+	getstim()	
+	changelevel()
+	
+
+	//==== Reset all vec and var 
+	N_run=10
+         
+        v_voltage=new Vector(tstop/dt+2,0)
+        v_voltage.record(&gen0[0].x,dt) //record the event
+    
+        mean_T=new Vector(N_run,0)
+  	sq_T=new Vector(N_run,0)
+  	num_T=new Vector(N_run,0)
+
+	T=gen0[0].interval
+	bin_ISI=bin
+  	topbin=int(10*T/bin_ISI)*bin_ISI  // 10 intervals and 1 for OVF
+  	OVFbin=int(1*T/bin_ISI)*bin_ISI
+  	  	
+  	mean_ISI2=0
+   	std_ISI=0
+   	sq_sum=0
+   	CV_ISI=0
+   	
+  	mean_ISI=0
+  	OVF=0
+  	EI=0
+  	
+  	hist_ISI=new Vector((topbin+OVFbin)/bin_ISI+1,0)
+	v_all=new Vector(tstop/dt+2,0)
+  	v_all.rebin(v_all,bin/dt)
+  
+  
+  }
+  
+  //---------- end reset() ----------------
+  //---------------------------------------
+  
+
+
+
+
+  
+  //----- Main simulation function --------------------
+  // Measure model response at one ITD condition  
+  //---------------------------------------------------
+
+
+  proc rerun() {
+        
+	 reset()
+
+  	//*********************************
+            while (N_run>0) {
+                v_ISI=new Vector()
+  		v_voltage=new Vector(tstop/dt+2,0)
+
+            	  //==== for input ISI and PS
+  		netcon0[0].record(v_ISI)//record input event time
+  		v_voltage.record(&gen0[0].x,dt) //record the event
+  		  		
+  		finitialize(v_init) 
+  		fcurrent()
+
+  		integrate()
+  		get_pth()
+  		get_ISI()
+               
+		N_run=N_run-1 	
+	  	} 
+
+   
+   N_run=10
+   
+   g_ISI.erase_all()
+   g_pth.erase_all()
+   ISI()
+   pth()
+  
+   print "----------- END -----------------------"
+}
+  
+  
+ proc integrate() {
+    
+    	while (t < tstop) {
+    		fadvance()  // advance solution 
+ 		 }
+  print "N_run=",N_run   // show run time 
+  print "-------------------------"
+  print "Number_input_E=",v_ISI.size
+  
+  
+  if(v_ISI.size<=2) { v_ISI=new Vector(3,0) }
+  L_ISI=v_ISI.size
+  //print "L_ISI=",L_ISI
+  L_PTH=v_voltage.size		
+  //print "L_PTH=",L_PTH	
+}
+
+
+
+
+//******* update synapse and input information   ****
+
+
+  proc changelevel() {
+	
+	x_thres_contra=1-Rave_E_contra*gen0[0].interval/1000        //   prob(x>1-p)=p
+    	x_thres_ipsi=1-Rave_E_ipsi*gen1[0].interval/1000        
+
+	
+      
+	if (x_thres_contra<=0) {x_thres_contra=1e-4} 
+        if (x_thres_ipsi<=0) {x_thres_ipsi=1e-4} 
+        
+
+  	for i=0,9  { 
+    	netcon0[i].threshold=x_thres_contra  
+	netcon1[i].threshold=x_thres_ipsi
+       }
+      
+  }
+
+  proc getcon() {
+		for i=0,9 {
+		syn0[i].con=syn0[0].con    //umho
+		syn1[i].con=syn1[0].con    //umho
+		
+		}
+		print "*******************"
+		print "syn0 con =",syn0[1].con
+		print "syn1 con =",syn1[1].con
+		
+		print "--------------------------"
+  }
+
+  proc getstim() {
+		//------- EE -------------
+		for i=0, 9{
+                gen0[i].interval=gen0[0].interval   // 
+  		gen0[i].start=gen0[0].start	  //  identify ITD
+		gen0[i].factor=gen0[0].factor	
+				
+		gen1[i].interval=gen1[0].interval   
+  		gen1[i].start=gen1[0].start	  //  identify ITD
+		gen1[i].factor=gen1[0].factor		
+  	        } 
+
+		
+       		
+  	        print "E0 start =",gen0[1].start
+		print "E0 interval =",gen0[1].interval
+		print "E1 start =",gen1[1].start
+		print "E1 interval =",gen1[1].interval  
+		
+		print "*******************"
+				
+ }
+  //----------- END rerun() -----------------------------------
+  //-----------------------------------------------------------
+
+
+
+  //**** GUI  "synapse and stimuli"  ****
+  //---------------------------------------------
+  
+  objref syn_con, gen_start
+  
+  syn_con = new HBox()
+  syn_con.intercept(1)       
+  
+  xpanel("")
+  xlabel("E0 syn contra")
+  xvalue("G [nS]","syn0[0].con")
+  xpanel()
+  
+  
+  syn_con.intercept(0)       
+  syn_con.map("SYN",0,150,-1,50)              
+
+
+ //--------------------------------------------------------------------------
+
+  gen_start = new HBox()
+  gen_start.intercept(1)       
+  
+  xpanel("")
+  xlabel("E0 input contra")
+  xvalue("delay [ms]","gen0[0].start")
+  xvalue("period","gen0[0].interval")
+  xvalue("synchrony factor","gen0[0].factor")
+  xvalue("rate [spikes/sec]","Rave_E_contra")
+  xpanel()
+  
+  
+ 
+  gen_start.intercept(0)       
+  gen_start.map("GEN",0,350,-1,50)              
+
+  
+//--------------------------------------------------------------------------
+
+
+  objref boxv
+  boxv=new VBox()
+  boxv.intercept(1)
+  xpanel("")
+  xbutton("calculate VS of E0 input","VS()")
+  xvalue("vs")
+  xbutton("Reset","reset()")
+  xbutton("Single Run", "rerun()")
+  xbutton("Quit", "quit()")
+  
+  g_ISI=new Graph()
+  g_pth=new Graph()
+  xpanel()
+  boxv.intercept(0)
+  boxv.map("CONTROL",600,50,-1,50)
+  //----------------------------------------------------------------
+
+
+
+