% This script experiments with subtle patterns driving Hodgkin-Huxley 
% neurons more accurately via dynamic effects

clear all

global HH_T HH_I

dt = .0002;
T = .6;
psc = PSC(dt);

% two PSC time constants will be compared 
pte = 0:dt:.2;
pte2 = 0:dt:.05;
psce = exp(-pte/.02); psce = psce/sum(psce);
psce2 = exp(-pte2/.005); psce2 = psce2/sum(psce2);

trials = 100;
for i = 1:trials
    % Assume a pool of 50,000 neurons each of which fires at a low rate over
    % the course of 1/2 second, except for one 5ms bin in which it fires at a
    % higher rate. There is an even distribution of such bins over the .5s. We
    % choose the 10,000 of these neurons with elevated bins in the 200-300 ms
    % range (again evenly distributed). Poisson firing all around. Groups of
    % 500 for each of the 20 bins from 200-300ms. 
    lowRate = 200 * .1;
    highRate = 200 * .1799; %876 trials needed with elevated rate 200*0.1799 (error: 12.2/20 spikes/s) (from exp_power results)
    n = 10000;
    groupSize = 500;
    
    % generate firing patterns for lead-in of .1s plus .5s of simulation 
    [spikes, cov] = genUncorrelated(n, T, .0002, lowRate, [1 0 0], struct('RT', 0));
    
    % generate additional elevated-rate spikes and add to above patterns at appropriate times  
    [elevSpikes, cov] = genUncorrelated(n, .05, .0002, highRate-lowRate, [1 0 0], struct('RT', 0));
    add = T - .3 + .005 * floor((0:n-1)/groupSize)' * ones(1,size(elevSpikes,2));
    elevSpikes = elevSpikes + add;
    elevSpikes(find(elevSpikes-add > .005)) = 0;
    spikes = [spikes elevSpikes]; %ordering doesn't matter here

    spikeIndices = round(spikes/dt);
    len = round(T/dt);
    spikeSignal = zeros(1,len);
    for j = 1:size(spikes,1)
        ind = spikeIndices(j,:);
        ind = ind(find(ind>0 & ind<=len));
        spikeSignal(ind) = spikeSignal(ind) + 1;
    end
    
    % two alternative excitatory inputs calculated from different PSC models 
    excitation = conv(spikeSignal, psce); excitation = excitation(1:len);
    excitation = excitation * 2;
    excitation = max(74, excitation); %clipped to reduce startup artefact

    excitation2 = conv(spikeSignal, psce2); excitation2 = excitation2(1:len);        
    excitation2 = excitation2 * 2;
    excitation2 = max(74, excitation2);

    HH_T = (dt*1000) * (1:length(excitation));
    
    % we test a range of balancing inhibition and report the best
    % results (i.e. the ones that result in the most selective firing)
    inhibition = 70:1:95;

    for j = 1:length(inhibition)
        % Hodgkin-Huxley simulations 
        HH_I = excitation - inhibition(j);        
        [odet1,r] = ode45('hh', [1 T*1000], [0 0 0 0]);
        V1 = r(:,1)';
        indices = find([0 diff(V1 > 50)] > .5);
        spikeTimes = odet1(indices) / 1000;
	
        HH_I = excitation2 - inhibition(j);
        [odet2,r] = ode45('hh', [1 T*1000], [0 0 0 0]);
        V2 = r(:,1)';
        indices = find([0 diff(V2 > 50)] > .5);
        spikeTimes2 = odet2(indices) / 1000;
        
        layer1(i,j,1:size(spikes,2)) = spikes(j,:); %example presynaptic (near-random) neuron
        layer2a(i,j,1:length(spikeTimes)) = spikeTimes; %with slow-PSC excitation 
        layer2b(i,j,1:length(spikeTimes2))= spikeTimes2; %with fast-PSC excitation 
	end
    
    save data_subtleDrive.mat    
end