function II_array=EPSC_excitation_response(v,c,excitation,mag_mult)
tic
% global n_k_hh m_na_hh h_na_hh
global a_k_rm b_k_rm c_k_rm m_na_rm h_na_rm n_htk_rm p_htk_rm w_ltk_rm z_ltk_rm w_ltk_kv7 z_ltk_kv7 w_ltk_kvA z_ltk_kvA r_h_rm
global dt dur gb_syn VV plot_syn time EPSC_shape Isyn %(added Isyn on 3/9/2016
tau=excitation; % the average II time
% For calculating the size of the array containing the synaptic excitation
if EPSC_shape==1
EPSC_length=ceil(tau*9.5/dt);
elseif EPSC_shape==2
EPSC_length=25/dt;
elseif EPSC_shape==3
EPSC_length=52/dt;
end
% Preallocating large arrays
VV=zeros(1,dur/dt);
II=zeros(1,dur/5);
% Initial values to initiate backwards difference equations
V(1:2)=v;
[mi_rm variable]=inf_tau_m_rm(V(1));
[hi_rm variable]=inf_tau_h_rm(V(1));
[ai_rm variable]=inf_tau_a_rm(V(1));
[bi_rm variable]=inf_tau_b_rm(V(1));
[ci_rm variable]=inf_tau_c_rm(V(1));
[ni_htk_rm variable]=inf_tau_n_htk_rm(V(1));
[pi_htk_rm variable]=inf_tau_p_htk_rm(V(1));
[wi_ltk_rm variable]=inf_tau_w_ltk_rm(V(1));
[zi_ltk_rm variable]=inf_tau_z_ltk_rm(V(1));
[ri_h_rm variable]=inf_tau_r_rm(V(1));
[wi_ltk_kv7 variable]=inf_tau_w_ltkcnq_rm(V(1)); % initialize these values if you want to have a complex representation of IKL
[zi_ltk_kv7 variable]=inf_tau_z_ltkcnq_rm(V(1));
[wi_ltk_kvA variable]=inf_tau_w_ltk_rm(V(1));
[zi_ltk_kvA variable]=inf_tau_z_ltkA_rm(V(1));
m_na_rm(1:2)=mi_rm;
h_na_rm(1:2)=hi_rm;
a_k_rm(1:2)=ai_rm;
b_k_rm(1:2)=bi_rm;
c_k_rm(1:2)=ci_rm;
n_htk_rm(1:2)=ni_htk_rm;
p_htk_rm(1:2)=pi_htk_rm;
w_ltk_rm(1:2)=wi_ltk_rm;
z_ltk_rm(1:2)=zi_ltk_rm;
r_h_rm(1:2)=ri_h_rm;
w_ltk_kv7(1:2)=wi_ltk_kv7;
z_ltk_kv7(1:2)=zi_ltk_kv7;
w_ltk_kvA(1:2)=wi_ltk_kvA;
z_ltk_kvA(1:2)=zi_ltk_kvA;
gb_syn=0;
Il(3)=I_l(V(2));
Isyn(3)=I_syn(V(2));
Ina_rm=I_na_rm(V(2));
Ik_rm=I_k_rm(V(2));
Ihtk_rm=I_htk_rm(V(2));
Iltk_rm=I_ltk_rm(V(2));
Ih_rm=I_h_rm(V(2));
Il=I_l(V(2));
Isyn=I_syn(V(2));
V(3)=(-2*dt/c*(Isyn+Ina_rm+Ik_rm+Ihtk_rm+Iltk_rm+Ih_rm+Il)+4*V(2)-V(1))*1/3;
V(1)=V(2);
V(2)=V(3);
% EPSC_II=100; % Need to calculate how long it takes for single
% % compartment to reach steady steat resting potential
% % This is the "time" at wich EPSC spike waveforms
% % Will start to be calculated/created
% II_waveform=zeros(1,EPSC_II);
% % for first few steps, the synaptic conductance is zero
% EPSC_II_count=1; % Start countdown to next EPSC spike
% filter_count=1;
% overflow=0;
% H=[1 0.31135116787317596 1; 1 -1.9669151776947449 0.96745369748177656];
% % Coefficents of filter
%
% spike_count=1; % Needed to start the arry of II times for spikes from
% % the single compartment model
% ii_count=1; % Needed to start the counter for interspike intervals
VV=(V);
% Generating the stimulus array
EPSC_train=generate_EPSC_train(EPSC_length,excitation);
plot_syn=EPSC_train/10*mag_mult/97; %EPSC converted to pA (/10), converted to scale factor (*mag_mult), and converted to conductance (/97 mV holding potential of recording)
for zz=4:dur/dt
% % % % % % % Generates the EPSC spikes (adds a new one once the intervent
% % % % % % % interval "runs out", and then calculates a new intervent interval)
% % % % % % if EPSC_II_count>=EPSC_II
% % % % % % if EPSC_II_count < EPSC_length % Still possibility of residual conductance from previous EPSC... This value needs to be changed if alpha is also changed
% % % % % % [EPSC_II EPSC_A]=rnd_spike(excitation);
% % % % % % buff_length=length(II_waveform)-EPSC_II_count; % Calculates how the long the buffer is that keeps information about residual conductances from previous EPSC spikes
% % % % % % buffer=II_waveform(EPSC_II_count:length(II_waveform)); % Buffer that keeps track of any residual conductances from previous EPSC spikes
% % % % % % new_waveform=zeros(1,round(EPSC_length)); % 500? pick appropriate #
% % % % % % new_waveform(1)=EPSC_A; % Amplitude of new EPSC spike
% % % % % % if EPSC_shape==1
% % % % % % new_waveform=alpha_func(new_waveform)/10; % Calculate the waveform for the calculated current (pA) clamped @ -94 mV for new EPSC spike
% % % % % % elseif EPSC_shape==2
% % % % % % new_waveform=alpha_func_elongated(new_waveform)/10;
% % % % % % elseif EPSC_shape==3
% % % % % % new_waveform=alpha_func_blunted(new_waveform)/10;
% % % % % % end
% % % % % % II_waveform=new_waveform;
% % % % % % II_waveform(1:buff_length)=II_waveform(1:buff_length)+buffer(1:buff_length); % Superimpose new EPSC spike over any residual conductances from previous spikes
% % % % % % EPSC_II_count=1;
% % % % % % else % No residual conductance from previous EPSC
% % % % % % [EPSC_II EPSC_A]=rnd_spike(excitation);
% % % % % % new_waveform=zeros(1,EPSC_length); % 500? pick appropriate #
% % % % % % new_waveform(1)=EPSC_A; % Amplitude of new EPSC spike
% % % % % % new_waveform=alpha_func(new_waveform)/10; % Calculate the waveform for the calculated current (pA) clamped @ -94 mV for new EPSC spike
% % % % % % II_waveform=new_waveform;
% % % % % % EPSC_II_count=1;
% % % % % % end
% % % % % % end
% Updating g_syn
% if EPSC_II_count < EPSC_length
% gb_syn=EPSC_train(zz)/97*mag_mult;
% plot_syn(zz)=gb_syn;
% else
% gb_syn=0;
% plot_syn(zz)=gb_syn;
% end
% EPSC_II_count=EPSC_II_count+1;
gb_syn=EPSC_train(zz)/10*mag_mult/97;
% Updating currents from each type of channel
Ina_rm=I_na_rm(V(2));
Ik_rm=I_k_rm(V(2));
Ihtk_rm=I_htk_rm(V(2));
Iltk_rm=I_ltk_rm(V(2));
Ih_rm=I_h_rm(V(2));
% Updating currents from leak channel and synaptic conductances
Il=I_l(V(2));
Isyn=I_syn(V(2));
% Calculating the new voltage value given the calculated currents
V(3)=(-2*dt/c*(Isyn+Ina_rm+Ik_rm+Ihtk_rm+Iltk_rm+Ih_rm+Il)+4*V(2)-V(1))*1/3;
% % Spike detection and calculation of Interspike Interval
% ii_count=ii_count+1;
% if V(3) > -30
% if V(3) < V(2) && V(2) > V(1)
% II(spike_count)=ii_count;
% ii_count=1;
% spike_count=spike_count+1;
% end
% end
% Updating new values for voltage buffer
V(1)=V(2);
V(2)=V(3);
VV(zz)=V(3); % for recording voltages
end
%% Spike Detection
II_array=spike_detection(VV);
figure(69)
subplot(3,1,1)
plot(dt:dt:dur,VV); hold on % Plotting Voltages
title('Model Response to Synaptic Excitation')
xlabel('time (ms)')
ylabel('mV')
ylim([-100 100])
subplot(2,1,2)
plot(dt:dt:dur,plot_syn(1:dur/dt)); hold on
title('Synaptic Excitation')
xlabel('time (ms)')
ylabel('mV')
ylim([0 1.5])
% II_array=II(2:length(II));
time=toc;