% Simplified version of synfire chain in cortical network à la Abeles, Bienenstock, Diesmann, etc.  


clear all
close all

P = 10;                         % length of network
N = 10;                         % width of network
epoch = .11;
bin_length = .0001;
V_rest = -.07;
V_th = -.06;
V_start = -.061                 % initial potential
tau_m = 0.01;
input_strength = .003;          % external input to set chain in motion
input_jitter = .001;
input_length = .003;            % duration of external input
excitatory_noise = .0002;       % strength of excitatory background noise
inhibitory_noise = .000105;     % strength of inhibitory background noise
excitatory_jitter = .00003;     % variability of excitatory background noise
inhibitory_jitter = .00003;     % variability of inhibitory background noise
conduction_time = .01;
conduction_jitter = .1;
convergence_jitter = .000002;
amplification_factor = .4;      % amplification of convergent input

time_axis = 0:bin_length:epoch;
bin_number = length(time_axis);
input_onset = round(bin_number/10);
input_bins = round(input_length/bin_length);            % duration of external input (bin units)
conduction_bins = round(conduction_time/bin_length);    % conduction time (bin units)
input_offset = input_onset + input_bins;

latency = [];                       % latency of first spike in pool
mean_latency = [];
standard_deviation = [];

crossco_time = [-epoch:bin_length:epoch]; 

p = 1;     % first pool (receives external input)
V = V_start*ones(N,bin_number);
excitatory_background_noise = max(0,normrnd(excitatory_noise,excitatory_jitter,N,bin_number));
inhibitory_background_noise = max(0,normrnd(inhibitory_noise,inhibitory_jitter,N,bin_number));
V = V + excitatory_background_noise - inhibitory_background_noise;


for n = 1:N    % loop for all neurons in first pool

    occurence_spike = 0;

    for bin = 2:bin_number
        
        if bin >= input_onset & bin <= input_offset
            
            input = max(0,normrnd(input_strength,input_jitter,1,1));
            
        else
            
            input = 0;
            
        end;
        
        if V(n,bin-1) < V_th
            
            V_derivative = (V_rest - V(n,bin-1))/tau_m;
            V(n,bin) = V(n,bin-1) + bin_length*V_derivative + input + excitatory_background_noise(n,bin) - inhibitory_background_noise(n,bin);   % membrane potential
            
        else
            
            V(n,bin-1) = .03;     % EPSP
            V(n,bin) = V_rest;
            
            if occurence_spike == 0
                latency = [latency bin];
            end
            
            occurence_spike = 1;
            
         end
            
     end
     
     figure(1)
     subplot(P,1,p), plot(time_axis,V), axis([0 epoch -.1 .1]), title('Membrane Potentials of Neurons in Various Pool Levels'), ylabel('V_m (V)')     % membrane potential over time (first pool)
     
end  

V_1 = V(1,:);
V_2 = V(2,:);
for bin = 1:bin_number
    if V_1(bin) == .03
        V_1(bin) = 1;
    else 
        V_1(bin) = 0;   
    end
    if V_2(bin) == .03
        V_2(bin) = 1;
    else 
        V_2(bin) = 0;
    end
end
crossco = xcov(V_1,V_2);
figure(2)
subplot(P,1,p), plot(crossco_time,crossco), axis([-epoch epoch -.5 Inf]), title('Cross Correlograms')   % cross correlogram of two first neurons of first pool
    
mean_latency = [mean_latency mean(latency)];
standard_deviation = [standard_deviation std(latency)];
latency =[];

level = input_offset;   % beginning of second pool

for p = 2:P     % internally driven pools
    
    clear V_1, clear V_2
    
    potential_cache(:,:) = V(:,:);     % store potential matrix (neurons, bins)

    excitatory_background_noise = max(0,normrnd(excitatory_noise,excitatory_jitter,N,bin_number));
    inhibitory_background_noise = max(0,normrnd(inhibitory_noise,inhibitory_jitter,N,bin_number));
    
    V = V_start*ones(N,bin_number);
    V = V + excitatory_background_noise - inhibitory_background_noise;
   
    input = zeros(N,bin_number);

    for n = 1:N
        
        occurence_spike = 0;
        
        for bin = level:bin_number
        
            potential_transform(n,bin) = potential_cache(n,bin-conduction_bins+round(normrnd(0,conduction_jitter,1,1)));     % introduce conduction delays and conduction jitter
            
        end
        
        for bin = level:bin_number
        
            individual_input(n,bin) = max(0,potential_transform(n,bin));
    
        end
    
        input = sum(individual_input,1);                        % add all imputs from previous pool
        amplified_input = input*amplification_factor;           % adjust (amplify) inputs (e.g., to keep model stable)
        input_rep = repmat(amplified_input,N,1);
        input_rep_var = normrnd(input_rep,convergence_jitter);

        for bin = 2:bin_number
        
            if V(n,bin-1) < V_th
            
                V_derivative = (V_rest - V(n,bin-1))/tau_m;
                
                if bin >= level

                    V(n,bin) = V(n,bin-1) + bin_length*V_derivative + input_rep_var(n,bin) + excitatory_background_noise(n,bin) - inhibitory_background_noise(n,bin);   % membrane potential (after onset of volley)    
            
                else
                    
                    V(n,bin) = V(n,bin-1) + bin_length*V_derivative + excitatory_background_noise(n,bin) - inhibitory_background_noise(n,bin);                          % membrane potential (before onset of volley)
                end
                
             else
            
                 V(n,bin-1) = .03;   % EPSP
                 V(n,bin) = V_rest;
                
                  if occurence_spike == 0
                    
                      latency = [latency bin];
                    
                  end
            
                  occurence_spike = 1;
            
             end
        
        end
   
        figure(1), subplot(P,1,p), plot(time_axis,V), axis([0 epoch -.1 .1]), ylabel('V_m (V)')    % membrane potential over time (internally driven pools)
        if p == P
        xlabel('time (s)')
        end 
    
    end
    
    mean_latency = [mean_latency mean(latency)];
    standard_deviation = [standard_deviation std(latency)];
    latency =[];
    
    level = level + conduction_bins;   % update pool level
     
    V_1 = V(1,:);
    V_2 = V(2,:);
    for bin = 1:bin_number
        if V_1(bin) == .03
            V_1(bin) = 1;
        else 
            V_1(bin) = 0;
        end
        if V_2(bin) == .03
            V_2(bin) = 1;
        else 
            V_2(bin) = 0;
        end
    end
    crossco = xcov(V_1,V_2);
    figure(2)
    subplot(P,1,p), plot(crossco_time,crossco), axis([-epoch epoch -.5 Inf])    % cross correlograms of two first neurons of internally driven pools
    if p == P
    xlabel('time (s)')
    end

end

mean_latency = mean_latency*bin_length                  % print vector of mean latencies of first spikes of pools
standard_deviation = standard_deviation*bin_length      % print vector of standdard deviations of first spikes of pools

figure(3)
plot(mean_latency,standard_deviation,'bs','MarkerFaceColor','b','MarkerSize',6)  % standard deviation vs mean latencies of first spikes of pools
xlabel('mean latency (s)'), ylabel('standard deviation')
lsline
cc = corrcoef(mean_latency,standard_deviation);
correlation_coefficient = cc(1,2)