%% Full model of LIF with inhibitory and electrical connections
%  Individual applied current for each neuron (I_1 and I_2)
%  What happens when you adjust I_2 to 1.2, 1.25, 1.3?
%  Explore effects of changing parameters: alpha, g_s, beta, g_c, I_1 and I_2
clear;
close all;

%% Change these parameters!

% Applied current for each cell
I_1=1.1;
I_2=1.1;

% parameters for the inhibitory connection
alpha=3;
g_s=0.2;

% parameters for the electrical connection; default is no coupling
beta=0.2;
g_c=0;

%% Numerical code begins for solving differential equation
time_step=0.0001;
total_time=50;

v_1 = -1*ones(round(total_time/time_step),1);
v_2 = -1*ones(round(total_time/time_step),1);

t_sim=[1:length(v_1)]*time_step;

t_spike1=0;
t_spike2=0;

% initial conditions
v_1(1)=0.59;
v_2(1)=0;

I_syn2=zeros(round(total_time/time_step),1);
I_syn1=zeros(round(total_time/time_step),1);

spike1_count=0;
spike2_count=0;
for t = 2:round(total_time/time_step)
    
    % spike if cell 1 is above threshold
    if v_1(t-1)>1 && v_1(t-1)~=3
        v_1(t)=3;
        
        spike1_count=spike1_count+1;
        t_spike1(spike1_count)=t_sim(t-1); 
    end
    
    % spike if cell 2 is above threshold
    if v_2(t-1)>1 && v_2(t-1)~=3
        v_2(t) = 3;
       
        spike2_count=spike2_count+1;
        t_spike2(spike2_count)=t_sim(t-1);  
    end
    
    % total synaptic input into cell 2 from cell 1
    if spike1_count~=0
        for i=1:length(t_spike1)
            I_syn2(t)=I_syn2(t)+alpha^2*(t_sim(t)-t_spike1(i))*exp(-alpha*(t_sim(t)-t_spike1(i)));
        end
    end
    
    % total synaptic input into cell 1 from cell 2
    if spike2_count~=0
        for i=1:length(t_spike2)
            I_syn1(t)=I_syn1(t)+alpha^2*(t_sim(t)-t_spike2(i))*exp(-alpha*(t_sim(t)-t_spike2(i)));
        end
    end
    
    % Perform Euler integration for cell 2
    if v_2(t) ~= 3
        if v_2(t-1)==3 && v_1(t-1)~=3
            v_2(t) = 0 + time_step*(-v_2(t-1)+I_2-g_s*I_syn2(t-1)...
                + g_c*(v_1(t-1)-v_2(t-1)));
        elseif v_2(t-1)~=3 && v_1(t-1)~=3
            v_2(t) = v_2(t-1) + time_step*(-v_2(t-1)+I_2-g_s*I_syn2(t-1)...
                + g_c*(v_1(t-1)-v_2(t-1)));
        else
            v_2(t) = v_2(t-1) + time_step*(-v_2(t-1)+I_2-g_s*I_syn2(t-1)...
                + g_c*(0-v_2(t-1)));
        end
        
        
        if v_1(t-1)==3
            v_2(t) = v_2(t)+g_c*beta;
        end
    end
    
    % Perform Euler integration for cell 2
    if v_1(t) ~=3
        
        if v_1(t-1)==3 && v_2(t-1)~=3
            v_1(t) = 0 + time_step*(-v_1(t-1)+I_1-g_s*I_syn1(t-1));
        elseif v_1(t-1)~=3 && v_2(t-1)~=3
            v_1(t) = v_1(t-1) + time_step*(-v_1(t-1)+I_1-g_s*I_syn1(t-1)...
                + g_c*(v_2(t-1)-v_1(t-1)));
        else
            v_1(t) = v_1(t-1) + time_step*(-v_1(t-1)+I_1-g_s*I_syn1(t-1)...
                + g_c*(0-v_1(t-1)));
        end
        
        
        if v_2(t-1)==3
            v_1(t) = v_1(t)+g_c*beta;
        end
    end

    
end


%% Plotting begins

subplot(2,1,1)
plot(t_sim,v_1,'red','linewidth',1.5)
hold on;
plot(t_sim,v_2,'linewidth',1.5);
axis([10 24 -.1 3.1])
 
legend('Cell 1', 'Cell 2')

ylabel('v','fontsize',16)
xlabel('Time','fontsize',16)
set(gca,'fontsize',16)

subplot(2,1,2)
plot(t_sim,-I_syn1*g_s,'red','linewidth',1.5)
hold on;
plot(t_sim,-I_syn2*g_s,'blue','linewidth',1.5)

legend('Synaptic Input for Cell 1', 'Synaptic Input for Cell 2')

axis([10 24 -0.3 0.1])
    
ylabel('I_{syn}','fontsize',16)
xlabel('Time','fontsize',16)
set(gca,'fontsize',16)

