% This file is licensed under the Allen Institute Terms of Use, which can be found at:
%http://www.alleninstitute.org/Media/policies/terms_of_use_content.html

%------------------------------------------------------------------------------------

%Script to generate data for Fig.1h (EIF with current-based synapses)
%-----------------------------------------------------------------------------

%NOTE: Data files for spike-times and voltage records from NEST simulations have 
%also been made available and can be used to reproduce figures exactly.

clear variables;

%Initializing parameters
lk=1/(20*10^-3);% (s^-1) leak constant 1/RC
dt0=0.5*10^-3;% (s) time step
tf=200*10^-3; % (s) final time
tvec=0:dt0:tf;% (s) defines grid for time
dv0=0.5*10^-3;% (V) defines max step in potential
v_low = 0*10^-3;% (V) Resting potential
th=20*10^-3;% (V) Threshold
v0=(v_low:dv0:th);% (V) defines the grid for the membrane potential 


%Initial probability density

%delta function probability starting at 0 mV
p0=zeros(size(v0));
p0(v0==0)=1;

flag_poisson=1; %use Poisson input train 
fr = 100;       %Input frequency
voltage_fraction=1; % measure of residual current after spike (set to 1 since tau_syn = 0.0ms (one time step)
                                      %for conductance-based simulations
neur_type = 'EIF'; 

%%%%DiPDE

%%Conductance-based synapses

vrev_exc = 100*10^-3;         %V   Excitatory reversal potential
vrev_inh = -10*10^-3;         %V   Inhibitory reversal potential
vrevs = [vrev_exc];           %    Only excitatory synapses for Fig.1

wts_cond = [4.805*10^-3];
flag_cond_syn=1;             % conductance-based synapses
dwts = ones(size(wts_cond)); 
qwts = ones(size(wts_cond));
qdw_cond=[wts_cond./(vrevs); (vrevs).*ones(size(wts_cond)); qwts.*dwts];

                                      
[cellpop connect]=init_cellpop_no_dep;

cellpop(1).v0=v0;
cellpop(1).p0=p0;
cellpop(1).external_fr=[fr*ones(size(tvec))];

cellpop(2).v0=v0;
cellpop(2).p0=p0;
cellpop(2).external_fr=zeros(size(tvec));

connect(1).qdw=qdw_cond;
connect(1).flag_cond_syn=flag_cond_syn;
connect(1).flag_poisson=flag_poisson;
connect(1).voltage_frac=voltage_fraction;

[cellpop connect t_vec]=dipde_multi_pop_fig1(cellpop,connect,dt0,tf,neur_type);

t_cond_dipde_eif = t_vec;
dt_cond_dipde_eif = t_vec(2)-t_vec(1);
v_cond_dipde_eif = cellpop(2).v;
prob_unif_cond_dipde_eif = cellpop(2).prob_unif;
fout_cond_dipde_eif = cellpop(2).fout;
fout_cond_up_eif =  cellpop(2).fout_up;
fout_cond_low_eif =  cellpop(2).fout_low;

%--------------------------------------------------------------------------%
  
%%%%SIMULATIONS
 
dv_sim = 0.5*10^-3;
dt_sim = 0.05*10^-3;
vsim1 = (v_low+dv_sim/2:dv_sim:th-dv_sim/2);    %voltage vector to be used for simulation heat-maps
tsim1 = 0:dt_sim:tf;
nn=10000;

dtbin = 2*10^-3;
tbin = 0:dtbin:tf;
tplot = 0:dtbin:tbin(end)-dtbin;

%Conductance-based synapses
load EIFCO_5mV100HzFin1920p05_st_nest;

f3 = EIFCO_5mV100HzFin1920p05_st_nest;
f3 = [f3(f3(:,2)<10001,1) f3(f3(:,2)<10001,1)]; %to count from first 10000 neurons

s_cond = f3(:,1);
s1_cond = s_cond-1.2;
s2_cond= s1_cond(s1_cond<=200);
simfrate_cond = (1/nn)*((length(tbin)-1)/tbin(end))*hist(s2_cond*10^-3,tbin);

% %--------------------------------------------------------------------------%

save('Fig1h_EIF_data.mat','v0','tsim1','vsim1','t_cond_dipde_eif','prob_unif_cond_dipde_eif','tplot', 'simfrate_cond','fout_cond_dipde_eif','fout_cond_up_eif','fout_cond_low_eif','dt_cond_dipde_eif')

