% Derrick Roos
%
% LF_network.m - This function rules the Lansner-Fransen type simulations as
% detailed in their paper - "Modelling Hebbian cell assemblies comprised of
% cortical neurons". There are an equal number of excitory and inhibitory
% cells. Each excitory cell has a companion inhibitory cell which synapses
% onto it. The network is trained with different patterns loaded in. Then
% its calls a code which produces connection strengths. Positive strengths
% form Excitory to Excitory connections. Negative strengths form
% connections from Excitory to Inhibitory cells.
%
% v = LF_network(dt,Tcutoff,Tfin,I0,stim_cells,Patterns)
%      where: dt = time step
%        Tcutoff = time cutoff for stimulated cells
%           Tfin = time end simulation
%             I0 = stimulation size for cells
%             stim_cells = vector of cells recieving stimulus
%            Patterns = matrix where each row is a pattern the network is
%                       trained with
%
%      returns: v = Voltages for excitory cells, indexed as
%                           v(cell,timestep,compartment)
%
%
% ex.    p = pattern_maker(40,8,8)
%        v = LF_network(.01,25,45,1,[1 6 8 10 25],p);

function [v,Fsv] = LF_network(dt,Tcutoff,Tfin,I0,stim_cells,Patterns)

% find out number of patterns and number of cells from Pattern matrix
[n_patterns,n_cells] = size(Patterns);

% Patterns = pattern_maker(n_cells,2,2);

% Use Tony Kellems Code to produce Wieght matrix, each pattern is given equal signifigance
W = make_W(Patterns,ones(n_patterns,1));

% Use D.Roos code to produce G connection/conduction matrix ( splits into new network with
% excitory cells and inhibitory) hence all entries in G are positive. G(i,j) represents
% strength from  cell j to cell i. First half of rows are connections to
% excitory cells. Second half of rows are inhibitory cells. Constants
% entered represent the ideal average value. See make_G for more detail
G = make_G(W,0.01,0.007,0.03);

% time that excitory synapse stays on. Time and strength where tested on
% p_cell_syn to achieve goal in paper of about 3.6 mV EPSP and 20 ms half
% width
exc_syn_time = 3;
% time that inhibitory synapse stays on. Time and strength where tested on
% p_cell_syn to achieve goal in paper of about -5 mV IPSP and 10 ms half
% width
in_syn_time = 7;

% FS-cell parameters, FS - Fast spiking cells, the model used for all
% inhibitory cells.
Fs.comps = 2;
Fs.Gm = 0.0016;
Fs.Cm = 0.016;
Fs.Gm(2:Fs.comps) = 0.0096;
Fs.Cm(2:Fs.comps) = 0.288;
Fs.Gcore = 0.0638;
Fs.GNa = 1; % mewS
Fs.vNa = 50;
Fs.vK = -90;
Fs.GK = 1;
Fs.vCa = 150;
Fs.GCa = 0;
Fs.GKCA = 0.01;
Fs.Gleak = Fs.Gm;
Fs.vLeak = -70;
Fs.APinflux = 0.013;
Fs.APdecay = 20;
Fs.v1 = -70;

% p-cell parameters, Model used for all excititory cells
comps = 4;
Gm = 0.0032; %soma
Cm = 0.032; % soma
Gm(2:comps) = 0.0096;
Cm(2:comps) = 0.288;
Gcore = 0.04;
GNa = 1; % mewS
vNa = 40;
vK = -70;
GK = .5;
vCa = 150;
GCa = 0;
GKCA = 0.0017;
Gleak = Gm;
vLeak = -50;
APinflux = 4;
APdecay = 30;
v1 = -50;

% number of timesteps
N = ceil(Tfin/dt);

t = zeros(N,1);   % allocate space for long vectors

% FS - preallocate ion gating variables, voltage and more
Fs.v = zeros(n_cells,N,Fs.comps);
Fs.n = zeros(n_cells,N);
Fs.m = Fs.n;
Fs.h = Fs.n;
Fs.q = Fs.n;
Fs.CaAP = Fs.n;
Fs.ISyn = zeros(n_cells,Fs.comps);

% p-cells - preallocate ion gating variables, voltage and more
v = zeros(n_cells,N,comps);
n = zeros(n_cells,N);
m = n;
h = n;
q = n;
CaAP = n;
Istim = zeros(n_cells,1);
ISyn = zeros(n_cells,comps);

% find intial points for gating variables
Fs.n(:,1) = Fan(Fs.v1)/(Fan(Fs.v1)+Fbn(Fs.v1));
Fs.m(:,1) = Fam(Fs.v1)/(Fam(Fs.v1)+Fbm(Fs.v1));
Fs.h(:,1) = Fah(Fs.v1)/(Fah(Fs.v1)+Fbh(Fs.v1));
Fs.q(:,1) = Faq(Fs.v1)/(Faq(Fs.v1)+Fbq(Fs.v1));
Fs.CaAP(1:n_cells,1)=0;
Fs.v(1:n_cells,1,:) = Fs.v1;

n(1:n_cells,1) = an(v1)/(an(v1)+bn(v1));
m(1:n_cells,1) = am(v1)/(am(v1)+bm(v1));
h(1:n_cells,1) = ah(v1)/(ah(v1)+bh(v1));
q(1:n_cells,1) = aq(v1)/(aq(v1)+bq(v1));
CaAP(1:n_cells,1)=0;
v(1:n_cells,1,:) = v1;

% vector containing time remaining for cells synapse to remain on
S = zeros(2*n_cells,1);


%%% ENTER MAIN LOOP THAT TIMESTEPS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for j=2:N,
    t(j) = (j-1)*dt;

    %%%%%%%% DO FAST SPIKING/ INHIBITORY CELL UPDATE %%%%%%%%%%%%%

    % Update Soma For Fast Spiking Cells, use hybrid euler scheme
    Fs.n(:,j) = ( Fs.n(:,j-1) + dt*Fan(Fs.v(:,j-1,1)) )./(1 + dt*(Fan(Fs.v(:,j-1,1))+Fbn(Fs.v(:,j-1,1))) );
    Fs.m(:,j) = ( Fs.m(:,j-1) + dt*Fam(Fs.v(:,j-1,1)) )./(1 + dt*(Fam(Fs.v(:,j-1,1))+Fbm(Fs.v(:,j-1,1))) );
    Fs.h(:,j) = ( Fs.h(:,j-1) + dt*Fah(Fs.v(:,j-1,1)) )./(1 + dt*(Fah(Fs.v(:,j-1,1))+Fbh(Fs.v(:,j-1,1))) );
    Fs.q(:,j) = ( Fs.q(:,j-1) + dt*Faq(Fs.v(:,j-1,1)) )./(1 + dt*(Faq(Fs.v(:,j-1,1))+Fbq(Fs.v(:,j-1,1))) );
    Fs.CaAP(:,j) = ( dt*(Fs.vCa - Fs.v(:,j-1,1)).*Fs.APinflux.*Fs.q(:,j).^5 + Fs.CaAP(:,j-1,1) )./( 1+dt*Fs.APdecay );

    Fs.top = Fs.Cm(1)*Fs.v(:,j-1,1)+dt*(Fs.Gleak(1)*Fs.vLeak + Fs.GNa*Fs.m(:,j).^3.*Fs.h(:,j)*Fs.vNa + Fs.GK*Fs.n(:,j).^4*Fs.vK + Fs.GCa*Fs.q(:,j).^5*Fs.vCa +Fs.GKCA*Fs.vK*Fs.CaAP(:,j) + Fs.ISyn(:,1));

    Fs.bot = Fs.Cm(1) + dt*(Fs.Gleak(1) + Fs.GNa*Fs.m(:,j).^3.*Fs.h(:,j) + Fs.GK*Fs.n(:,j).^4 + Fs.GCa*Fs.q(:,j).^5 + Fs.GKCA*Fs.CaAP(:,j));

    % Add compartmental current
    Fs.top = Fs.top + dt*Fs.Gcore*Fs.v(:,j-1,2);
    Fs.bot = Fs.bot + dt*Fs.Gcore;
    Fs.v(:,j,1) = Fs.top./Fs.bot;

    for k = 2:Fs.comps

        Fs.top = Fs.Cm(k)*Fs.v(:,j-1,k)+dt*(Fs.Gleak(k)*Fs.vLeak + Fs.ISyn(:,k));
        Fs.bot = Fs.Cm(k) + dt*Fs.Gleak(k);

        % check for right neighbor
        if k < Fs.comps
            Fs.top = Fs.top + dt*Fs.Gcore*Fs.v(:,j-1,k+1);
            Fs.bot = Fs.bot + dt*Fs.Gcore;
        end
        % left neighbor compartment
        Fs.top = Fs.top + dt*Fs.Gcore*Fs.v(:,j-1,k-1);
        Fs.bot = Fs.bot + dt*Fs.Gcore;


        Fs.v(:,j,k) = Fs.top./Fs.bot;

    end

    %%%%%%%% DO P-CELL/ EXCITITORY CELL UPDATE %%%%%%%%%%%%%

    % update injected stimulus
    Istim(stim_cells) = I0*(t(j)>2)*(t(j)<Tcutoff);

    % Update Soma
    n(:,j) = ( n(:,j-1) + dt*an(v(:,j-1,1)) )./(1 + dt*(an(v(:,j-1,1))+bn(v(:,j-1,1))) );
    m(:,j) = ( m(:,j-1) + dt*am(v(:,j-1,1)) )./(1 + dt*(am(v(:,j-1,1))+bm(v(:,j-1,1))) );
    h(:,j) = ( h(:,j-1) + dt*ah(v(:,j-1,1)) )./(1 + dt*(ah(v(:,j-1,1))+bh(v(:,j-1,1))) );
    q(:,j) = ( q(:,j-1) + dt*aq(v(:,j-1,1)) )./(1 + dt*(aq(v(:,j-1,1))+bq(v(:,j-1,1))) );
    CaAP(:,j) = ( dt*(vCa - v(:,j-1,1)).*APinflux.*q(:,j).^5 + CaAP(:,j-1,1) )./( 1+dt*APdecay );

    top = Cm(1)*v(:,j-1,1)+dt*(Gleak(1)*vLeak + GNa*m(:,j).^3.*h(:,j)*vNa + GK*n(:,j).^4*vK + GCa*q(:,j).^5*vCa +GKCA*vK*CaAP(:,j)+ + ISyn(:,1) + Istim(:));

    bot = Cm(1) + dt*(Gleak(1) + GNa*m(:,j).^3.*h(:,j) + GK*n(:,j).^4 + GCa*q(:,j).^5 + GKCA*CaAP(:,j));

    % Add compartmental current
    top = top + dt*Gcore*v(:,j-1,2);
    bot = bot + dt*Gcore;
    v(:,j,1) = top./bot;

    for k = 2:comps

        top = Cm(k)*v(:,j-1,k)+dt*(Gleak(k)*vLeak + ISyn(:,k));
        bot = Cm(k) + dt*Gleak(k);

        % check for right neighbor
        if k < comps
            top = top + dt*Gcore*v(:,j-1,k+1);
            bot = bot + dt*Gcore;
        end
        % left neighbor compartment
        top = top + dt*Gcore*v(:,j-1,k-1);
        bot = bot + dt*Gcore;


        v(:,j,k) = top./bot;

    end

    %%%%% UPDATE SYNAPTIC CURRENTS %%%%%%%%%%%%%%%%%%%

    % shorten time remaining for active synapses by timestep
    S = max(0,S - dt);

    % find p_cells that have spiked (checks when crosses 0 voltage line so
    % techniquelly will reset when voltage passes 0 mV  on way up and way
    % back down).  Then reset time these excitory synapse will be active
    % for
    S(find(v(:,j,1).*v(:,j-1,1)< 0)) = exc_syn_time; % reset synapse release timer for fired cells

    % find fs cells that have spiked and do the same procedure as above for
    % the fs cells.
    S(find(Fs.v(:,j,1).*Fs.v(:,j-1,1)< 0) + n_cells) = in_syn_time;

    % normalize such that get vector with 1's for active cell synapses, 0
    % for inactive cell synapses.
    act_cells = S./max(S,0.01);


    % excitory cells synapse onto compartment 4. Update by using synapse
    % model of  conductance * active * (Syn_rev - voltage).  Here the
    % reversal potential for the excitory synapse is 0 mV.
    ISyn(:,3) = G(1:n_cells,1:n_cells)*act_cells(1:n_cells).*(-v(:,j,3));


    % inhibitory cells synapse onto compartment 1. Here same synapse model
    % is being used. Reversal potential for inhibitory synapse here is set
    % to -90mV.
    ISyn(:,1) = G(1:n_cells,n_cells+1:end)*act_cells(n_cells+1:end).*(-90-Fs.v(:,j,1));

    %find synaptic current going into each inhibitory cell. Again 0mV
    %synapse reversal potential is being used.
    Fs.ISyn(:,2) = G(n_cells+1:end,1:n_cells)*act_cells(1:n_cells).*(-v(:,j,2));


end

Fsv = Fs.v;

% FS-Cells Gating function needed
function val = Fam(v)
A = 0.2; B = -30; C = 1;
val = A*(v-B)./(1-exp((B-v)/C));

function val = Fbm(v)
A= 0.06 ; B = -39 ; C = 20;
val = A*(B-v)./(1-exp((v-B)/C));

function val = Fah(v)
A= 0.08 ; B = -30; C = 0.2;
val = A*(B-v)./(1-exp((v-B)/C));

function val = Fbh(v)
A= 0.4 ; B = -26; C = 0.2;
val = A./(1+exp((B-v)/C));

function val = Fan(v)
A= 0.02 ; B = -21; C = 0.2;
val = A*(v-B)./(1-exp((B-v)/C));

function val = Fbn(v)
A= 20 ; B = -18; C = 0.2;
A= .02 ; B = -18; C = 0.2;
val = A*(B-v)./(1-exp((v-B)/C));

function val = Faq(v)
% A=  0.08; B = -10; C = 11;
A = 0.08; B = -15; C = 1;
val = A*(v-B)./(1-exp((B-v)/C));

function val = Fbq(v)
% A= 0.001 ; B = -10; C = 0.5;
A = 0.005; B = -10; C = 20;

val = A*(B-v)./(1-exp((v-B)/C));




% functions for P-Cell function
function val = am(v)
A = 0.2; B = -40; C = 1;
val = A*(v-B)./(1-exp((B-v)/C));

function val = bm(v)
A= 0.06 ; B = -49 ; C = 20;
val = A*(B-v)./(1-exp((v-B)/C));

function val = ah(v)
A= 0.08 ; B = -40; C = 1;
val = A*(B-v)./(1-exp((v-B)/C));

function val = bh(v)
A= 0.4 ; B = -36; C = 2;
val = A./(1+exp((B-v)/C));

function val = an(v)
A= 0.02 ; B = -15; C = 0.8;
val = A*(v-B)./(1-exp((B-v)/C));

function val = bn(v)
A= 0.04 ; B = -40; C = 0.4;
val = A*(B-v)./(1-exp((v-B)/C));

function val = aq(v)
% A=  0.08; B = -10; C = 11;
A = 0.08; B = -25; C = 1;
val = A*(v-B)./(1-exp((B-v)/C));

function val = bq(v)
% A= 0.001 ; B = -10; C = 0.5;
A = 0.005; B = -20; C = 20;
val = A*(B-v)./(1-exp((v-B)/C));