% 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.1 panels (a)-(e) (current-based synapses 5mV 100 Hz)
%------------------------------------------------------------------------------------

%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
th=20*10^-3;%V Threshold
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; %time discretization vector

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


%%%DiPDE Numerical Solution
%--------------------------

%Current-based synapses

wts_curr = [4.867*10^-3];    % peak of post-synaptic potential in NEST simulations
dwts = ones(size(wts_curr)); 
qwts = ones(size(wts_curr));
qdw_curr=[wts_curr; ones(size(wts_curr)); qwts.*dwts]; 

fr=100;                     %Input frequency in Hz
flag_poisson=1;             %Poisson input   
flag_cond_syn=0;            %implies current-based synapses
neur_type = 'LIF';          

voltage_fraction=1-(1/(exp(1)));  % measure of residual current after spike (set to (1-1/e) since tau_syn = 0.1ms (one time step)
                                  %for current-based simulations

%Initial probability density

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

[cellpop connect]=init_cellpop_no_dep;  %initialize network

%Set parameters for current simulation
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;

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

%outputs to save for figures
t_curr_dipde = t_vec;
dt_curr_dipde = t_vec(2)-t_vec(1);
v_curr_dipde = cellpop(2).v;
prob_unif_curr_dipde = cellpop(2).prob_unif;
fout_curr_dipde = cellpop(2).fout;
fout_curr_up =  cellpop(2).fout_up;
fout_curr_low =  cellpop(2).fout_low;

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

%%%% NEST SIMULATIONS
%--------------------

dv_sim = 0.5*10^-3;
dt_sim = 0.1*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;           %number of neurons used in NEST simulations

dtbin = 2*10^-3;   %bin-size for binning spike-times
tbin = 0:dtbin:tf; 
tplot = 0:dtbin:tbin(end)-dtbin;

%Current-based synapses (load spiketimes and voltage data from NEST simulations)
load IFCE_OGUniMoSprProb_5mV_100Hz_st_nest;
load IFCE_OGUniMoSprProb_5mV_100Hz_v_nest;

f1 = IFCE_OGUniMoSprProb_5mV_100Hz_st_nest;
f1 = [f1(f1(:,2)<10001,1) f1(f1(:,2)<10001,1)]; %to count from first 10000 neurons
f2 = IFCE_OGUniMoSprProb_5mV_100Hz_v_nest;

s_curr = f1(:,1);         % (ms) extracting spike-times
s1_curr = s_curr-1.2;     % subtracting offset due to inherent delay in NEST simulation inputs 
s2_curr = s1_curr(s1_curr<=200); % NEST simulation run for 205 ms to avoid binning problems with last bin
simfrate_curr = (1/nn)*((length(tbin)-1)/tbin(end))*hist(s2_curr*10^-3,tbin); % obtain histogram of spike-times

p_curr = f2(1:(length(tsim1)+12)*nn, 1); % (mV) extracting membrane potential records; 12 accounts for the 1.2 ms offset
p2_curr = reshape(p_curr*10^-3, nn, []); % reshape voltage records obtained from NEST
probsim_curr = (1/nn)*((length(vsim1)-1)/vsim1(end))*hist(p2_curr(:,13:end),vsim1); % histogram for heat-map in Fig. 1(a)

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

%%%%FOKKER_PLANCK
%-----------------

%Drift and diffusion to be input for equivalent Fokker-Planck numerical solution
mu = wts_curr*fr*(1/lk); %(V) drift
sig2 = wts_curr.^2*fr;   % not actual diffusion, defined input for FP numerics
sigma2 = (sig2/2)*(1/lk); %(V^2/s) actual diffusion

%Use the adjusted bin-size obtained from DiPDE as input to FP numerics
dv_curr_dipde = (v_curr_dipde(2) - v_curr_dipde(1));
nbins_new = floor(th/dv_curr_dipde)-1;
dv_new = th/nbins_new;
v_new = 0:dv_new:th;

[v_fp t_fp fout_fp ren_prob_fp ren_prob_fp_full raw_prob_fp dt_fp]=fp_lv_leak_curr_ftcs_new(th,lk,dt0,tf,v_new,p0,mu,sig2);

%Re-binning the FP prob distribution according to original v0 bins
prob_unif_fp_slice = zeros(1,length(v_new));
 for ii = 1:length(t_fp)
     prob_unif_fp_slice = raw_prob_fp(ii,:);
     temp_fp = rebin_uniform_vector(v_new,v0,prob_unif_fp_slice);
     prob_unif_fp0(ii,:) = temp_fp;
     prob_unif_fp0(ii,:)=(temp_fp/(v0(end)-v0(1)))*(length(v0)-1);
 end     

%Save data
save('Fig1data.mat','v0','tsim1','vsim1','probsim_curr','t_curr_dipde','prob_unif_curr_dipde','t_fp','prob_unif_fp0','tplot', 'simfrate_curr','fout_curr_dipde','fout_curr_up','fout_curr_low','fout_fp','dt_curr_dipde','dt_fp')


