% 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 Supp. Fig. S2 showing time-step-independence 
%(0.1,0.5 and 1.0 ms) of DiPDE

clear variables;

%Initializing parameters
lk=1/(20*10^-3);% (s^-1) leak constant 1/RC
taum = 20*10^-3; % (s) membrane time constant
taus = 0.1*10^-3;  % (s) synaptic time constant


tf=200*10^-3; % (s) final time
dv0=0.5*10^-3;% (V) defines max step in potential
th=20*10^-3;% (V) Threshold
v_low = 0*10^-3;% (V) Resting potential
v0=(v_low:dv0:th);% (V) defines the grid for the membrane potential 

%%%% w in curr-based DiPDE

vmax = 5.0*10^-3; % (V) value of (J_syn*tau_syn/C_m) used in simulations 
v_max_dep = vmax*(taus/taum)^(taus/(taum-taus)); % (V) maximum depolarization due to expo decaying current-based synapse  

%%Current-based synapses

wts_curr = v_max_dep; %(V) Setting w in DiPDE to be equal to the maximum depolarization
dwts = ones(size(wts_curr)); 
qwts = ones(size(wts_curr));
qdw_curr=[wts_curr; ones(size(wts_curr)); qwts.*dwts]; 

flag_cond_syn=0;            % current-based synapses

flag_poisson=1; %use Poisson input train 
fr = 100;       %Input frequency (Hz)

neur_type = 'LIF'; 

%Initial probability density

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%DiPDE with dt = 0.1 ms

dt0 = 0.1*10^-3;
tau_ref = dt0;
tvec=0:dt0:tf;% (s) defines grid for time
voltage_fraction=(1-exp(-(tau_ref/taus))); % measure of residual current after spike

[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_curr;
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);                               

%Firing rate with 95% confidence-intervals
frate_orig_dt0p1 = cellpop(2).fout;
frate_orig_dt0p1_up = cellpop(2).fout_up;
frate_orig_dt0p1_low = cellpop(2).fout_low;
t_orig_dt0p1 = t_vec;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%DiPDE with dt = 0.5 ms

dt0 = 0.5*10^-3;
tau_ref = dt0;
tvec=0:dt0:tf;% (s) defines grid for time
voltage_fraction=(1-exp(-(tau_ref/taus))); % measure of residual current after spike

[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_curr;
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);                               

%Firing rate with 95% confidence-intervals
frate_orig_dt0p5 = cellpop(2).fout;
frate_orig_dt0p5_up = cellpop(2).fout_up;
frate_orig_dt0p5_low = cellpop(2).fout_low;
t_orig_dt0p5 = t_vec;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%DiPDE with dt = 1.0 ms

dt0 = 1*10^-3;
tau_ref = dt0;
tvec=0:dt0:tf;  % (s) defines grid for time
voltage_fraction=(1-exp(-(tau_ref/taus)));  % measure of residual current after spike

[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_curr;
cellpop(1).external_fr=[fr*ones(size(tvec))];
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);                               

%Firing rate with 95% confidence-intervals
frate_orig_dt1 = cellpop(2).fout;
frate_orig_dt1_up = cellpop(2).fout_up;
frate_orig_dt1_low = cellpop(2).fout_low;
t_orig_dt1 = t_vec;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Save data
save('SuppFigS2','t_orig_dt0p1','frate_orig_dt0p1','frate_orig_dt0p1_up','frate_orig_dt0p1_low','t_orig_dt0p5','frate_orig_dt0p5','frate_orig_dt0p5_up','frate_orig_dt0p5_low','t_orig_dt1','frate_orig_dt1','frate_orig_dt1_up','frate_orig_dt1_low');
