% 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 Figure 3 (input-output curves for a given
%%%%set of distributions. Please refer to 
%METHODS: Input-Output Curves in the paper for the implementation.

clear variables;

th=20*10^-3;%V Threshold
lk=1/(20*10^-3);%s^-1 leak constant 1/RC

%time discretization
dt=1*10^-3;%s time step
tf=500*10^-3; %s final time
tvec=0:dt:tf;

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

%initial probability distribution
p0=zeros(size(v0));
p0(1)=1;      
voltage_fraction=1;

baserates=0.5:0.25:2;    %Set of input rates to plot (will be multiplied by firing rate)
fractions=0.5:0.1:1;     %Set of fractions of firing rates at which to take equilibrium distributions to identify peak transients


%%%%%Generate synaptic weight distributions
match_drift_only=1;          %distributions get matched for mean input rate and drift
match_4_distributions_1_var;

qdw_gauss=[wts_gauss;ones(size(wts_gauss));qwts_gauss];
qdw_hyp=[wts_hyp; ones(size(wts_hyp)); qwts_hyp]; 
qdw_exp=[wts_exp; ones(size(wts_exp)); qwts_exp]; 
qdw_logn=[wts_logn; ones(size(wts_logn)); qwts_logn]; 

if match_drift_only==1
    qdw_del=[0.001; 1; 1];
    qdw_bidel=[0.0005 0.015;1 1;14/14.5 1-(14/14.5)];
    fr_del=fr_gauss;
    fr_bidel=fr_gauss;
else
    qdw_del=[0.0018; 1; 1];
    fr_del=1/qdw_del(1,1);
    qdw_bidel=[0.0005 0.015;1 1;0.9967 0.0033];
    fr_bidel=1825.33;
end;

%%% DELTA function %%%%%%%%%%%%%%%%%%%%
frates=baserates*fr_del;
qdw=qdw_del;
transient_frates=[0 0 0];
fileID = fopen('delta_fractions.dat','w');
firing_rate_store=cell(length(fractions),length(frates));
tvec_store=cell(length(fractions),length(frates));

%Find peak firing rates
for ii=1:length(fractions)
    startfires=frates*fractions(ii);
    for jj=1:length(frates)
        [cellpop connect]=init_cellpop_no_dep;
        connect(1).qdw=qdw;
        cellpop(1).external_fr=[startfires(jj)*ones(size(tvec)) frates(jj)*ones(size(tvec))];
        [cellpop connect t_vec]=dipde_multi_pop2(cellpop,connect,dt,tf*2);
        transient_frates=[transient_frates;fractions(ii) frates(jj) max(cellpop(2).fout(t_vec>tf))];
        fprintf(fileID,'%d %6.2f %6.2f\n',fractions(ii),frates(jj),max(cellpop(2).fout(t_vec>tf)));
        firing_rate_store{ii,jj}=cellpop(2).fout;
        tvec_store{ii,jj}=t_vec;
    end;
end;
fclose(fileID);
save('delta_all_firing.mat','firing_rate_store','tvec_store'); %delta-function data for Fig.3 panels a),b),c)
save('delta_transient_frates','transient_frates');

%%%%BIMODAL distribution%%%%%%%%%%%%%%%%%%%%
frates=baserates*fr_bidel;
qdw=qdw_bidel;
transient_frates=[0 0 0];
fileID = fopen('bidel_fractions.dat','w');

%Find peak firing rates
for ii=1:length(fractions)
    startfires=frates*fractions(ii);
    for jj=1:length(frates)
        [cellpop connect]=init_cellpop_no_dep;
        connect(1).qdw=qdw;
        cellpop(1).external_fr=[startfires(jj)*ones(size(tvec)) frates(jj)*ones(size(tvec))];
        [cellpop connect t_vec]=dipde_multi_pop2(cellpop,connect,dt,tf*2);
        transient_frates=[transient_frates;fractions(ii) frates(jj) max(cellpop(2).fout(t_vec>tf))];
        fprintf(fileID,'%d %6.2f %6.2f\n',fractions(ii),frates(jj),max(cellpop(2).fout(t_vec>tf)));
        tvec_store{ii,jj}=t_vec;
    end;
end;
fclose(fileID);
save('bidel_transient_frates','transient_frates');

%%%%GAUSSIAN distribution%%%%%%%%%%%%%%%%%%%%
frates=baserates*fr_gauss;
qdw=qdw_gauss;
transient_frates=[0 0 0];
fileID = fopen('gauss_fractions.dat','w');

%Find peak firing rates
for ii=1:length(fractions)
    startfires=frates*fractions(ii);
    for jj=1:length(frates)
        [cellpop connect]=init_cellpop_no_dep;
        connect(1).qdw=qdw;
        cellpop(1).external_fr=[startfires(jj)*ones(size(tvec)) frates(jj)*ones(size(tvec))];
        [cellpop connect t_vec]=dipde_multi_pop2(cellpop,connect,dt,tf*2);
        transient_frates=[transient_frates;fractions(ii) frates(jj) max(cellpop(2).fout(t_vec>tf))];
        fprintf(fileID,'%d %6.2f %6.2f\n',fractions(ii),frates(jj),max(cellpop(2).fout(t_vec>tf)));
        tvec_store{ii,jj}=t_vec;
    end;
end;
fclose(fileID);
save('gauss_transient_frates','transient_frates');

%%%%HYPERBOLIC distribution%%%%%%%%%%%%%%%%%%%%
frates=baserates*fr_hyp;
qdw=qdw_hyp;
transient_frates=[0 0 0];
fileID = fopen('hyper_fractions.dat','w');


%Find peak firing rates
for ii=1:length(fractions)
    startfires=frates*fractions(ii);
    for jj=1:length(frates)
        [cellpop connect]=init_cellpop_no_dep;
        connect(1).qdw=qdw;
        cellpop(1).external_fr=[startfires(jj)*ones(size(tvec)) frates(jj)*ones(size(tvec))];
        [cellpop connect t_vec]=dipde_multi_pop2(cellpop,connect,dt,tf*2);
        transient_frates=[transient_frates;fractions(ii) frates(jj) max(cellpop(2).fout(t_vec>tf))];
        fprintf(fileID,'%d %6.2f %6.2f\n',fractions(ii),frates(jj),max(cellpop(2).fout(t_vec>tf)));
        tvec_store{ii,jj}=t_vec;
    end;
end;
fclose(fileID);
save('hyp_transient_frates','transient_frates');

%%%%EXPONENTIAL distribution%%%%%%%%%%%%%%%%%%%%
frates=baserates*fr_exp;
qdw=qdw_exp;
transient_frates=[0 0 0];
fileID = fopen('exp_fractions.dat','w');

%Find peak firing rates
for ii=1:length(fractions)
    startfires=frates*fractions(ii);
    for jj=1:length(frates)
        fractions(ii)
        frates(jj)
        [cellpop connect]=init_cellpop_no_dep;
        connect(1).qdw=qdw;
        cellpop(1).external_fr=[startfires(jj)*ones(size(tvec)) frates(jj)*ones(size(tvec))];
        [cellpop connect t_vec]=dipde_multi_pop2(cellpop,connect,dt,tf*2);
        transient_frates=[transient_frates;fractions(ii) frates(jj) max(cellpop(2).fout(t_vec>tf))];
        fprintf(fileID,'%d %6.2f %6.2f\n',fractions(ii),frates(jj),max(cellpop(2).fout(t_vec>tf)));
        tvec_store{ii,jj}=t_vec;
    end;
end;
fclose(fileID);
save('exp_transient_frates','transient_frates');

%%%%LOGNORMAL distribution%%%%%%%%%%%%%%%%%%%%
frates=baserates*fr_logn;
qdw=qdw_logn;
transient_frates=[0 0 0];
fileID = fopen('logn_fractions.dat','w');

%Find peak firing rates
for ii=1:length(fractions)
    startfires=frates*fractions(ii);
    for jj=1:length(frates)
        [cellpop connect]=init_cellpop_no_dep;
        connect(1).qdw=qdw;
        cellpop(1).external_fr=[startfires(jj)*ones(size(tvec)) frates(jj)*ones(size(tvec))];
        [cellpop connect t_vec]=dipde_multi_pop2(cellpop,connect,dt,tf*2);
        transient_frates=[transient_frates;fractions(ii) frates(jj) max(cellpop(2).fout(t_vec>tf))];
        fprintf(fileID,'%d %6.2f %6.2f\n',fractions(ii),frates(jj),max(cellpop(2).fout(t_vec>tf)));
        tvec_store{ii,jj}=t_vec;
    end;
end;
fclose(fileID);
save('logn_transient_frates','transient_frates');

%%%%INPUT-OUTPUT CURVES WITH SYNAPTIC DEPRESSION%%%%%%%%%%%%%%%%%%%%

%%%%Set up multi-population run for synaptic depression
%external population/input

[cellpop connect]=init_cellpop;

transient_frates_dep=[-1 0 0];
new_baserates=fliplr(logspace(-1,log10(100),15));
fileID = fopen('depression_rates_025_1msdt.dat','w');

for jj=1:length(new_baserates)
    
    %delta distribution
    [cellpop connect]=init_cellpop;
    connect(1).qdw=qdw_del;
    frates=new_baserates(jj)*fr_del;
    cellpop(1).external_fr=frates*ones(size(tvec));
    [cellpop connect tvec]=dipde_multi_pop2(cellpop,connect,dt,tf);
    transient_frates_dep=[transient_frates_dep;0 frates cellpop(2).fout(end)];
    fprintf(fileID,'%d %6.2f %6.2f\n',0,frates,cellpop(2).fout(end));
    
    %bimodal distribution
    [cellpop connect]=init_cellpop;
    connect(1).qdw=qdw_bidel;
    frates=new_baserates(jj)*fr_bidel;
    cellpop(1).external_fr=frates*ones(size(tvec));
    [cellpop connect tvec]=dipde_multi_pop2(cellpop,connect,dt,tf);
    transient_frates_dep=[transient_frates_dep;1 frates cellpop(2).fout(end)];
    fprintf(fileID,'%d %6.2f %6.2f\n',1,frates,cellpop(2).fout(end));
    
    %gaussian distribution
    [cellpop connect]=init_cellpop;
    connect(1).qdw=qdw_gauss;
    frates=new_baserates(jj)*fr_gauss;
    cellpop(1).external_fr=frates*ones(size(tvec));
    [cellpop connect tvec]=dipde_multi_pop2(cellpop,connect,dt,tf);
    transient_frates_dep=[transient_frates_dep;2 frates cellpop(2).fout(end)];
    fprintf(fileID,'%d %6.2f %6.2f\n',2,frates,cellpop(2).fout(end));
    
    %power law distribution
    [cellpop connect]=init_cellpop;
    connect(1).qdw=qdw_hyp;
    frates=new_baserates(jj)*fr_hyp;
    cellpop(1).external_fr=frates*ones(size(tvec));
    [cellpop connect tvec]=dipde_multi_pop2(cellpop,connect,dt,tf);
    transient_frates_dep=[transient_frates_dep;3 frates cellpop(2).fout(end)];
    fprintf(fileID,'%d %6.2f %6.2f\n',3,frates,cellpop(2).fout(end));
    
    %exponential distribution
    [cellpop connect]=init_cellpop;
    connect(1).qdw=qdw_exp;
    frates=new_baserates(jj)*fr_exp;
    cellpop(1).external_fr=frates*ones(size(tvec));
    [cellpop connect tvec]=dipde_multi_pop2(cellpop,connect,dt,tf);
    transient_frates_dep=[transient_frates_dep;4 frates cellpop(2).fout(end)];
    fprintf(fileID,'%d %6.2f %6.2f\n',4,frates,cellpop(2).fout(end));
    
    %lognormal distribution
    [cellpop connect]=init_cellpop;
    connect(1).qdw=qdw_logn;
    frates=new_baserates(jj)*fr_logn;
    cellpop(1).external_fr=frates*ones(size(tvec));
    [cellpop connect tvec]=dipde_multi_pop2(cellpop,connect,dt,tf);
    transient_frates_dep=[transient_frates_dep;5 frates cellpop(2).fout(end)];
    fprintf(fileID,'%d %6.2f %6.2f\n',5,frates,cellpop(2).fout(end));
    
end;
fclose(fileID);
save('depression_rates_025_1msdt.mat','transient_frates_dep');
