% 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. S4 showing  bounds that can obtained for
% non-instantaneous synapses using DiPDE instantaneous synapse formalism. 
%Refer to Methods: Non-instantaneous synapses in the paper

%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;

%%%DiPDE settings
dt0 = 0.1*10^-3; %(ms) time-step used
tau_ref = dt0; %parameter required for wrap_around_frac in DiPDE
tf=200*10^-3; % (s) final time
tvec = 0:dt0:tf; % time_discretization vector
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 

%Initializing parameters
lk=1/(20*10^-3);% (s^-1) leak constant 1/RC
taum = 20*10^-3; % (s) membrane time constant
taus = 1*10^-3;  % (s) synaptic time constant
tmax = (taum*taus)/(taum-taus)*log(taum/taus); % (s) time taken to reach maximum depolarization
tmin = -taum*log(taum/(taum-taus));  % (s) parameter used for obtaining an upper bound on output firing rate

%%%% Possible choices of w in current-based DiPDE numerical solutions

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  
v_max_weak = vmax*((taum)/(taum-taus)); % (V) for the weak upper bound

%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

neur_type = 'LIF'; 

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

%Equating total charge to get an estimate of output firing rate

tau_d = 0*10^-3;
V0_lst = vmax*exp(-tau_d/taum); % (V)

wts_curr = V0_lst;
dwts = ones(size(wts_curr)); 
qwts = ones(size(wts_curr));
flag_cond_syn=0;            % current-based synapses
qdw_curr=[wts_curr; ones(size(wts_curr)); qwts.*dwts]; 
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_taud0 = cellpop(2).fout;
frate_orig_taud0_up = cellpop(2).fout_up;
frate_orig_taud0_low = cellpop(2).fout_low;
t_orig_taud0 = t_vec;

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

%Obtaining an upper bound on output firing rate

tau_d = tmin;
V0_lst = vmax*exp(-tau_d/taum);

wts_curr = V0_lst;
dwts = ones(size(wts_curr)); 
qwts = ones(size(wts_curr));
flag_cond_syn=0;            % current-based synapses
qdw_curr=[wts_curr; ones(size(wts_curr)); qwts.*dwts]; 
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_taudtmin = cellpop(2).fout;
frate_orig_taudtmin_up = cellpop(2).fout_up;
frate_orig_taudtmin_low = cellpop(2).fout_low;
t_orig_taudtmin = t_vec;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Obtaining an intermediate estimate on the output firing rate

tau_d = taus;
V0_lst = vmax*exp(-tau_d/taum); % (V) 

wts_curr = V0_lst;
dwts = ones(size(wts_curr)); 
qwts = ones(size(wts_curr));
flag_cond_syn=0;            % current-based synapses
qdw_curr=[wts_curr; ones(size(wts_curr)); qwts.*dwts]; 
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_taudts = cellpop(2).fout;
frate_orig_taudts_up = cellpop(2).fout_up;
frate_orig_taudts_low = cellpop(2).fout_low;
t_orig_taudts = t_vec;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Equating maximum depolarization to obtain a lower bound on output firing
%rate

tau_d = tmax;
V0_lst = vmax*exp(-tau_d/taum); % (V)

wts_curr = V0_lst;
dwts = ones(size(wts_curr)); 
qwts = ones(size(wts_curr));
flag_cond_syn=0;            % current-based synapses
qdw_curr=[wts_curr; ones(size(wts_curr)); qwts.*dwts]; 
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_taudtmax = cellpop(2).fout;
frate_orig_taudtmax_up = cellpop(2).fout_up;
frate_orig_taudtmax_low = cellpop(2).fout_low;
t_orig_taudtmax = t_vec;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%NEST SIMULATIONS

nn = 10000;

dtsim_b3 = 0.1; %(ms)
td_b3 = 1+2*dtsim_b3; %offset due to inherent delay in NEST simulation inputs 

dtbin_b3 = 30*dtsim_b3*10^-3;
tbin_b3 = 0:dtbin_b3:tf;
tplot_b3 = 0:dtbin_b3:tbin_b3(end);

load IFCE_OGUniMoSprProb_5mV_100Hz_dttref0p1_taus1_correct_st_nest;
  
fb3st = IFCE_OGUniMoSprProb_5mV_100Hz_dttref0p1_taus1_correct_st_nest; 
 
sb3 = fb3st(:,1);         % (ms) extracting spike-times
sb3 = sb3-td_b3;     % subtracting offset due to inherent delay in simulation inputs 
sb3 = sb3(sb3<=200); % simulation run for 205 ms to avoid binning problems with last bin
simfrate_dttref0p1_1 = (1/(nn*dtbin_b3))*hist(sb3*10^-3,tbin_b3); % obtain histogram of spike-times

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

%Save data
save('SuppFigS4.mat','tplot_b3','simfrate_dttref0p1_1','t_orig_taud0','frate_orig_taud0','t_orig_taudtmin','frate_orig_taudtmin','t_orig_taudtmax','frate_orig_taudtmax','t_orig_taudts','frate_orig_taudts');
save('SuppFigS4_up_rates.mat','t_orig_taud0','frate_orig_taud0_up','t_orig_taudtmin','frate_orig_taudtmin_up','t_orig_taudtmax','frate_orig_taudtmax_up','t_orig_taudts','frate_orig_taudts_up');
save('SuppFigS4_low_rates.mat','t_orig_taud0','frate_orig_taud0_low','t_orig_taudtmin','frate_orig_taudtmin_low','t_orig_taudtmax','frate_orig_taudtmax_low','t_orig_taudts','frate_orig_taudts_low');
