% Jim Wang
% 9/9/08
%
% Katie Ward
% 9/23/08
%
% Mingbo Cai
% 10/21/08
% slightly modify Katie's program to draw activity profile
%
% 
%
% dt = time step (ms)
% T = total time of simulation (ms)
% theta0 = input PO
% eps = degree of tuning of LGN input, 0 <= eps <= 1
% lambda_e = space constant of excitatory neuron, in degree
% lambda_i = space constant of inhibitory neuron, in degree
% fbar = total firing rate of LGN input
% Ginp = maximal conductance of input synapses (mS/cm^2)
% NGbar_e = maximal conductance of excitatory neuron multiplied by neuron number(mS/cm^2)
% NGbar_i = maximal conductance of inhibitory neuron multiplied by neuron number(mS/cm^2)
% tau = decay constant for the conductance of all the neuron (ms)
% Ne_net = number of excitatory neurons
% Ni_net = number of inhibitory neurons
%
%
% Example: in order to generate Figure 4, use
% spkfeq = heas3(0.1,500,0,0,11.5,43,2700,0.0025,0.143,0.333,3,160,160), but
% it cannot repeat the result of Fig 4.(Running with 1600 neurons is
% usually too slow on PC, so here we use 160)
% For Figure 5, use the parameter
% spkfeq = heas3(0.1,500,0,0,11.5,43,2700,0.0025,0.2,1.33,3,160,160), but still
% not successful


function spkfreq = heas3(dt,T,theta0,eps,lambda_e,lambda_i,fbar,Ginp,NGbar_e,NGbar_i,tau,Ne_net,Ni_net)

close all

theta0 = theta0 * pi / 180;  %orientation of stimulus (radians)

v0 = -73;  %mV

Cm = 1;         % �F/cm^2

gbar_na = 100;  % mS/cm^2
gbar_k = 40;    % mS/cm^2
g_L = .05;

E_na = 55;      % mV
E_k = -80;      % mV
EL = -65;      % mV


phi = 0.1;

tau_A = 20;  %ms
gbar_A = 20;  %mS/cm^2


Ein = -80;  %mV
Eex = 0;  %mV
tau_inp = tau;  %ms   %%%VALUE?
Einp = 0;   %%%SAME AS Eex?


% initializing
%Ne_net = 100;
%Ni_net = 80;   %%%Don't know why, but don't use Ni_net = 1.  Other values seem to work.

Gbar_e = NGbar_e / Ne_net;
Gbar_i = NGbar_i / Ni_net;

ee = zeros(Ne_net,1);
Mee = zeros(Ne_net,Ne_net);
Mei = zeros(Ne_net,Ni_net);

%G = .01;  %mS/cm^2
profile_E = zeros(ceil(T/dt),Ne_net);
profile_I = zeros(ceil(T/dt),Ni_net);
net_E = struct('v', ee+v0, ...
    'Esyn', Eex, ...
    'gI', Mei, ...
    'gE', Mee, ...
    'gElgn', ee, ...
    'tau_gE', Mee + tau, ...
    'tau_gI', Mei + tau, ...
    'spks', [], ...
    'N', Ne_net, ...
    'h', ee, ...
    'n', ee, ...
    'b', ee, ...
    'm', ee, ...
    'finp', ee, ...
    'theta', ee, ...
    'lambda', lambda_e*pi/180, ...
    'Gbar', Gbar_e);

net_E.theta = (-pi/2) + (pi/Ne_net)*([1:net_E.N]');  %PO's of network cells

%eps = .5;
C = .5;
%fbar = 1570;  %max firing rate (spikes/sec)
net_E.finp = fbar*C * ((1-eps)+eps*cos(2*(net_E.theta-theta0)));

trains = generate_spike_trains_vec(Ne_net,net_E.finp,T/1000);
trains(:,1) = trains(:,1) * 1000;  %convert to ms

ii = zeros(Ni_net,1);

net_I = struct('v', ii+v0, ...
    'Esyn', Ein, ...
    'gE', Mei', ...
    'tau_gE', Mei'+ tau, ...
    'spks', [], ...
    'N', Ni_net, ...
    'h', ii, ...
    'n', ii, ...
    'b', ii, ...
    'm', ii, ...
    'theta', ii, ...
    'lambda', lambda_i*pi/180, ...
    'Gbar', Gbar_i);

net_I.theta = (-pi/2) + (pi/Ni_net)*([1:net_I.N]');  %PO's of network cells

%Construct weight matrices
%Ginp = .01;  %mS/(cm^2)





GEtoE = Mee;
GItoE = Mei;
for i = 1:Ne_net
    dtheta = abs(net_E.theta(i) - net_E.theta);
    GEtoE(i,:) = ((net_E.Gbar/net_E.lambda) * exp(-dtheta/net_E.lambda))';

    if Ni_net ~= 0
        dtheta = abs(net_E.theta(i) - net_I.theta);
        GItoE(i,:) = ((net_I.Gbar/net_I.lambda) * exp(-dtheta/net_I.lambda))';
    end
end
GEtoE = GEtoE - spdiags(diag(GEtoE),0,Ne_net,Ne_net);   %no cell synapses onto itself

GEtoI = Mei';
for i = 1:Ni_net
    dtheta = abs(net_I.theta(i) - net_E.theta);
    GEtoI(i,:) = ((net_E.Gbar/net_E.lambda) * exp(-dtheta/net_E.lambda))';

end

% Visualize weight matrices
figure(1)
colormap(1-hot)
subplot(3,1,1)
imagesc(GEtoE)
title('GEtoE')
colorbar

xlabel('pres')
ylabel('posts')
subplot(3,1,2)
imagesc(GEtoI)
colorbar
title('GEtoI')

subplot(3,1,3)
imagesc(GItoE)
title('GItoE')
colorbar
drawnow




%Visualize Poisson spike trains
figure(2)
plot(trains(:,1), trains(:,2),'.')
xlabel('time (ms)')
ylabel('exc. neuron')
title('Spike trains feeding excitatory network neurons')

vthreshold = -20;    % threshold for firing (mV)
net_E.v(:) = v0;          % resting membrane potential (mV)
net_I.v(:) = v0;
net_E.gE = net_E.gE*0;
net_E.gElgn = net_E.gElgn * 0;
net_E.spks = [];
%net_E.Gbar = Gbar_e;

%Initiate values for external spikes
ct = 1;
spk1 = ee;

%Variables needed for spike counting
aboveE = ee;
belowE = ee;
aboveI = ii;
belowI = ii;
spkcnt = ee;
spkcnt1 = ii;
%timeseries = zeros(Ne_net,ceil(T/dt));
j = 1;
while j*dt < T

    %Find network neurons receiving lgn input
    ind = 0;
    spk1 = 0*spk1;
    while trains(ct,1) < j*dt
        ind = ind + 1;
        spikedNeur = trains(ct,2);
        spk1(ind) = spikedNeur;
        ct = ct + 1;
    end
    lgn_spks = spk1(1:ind);  %indeces of net neurons receiving lgn input




    %Values needed for I_active
    alpha_mE = -0.1*(net_E.v+30) ./ (exp(-0.1*(net_E.v+30))-1);
    alpha_mI = -0.1*(net_I.v+30) ./ (exp(-0.1*(net_I.v+30))-1);

    beta_mE = 4*exp(-(net_E.v+55) / 18);
    beta_mI = 4*exp(-(net_I.v+55) / 18);

    alpha_hE = 0.07*exp(-(net_E.v+44) / 20);
    alpha_hI = 0.07*exp(-(net_I.v+44) / 20);

    beta_hE = 1 ./ (exp(-0.1*(net_E.v+14))+1);
    beta_hI = 1 ./ (exp(-0.1*(net_I.v+14))+1);

    alpha_nE = -0.01*(net_E.v+34) ./ (exp(-0.1*(net_E.v+34))-1);
    alpha_nI = -0.01*(net_I.v+34) ./ (exp(-0.1*(net_I.v+34))-1);

    beta_nE = 0.125*exp(-(net_E.v+44) / 80);
    beta_nI = 0.125*exp(-(net_I.v+44) / 80);

    a_infE = 1 ./ (exp(-(net_E.v+50)/20)+1);
    a_infI = 1 ./ (exp(-(net_I.v+50)/20)+1);

    b_infE = 1 ./ (exp((net_E.v+80)/6)+1);
    b_infI = 1 ./ (exp((net_I.v+80)/6)+1);

    m_infE = alpha_mE ./ (alpha_mE+beta_mE);
    m_infI = alpha_mI ./ (alpha_mI+beta_mI);

    tau_hE = phi ./ (alpha_hE+beta_hE);
    tau_hI = phi ./ (alpha_hI+beta_hI);

    h_infE = alpha_hE ./ (alpha_hE+beta_hE);
    h_infI = alpha_hI ./ (alpha_hI+beta_hI);

    tau_nE = phi ./ (alpha_nE+beta_nE);
    tau_nI = phi ./ (alpha_nI+beta_nI);

    n_infE = alpha_nE ./ (alpha_nE+beta_nE);
    n_infI = alpha_nI ./ (alpha_nI+beta_nI);

    if i == 1  %initial conditions
        net_E.h = h_infE;
        net_I.h = h_infI;

        net_E.n = n_infE;
        net_I.n = n_infI;

        net_E.b = b_infE;
        net_I.b = b_infI;
    end

    net_E.m = m_infE;
    net_I.m = m_infI;

    net_E.h = (net_E.h .* tau_hE + dt*h_infE) ./ (tau_hE+dt);  % hybrid euler
    net_I.h = (net_I.h .* tau_hI + dt*h_infI) ./ (tau_hI+dt);

    net_E.n = (net_E.n .* tau_nE + dt*n_infE) ./ (tau_nE+dt);
    net_I.n = (net_I.n .* tau_nI + dt*n_infI) ./ (tau_nI+dt);

    net_E.b = (net_E.b*tau_A + dt*b_infE) / (tau_A+dt);
    net_I.b = (net_I.b*tau_A + dt*b_infI) / (tau_A+dt);




    %Conductances decay
    net_E.gE = net_E.gE ./ (1+(dt ./ net_E.tau_gE));
    net_E.gI = net_E.gI ./ (1+(dt ./ net_E.tau_gI));
    net_I.gE = net_I.gE ./ (1+(dt ./ net_I.tau_gE));

    net_E.gElgn = net_E.gElgn / (1+(dt / tau_inp));

    %Add synaptic input to conductances
    net_E.gE(:,net_E.spks) = net_E.gE(:,net_E.spks) + GEtoE(:,net_E.spks);
    net_E.gI(:,net_I.spks) = net_E.gI(:,net_I.spks) + GItoE(:,net_I.spks);
    net_I.gE(:,net_E.spks) = net_I.gE(:,net_E.spks) + GEtoI(:,net_E.spks);

    net_E.gElgn(lgn_spks) = net_E.gElgn(lgn_spks) + Ginp;



    % Update net_E.v
    tm1 = g_L*EL;  %I_leak tm

    tm2 = gbar_na * (m_infE .^ 3) .* net_E.h * E_na;  %I_active tms
    tm3 = gbar_k * (net_E.n .^ 4) * E_k;
    tm4 = gbar_A * (a_infE .^ 3) .*net_E.b * E_k;

    tm5 = (sum(net_E.gE,2) .* (Eex - net_E.v)) + (sum(net_E.gI,2) .* (Ein - net_E.v));  %I_net tm

    tm6 = net_E.gElgn .* (Einp - net_E.v);  %I_ext tm

    top = net_E.v + (dt/Cm) * (tm1 + tm2 + tm3 + tm4 + tm5 + tm6);
    
    
    tm1 = g_L;  %I_leak tm

    tm2 = gbar_na * (m_infE .^ 3) .* net_E.h;  %I_active tms
    tm3 = gbar_k * (net_E.n .^ 4);
    tm4 = gbar_A * (a_infE .^ 3) .* net_E.b;

    tm5 = sum(net_E.gE,2) + sum(net_E.gI,2);  %I_net tm

    tm6 = net_E.gElgn;  %I_ext tm

    bottom = 1 + (dt/Cm) * (tm1 + tm2 + tm3 + tm4 + tm5 + tm6);
    net_E.v = top ./ bottom;


    %Update net_I.v
    tm1 = g_L*EL;  %I_leak term

    tm2 = gbar_na * (m_infI .^ 3) .* net_I.h * E_na;  %I_active terms
    tm3 = gbar_k * (net_I.n .^ 4) * E_k;
    tm4 = gbar_A * (a_infI .^ 3) .*net_I.b * E_k;

    tm5 = (sum(net_I.gE,2) .* (Eex - net_I.v));  %I_net tm

    top = net_I.v + (dt/Cm) * (tm1 + tm2 + tm3 + tm4 + tm5);
    

    tm1 = g_L;  %I_ext term

    tm2 = gbar_na * (m_infI .^ 3) .* net_I.h;  %I_active terms
    tm3 = gbar_k * (net_I.n .^ 4);
    tm4 = gbar_A * (a_infI .^ 3) .* net_I.b;

    tm5 = sum(net_I.gE,2);  %I_net term

    bottom = 1 + (dt/Cm) * (tm1 + tm2 + tm3 + tm4 + tm5);

    net_I.v = top ./ bottom;

    %Check for spikes
    belowE = 1*(net_E.v < vthreshold);
    prod = 1*belowE.*aboveE;  %spk if i was above and is now below
    spkcnt = spkcnt + prod;
    net_E.spks = find(prod);
    aboveE = 1*(net_E.v >= vthreshold);
    profile_E(j,net_E.spks) = 1;

    belowI = 1*(net_I.v < vthreshold);
    prod = belowI.*aboveI;  %spk if i was above and is now below
    spkcnt1 = spkcnt1 + prod;
    net_I.spks = find(prod);
    aboveI = 1*(net_I.v >= vthreshold);
    profile_I(j,net_I.spks) = 1;
    %timeseries(:,j) = net_E.v;
    j = j + 1;
 
end
j = j - 1;

spkfreq = spkcnt/T * 1000;  %Hz
spkfreq1 = spkcnt1/T * 1000; 

figure(3)
subplot(2,1,1)
plot([1:Ne_net]',spkfreq,'.')
xlabel('neuron')
ylabel('spike freq (Hz)')
title('Activity profile of excitatory neurons')
subplot(2,1,2)
plot([1:Ni_net]',spkfreq1,'.')
xlabel('neuron')
ylabel('spike freq (Hz)')
title('Activity profile of inhibitory neurons')

figure(4)
colormap gray;
subplot(2,1,1)
imagesc(profile_E');
ylabel('neuron')
xlabel('time')
title('temporal-spatial pattern of excitatory network');

subplot(2,1,2)
imagesc(profile_I');
ylabel('neuron')
xlabel('time')
title('temporal-spatial pattern of inhibitory network');