function  [v, spk]  = liandfr4(t_vect, g_syn_ex, v_rev_ex, g_syn_in, v_rev_in)
%
% [v, spk] = liandfr4(t_vect, g_syn_ex, v_rev_ex, g_syn_in, v_rev_in)
%
%liandfr4: returns the membrane potential and spike train vector 
%for a leaky integrate and fire neuron receiving input 
%in terms of conductance changes (units: nS) at times t_vect. The 
%scalar v_rev is the reversal potential of the synapse. 

%convert the synaptic conductance from nS to microS for consistency with 
%other variables below
g_syn_ex = 1e-3*g_syn_ex;
g_syn_in = 1e-3*g_syn_in;

%variable initialization
t_step = t_vect(2)-t_vect(1);
n_t = length(t_vect);
v = zeros(1,n_t);
spk = zeros(1,n_t);

%LIF parameters
vr = -70; %in mV 
vth = -54; %in mV
tref = 2; %in msec
tau = 30; %in msec
rin = 20; %in MOhms

vth0 = vth - vr; %threshold with respect to rest
tref0 = round(tref/t_step); %refractory period in units of the time step

%save time computationally compute in advance the constants needed
%for forward Euler
f1 = 1 - t_step/tau;
f2 = ((rin*t_step)/tau)*g_syn_ex;
f3 = ((rin*t_step)/tau)*g_syn_in;

%backward Euler constant
b0 = tau/t_step;
b1ex = 1./(1 + (tau/t_step) + (rin*g_syn_ex));
b1in = 1./(1 + (tau/t_step) + (rin*g_syn_in));
b1 = b1ex + b1in;
b2ex = rin*v_rev_ex*g_syn_ex;
b2in = rin*v_rev_in*g_syn_in;
b2 = b2ex + b2in;

v(1) = 0; %initial condition
ref_time = 0; %no refractory period to start with

for i = 2:n_t
    if ( ref_time <= 0 )
        %refractory period over, integrate membrane
        %potential with forward or backward Euler
        
        %forward Euler
        v(i) = f1*v(i-1) - f2(i-1)*(v(i-1)-v_rev_ex) - f3(i-1)*(v(i-1)-v_rev_in); 
        %v(i) = b1*(b0*v(i-1) + b2(i)); %backward Euler
        
        if (v(i) > vth0)
            %spike reset membrane potential 
            spk(i) = 1;
            v(i) = 0;
            ref_time = tref0;
        end;
    else
        %decrement refractory time; wait until refractory period
        %is finished 
        ref_time = ref_time - 1;
    end;
end;