%%% This script is the top level script for controlling the modeling
%%% environment.  Note that some loops are now no longer working but are
%%% left in place temporarily while the top level code is cleaned up for
%%% public sharing. Original code written by Ariel Hight and edited by RK (7/24/2018)
%%% Code modified for IH conductances by C.M. Ventura (2017-2018)

clear all
close all

tic


%Set some global variables
%% global a_k_rm b_k_rm c_k_rm n_k_rm m_na_rm h_na_rm n_htk_rm p_htk_rm w_ltk_rm z_ltk_rm r_h_rm
global gb_k_rm gb_na_rm gb_ltk_rm gb_htk_rm gb_h_rm gl 
global Ena El Ek Eh Esyn
global dt dur VV I_print plot_syn time spike_rate_array spike_count EPSC_shape Ih_shift Isyn  %(added Isyn on 3/9/2016)



%Set the timing variables
%% time and time steps
dur=2500;                       % time duration (ms)
dt=0.01;                        % delta t (ms)
start=500;                      % time at onset of current clamp (ms)
c_duration=1500;                     % length of current clamp (ms)
%c_duration=5;
I_print=1;                      % Plot the individual ionic currents

nt=1;                           % time step (current t = nt*dt)

% calculate start and stop #'s
start=start/dt;
stop=c_duration/dt+start;


%% Experimental Conditions (e.g. reversal potential, etc.)

% Intristic Membrane Properties
c=0.9;                          % Membrane capacitance (microF/cm^2)

% Ionic Battery Values
Ek=-80.78;                         %K battery (mV)
% Ek=-60;                         %K battery (mV)
Ena=81.27;                         %Na battery (mV)
El=-65;                         %Leak battery (mV)
%El=0;
Eh=-46;                         %Ionic battery (mV)
Esyn=3;                         %Synaptic battery (mV)


%% Neuron Type
neuron_type=3;              %1 = Sustained;   2 = Transient

        if neuron_type==1         % Sustained
            V(1)=-60.9766;      %H&K               % ~Resting Potential 
            %V(1) = -54;
            %V(1) = -60.2;
            gb_k_rm=0;                      %gK bar (mS/cm^2)                           
            gb_na_rm=20;                    %gNa bar (mS/cm^2) %used 13 mS/cm2 in Hight and Kalluri and the larger 20 mS/cm2 in Ventura and Kalluri                         
            gb_htk_rm=2.8;                  %gK bar (high threshold K) (mS/cm^2)        
            gb_ltk_rm=0.0;                 %gK bar (low threshold K) (mS/cm^2)        
            gb_h_rm = 0.1;                  %V(1)=-60.2  %gh bar (hyper-pol act. cation) (mS/cm^2)   
            gl=0.03;                        %gLeak (mS/cm^2)
        elseif neuron_type== 2    % Transient
            V(1)=-72.7030;                     % ~Resting Potential 
            gb_k_rm=0;                      %gK bar (mS/cm^2)                           
            gb_na_rm=20;                    %gNa bar (mS/cm^2)                          
            gb_htk_rm=2.8;                %gK bar (high threshold K) (mS/cm^2)       
            gb_ltk_rm=1.1;                   %gK bar (low threshold K) (mS/cm^2)        
            gb_h_rm = 0.1;
            gl=0.03;                        %gLeak (mS/cm^2)
        elseif neuron_type==3     %Custom cell parameters        
            V(1)=-60.97;                     % ~Resting Potential (going to be variable, depending on composition of conductances) 
            gb_k_rm=0;                      %gK bar (mS/cm^2)                           
            gb_na_rm=20;                    %gNa bar (mS/cm^2)                         
            gb_htk_rm=2.8;                %gK bar (high threshold K) (mS/cm^2)        
            gb_ltk_rm=0.0;                   %gK bar (low threshold K) (mS/cm^2)         
            %gb_h_rm=0.0;                    %gh bar (hyper-pol act. cation) (mS/cm^2)   
            gb_h_rm = 0.91;
            gl=0.03;                        %gLeak (mS/cm^2)
        end

        
        
%% EPSC_shape
EPSC_shape=1;               %1 = alpha(0.4)     2 = alpha + elongation      3 = alpha blunted % The shape is defined in generat_EPSC_train.m
trial=neuron_type % appears to be redundant with neuron_type, edited by RK on 3/18/2013

%% Stimulation
epsc_amp_94=0;

% Conditions for Simulations ==> 
    % CS: The stimulation conditions and desired responses
    % Trial: The type of cell and its respective conductances

    cs=2;
 
    %% cs = 2 (CC voltage response, for step evoked currents); cs = 3 (EPSC II, for collecting II times); cs = 4(CC II)
    %see lines 214 onwards for controlling the batch processing for
    %step-evoked current-clamp simulations
    if cs==2
        I=zeros(1,dur/dt);           % Array representing current clamp
        I(start:stop)=5;                % current clamp (x10 pA)
     end
     
     
     
     %% EPSC_RESPONSE %%%
     if cs==3    %Calculating EPSC response
            % Batch for Intervent Interval Calculations
           b_switch=2;     % Batch on (change excitation variable) = 1; Batch on (change magnitude variable) = 2; Batch on (EPSC duration) = 3; Batch on (Dynamic Range Finder) = 4; Batch off = 0
           mag_array=10.^(-0.8:0.2:0.2);
           excite_array=[15 20 30 40 50 60 70 80 90 100 120 140 180 200 220 240 15 20 30 40 50 60 70 80 90 100 120 140 180 200 220 240 15 20 30 40 50 60 70 80 90 100 120 140 180 200 220 240];
           excite_array =[1000,500];
           
           if b_switch==1 %titrate the epsc interval
                mag_mult=0.1;
                b_length=length(excite_array);
                II_is=zeros(b_length,dur/2);    %Preallocation
                for batch=1:b_length
                    excitation=1000/excite_array(batch);                                          % mean EPSC intervent interval (ms)
                    rate_dispay=excite_array(batch)
                    temporary=EPSC_excitation_response(V(1),c,excitation,mag_mult);         %Iterspike Interval (ms)
%                     II_is(batch,1:length(temporary))=temporary;
                    %Save Data
%                         save_trial_data(neuron_type,temporary,V,mag_mult,excitation)
%                         append_excel_data(temporary,neuron_type,excitation,mag_mult)
                    %Save Data
                end
                
           elseif b_switch==2 %titrate the EPSC magnitute, keep excitation rate fixed.
                excitation =3;
                b_length=length(mag_array);
                II_is=zeros(b_length,dur/2);
                for batch=1:b_length
                   mag_mult=mag_array(batch);
                    temporary=EPSC_excitation_response(V(1),c,excitation,mag_mult);         %Iterspike Interval (ms)
                    II_is(batch,1:length(temporary))=temporary;
                    
                    %Save Data
%                         save_trial_data(trial,temporary,V,mag_mult,excitation)
%                         append_excel_data(temporary,trial,excitation,mag_mult)
                    %Save Data
                end
                
           elseif b_switch==3
                excitation=.4;
                mag_mult=1;
                b_length=21;
                II_is=zeros(b_length,dur/2);
                for batch=1:b_length
%                     spike_dur=1*batch;
                    temporary=EPSC_excitation_response(V(1),c,excitation,mag_mult)*dt;         %Iterspike Interval (ms)
                    II_is(batch,1:length(temporary))=temporary;
                    
                    %Save Data
%                         save_trial_data(trial,temporary,V,mag_mult,excitation)
%                         append_excel_data(temporary,trial,excitation,mag_mult)
%                     %Save Data
                    
                end
           elseif b_switch==0
                excitation=3;          % mean EPSC intervent interval (ms)
                mag_mult=1;             
%                 spike_dur=0.4;         % this does not refer to the length of the EPSC spike, but rather represents the alpha value (e.g. alpha = 1.4 corresponds to spike duration of ~ 17.5 ms)
                II_is=EPSC_excitation_response(V(1),c,excitation,mag_mult)*dt;       %Iterspike Interval (ms)
                
                %Save Data
%                     save_trial_data(trial,temporary,V,mag_mult,excitation)
%                     append_excel_data(temporary,trial,excitation,mag_mult)
                %Save Data
            
           elseif b_switch==4 % This code bit detects the boundaries of the response area.  Fig 9 in Hight and Kalluri 2018 was generated by repeated running this and saving data in an excel spreadsheet for subsequent plotting.
%                 spike_dur=0.4;
                array=boundary_detection_array(1);
                initial=[12 18];                % x,y coordinate
                x_direction=-1;             % -1 == left, 1== right
                y_direction=1;              % -1 == bottom, 1== top
                direction=[x_direction y_direction];
                boundary_detection(array,initial,direction,V(1),c,trial);
                for n=1:length(array(1,:,1))
                    visual(:,length(array(1,:,1))-n+1)=spike_rate_array(:,n);
                    visual_spike_count(:,length(array(1,:,1))-n+1)=spike_count(:,n);
                    visual_array_epsc(:,length(array(1,:,1))-n+1)=array(:,n,1);
                    visual_array_mag(:,length(array(1,:,1))-n+1)=array(:,n,2);
                end
                for n=1:length(array(1,:,1))
                    for m=1:length(array(:,1,1))
                        if visual_spike_count(m,n)==0
                            visual_trimmed(m,n)=0;
                        elseif visual_spike_count(m,n) < 0
                            visual_trimmed(m,n)=visual_spike_count(m,n);
                        else
                            visual_trimmed(m,n)=visual(m,n);
                        end
                    end
                end
                visual=visual';
                visual_spike_count=visual_spike_count';
                visual_trimmed=visual_trimmed';
           end            
     elseif cs == 4
%             t=dt:dur/dt;                    % array of time intervals
%             start=start/dt;c
%             stop=start+c_duration/dt;
%             
            II_cc=CC_response();
        


            
            
            
            
            
            
            
            %%%HERE BEGINS CONTROLS FOR BATCH SIMULATIONS in SIMPLE CURRENT
            %%%CLAMP
            %% CC_RESPONSE; Mostly for step-evoked response %%%%
     elseif cs == 1 || cs == 2
            b_switch=2;             % b_switch = 1              Batch on (change current clamp mag.); 
                                    % b_switch = 2              Batch on (change specific conductance); 
                                    % b_switch = 0              Batch off 
            plot_switch=3;          % plot_switch = 0           (independent plots); 
                                    % plot_switch = 1           (all on one plot, same figure); 
                                    % plot_switch = 2           (all on 1 figure, all subplots); 
                                    % plot_switch = 3           (no plots)
            v_rest_control=0;       % v_rest_control = 1        Control the resting potential using constant current
                                    % v_rest_control = 0        No control over resting potential
            threshold_switch=1;     % threshold_switch = 1   	Turn on code (detect the voltage and current threshold)
                                    % threshold_switch = 0      Turn off code (detecting the voltage and current threshold)
            save_switch=1;          % save_switch = 1           Save the responses & conditions of the stimulation
                                    % save switch = 0           don't save the responses
            
            % Batch for Current Clamp
            if b_switch==1 || b_switch==2
               %cc_mag_array=[-10,-5,-2,-1,0,1,2,5,7,10];                                       % CC magnitude array 1
               cc_mag_array = [5.5];
               %b_length=length(cc_mag_array); %for current clamp series
           
                
               gna_steps= ones(1,20)*20;
               gltk_steps= linspace(0.2,0.80,20); %titrating gltk
               gh_steps= ones(1,20)*0.91;
               b_length=length(gna_steps); %for conductances series
               V=ones(b_length,dur/dt)*V(1);
               legend_value=ones(1,b_length);
               %gna_steps=[8 13 18 23 28];
               %gltk_steps=[0 0 0 0 0];
               ss=zeros(1,length(cc_mag_array));       % preallocating the Vhold array
               I_threshold=zeros(1,length(cc_mag_array));
               for batch=1:b_length
                    if b_switch==1
                        cc_mag=cc_mag_array(batch);
                        I(start:stop)=cc_mag;
                        legend_value(batch)=cc_mag;
                    elseif b_switch==2
                        cc_mag=5.5; %fixed array
                        I(start:stop)=cc_mag;
                        % Select which conductance to chance (use ctl. R or
                        % T to toggle OFF or ON respectively)        
                         gb_na_rm=gna_steps(batch);
                         gb_h_rm=gh_steps(batch);
%                        gb_htk_rm=batch*1.5-1.5;
                         gb_ltk_rm=gltk_steps(batch);
                         batch;
                    end
                    if v_rest_control==1
                        if threshold_switch==1
                            % Where Threshold code starts
                            dur=2000;
                            I(1:dur/dt)=cc_mag_array(batch);
                            start2=400;
                            cc_mag_array_2=0.54;                               % CC magnitude Array 2
                            v_threshold=zeros(length(cc_mag_array_2),5);    % preallocating array for recording V thresholds with respect to the dV/dt threshold (first column will be current used, i.e. cc magnitude array 2)
                            I_thres_record(1:4)=0;                          % preallocating array for recording I thresholds with respect to the dV/dt threshold
                            dvdt=zeros(1,length(350/dt:(dur-2)/dt));     % preallocating array for dV/dt array
                            d2vdt2=zeros(1,length(350/dt:(dur-2)/dt));   % preallocating array for d2V/dt2 array
                            v_phase=zeros(1,length(350/dt:(dur-2)/dt));  % preallocating the voltage array corresponding to the dV/dt and d2V/dt2 arrays
                            v_max=zeros(1,length(cc_mag_array_2)+1);        % preallocating the v_max array (+1 so that v_max(1) can always be zero in the case that the first stimulation generates an action potential)
                            v_max(1)=-999;
                            % Preallocating various arrays depending on the
                            % threshold detection (e.g. sustained v.
                            % transient
                            if trial==1
                                v_record=zeros(1,2);
                                v_threshold=zeros(length(cc_mag_array_2),3);
                            else
                                v_record=zeros(1,4);
                                v_theshold=zeros(length(cc_mag_array_2),5);
                                I_thres_record=zeros(1,4);
                            end
                            for batch_2=1:length(cc_mag_array_2)
                                cc_mag_2=cc_mag_array_2(batch_2);
                                I(start2/dt:dur/dt)=cc_mag+cc_mag_2;
                                [V,Ioutput]=single_compart_second_rm(V(1),I,c,cs,trial,epsc_amp_94);
                                ss(batch)=V(400/dt-1);
                                v_hold=ss(batch);
                                if trial==1
                                else
                                    v_record(1:length(v_record))=0;             % counter for recording the voltage (i.e. once the voltage threshold is recorded, the loop stops looking for the voltage threshold)
                                end
                                v_threshold(batch_2,1)=cc_mag_2*10;         % Current step used (for determining V_Threshold)
                                for rr=0:length(V(350/dt:(dur-2)/dt))
                                    dvdt(rr+1)=(-V(350/dt+2+(rr+1))+8*V(350/dt+1+(rr+1))-8*V(350/dt-1+(rr+1))+V(350/dt-2+(rr+1)))/(12*dt);
                                    d2vdt2(rr+1)=(-V(350/dt-2+(rr+1))+16*V(350/dt-1+(rr+1))-30*V(350/dt+(rr+1))+16*V(350/dt+1+(rr+1))-V(350/dt+2+(rr+1)))/(12*dt^2);
                                    v_phase(rr+1)=V(350/dt+rr);
                                    if trial == 1
%                                         if v_phase(rr+1) > -60 && v_phase(rr+1) < -50 && d2vdt2(rr+1) > 0 && v_record(1)==0
%                                             potential_v_threshold=v_phase(rr+1);
%                                             v_record(1)=1;
                                        if v_phase(rr+1) > -54 && v_phase(rr+1) < -50 && d2vdt2(rr+1) > 0 && v_record(2)==0
%                                             v_threshold(batch_2,2)=potential_v_threshold;
                                            v_threshold(batch_2,2)=v_phase(rr+1);
                                            v_threshold(batch_2,3)=v_max(batch_2);
                                            I_threshold(batch)=cc_mag_2*10;
                                            v_record(2)=1;
                                        end
                                    else
                                        if dvdt(rr+1) >= 7.5 && v_record(1)==0;
                                            v_threshold(batch_2,1+1)=v_phase(rr+1);   % V_Threshold(dv/dt=0.00075)
                                            v_record(1)=1;
                                        end
                                        if dvdt(rr+1) >= 7.5 && I_thres_record(1)==0;
                                            I_threshold(batch,1)=cc_mag_2*10;
                                            I_thres_record(1)=1;
                                        end
                                        if dvdt(rr+1) >= 10 && v_record(2)==0;
                                            v_threshold(batch_2,2+1)=v_phase(rr+1);   % V_Threshold(dv/dt=0.001)
                                            v_record(2)=1;
                                        end
                                        if dvdt(rr+1) >= 10 && I_thres_record(2)==0;
                                            I_threshold(batch,2)=cc_mag_2*10;
                                            I_thres_record(2)=1;
                                        end
                                        if dvdt(rr+1) >= 12.5 && v_record(3)==0;
                                            v_threshold(batch_2,3+1)=v_phase(rr+1);   % V_Threshold(dv/dt=0.00125)
                                            v_record(3)=1;
                                        end
                                        if dvdt(rr+1) >= 12.5 && I_thres_record(3)==0;
                                            I_threshold(batch,3)=cc_mag_2*10;
                                            I_thres_record(3)=1;
                                        end
                                        if dvdt(rr+1) >= 15 && v_record(4)==0;
                                            v_threshold(batch_2,4+1)=v_phase(rr+1);   % V_Threshold(dv/dt=0.0015)
                                            v_record(4)=1;
                                        end
                                        if dvdt(rr+1) >= 15 && I_thres_record(4)==0;
                                            I_threshold(batch,4)=cc_mag_2*10;
                                            I_thres_record(4)=1;
                                        end
                                    end
                                end
                                v_max(batch_2+1)=max(v_phase);
                                if plot_switch==0
                                    figure(batch_2)
                                    plot(dt:dt:dur,V)
                                    ylim([-80 80])
                                    xlabel('time (ms)')
                                    ylabel('Vm (mV)')
                                    figure(batch_2+100)
                                    plot(v_phase,dvdt)
                                    xlabel('Vm (mV)')
                                    ylabel('dV/dt')
                                    color_code=['r' 'y' 'g' 'c' 'm' 'b' 'k' 'k' 'k' 'k' 'k' 'k' 'k'];
                                    figure(200)
                                    hold on
                                    plot(v_phase,d2vdt2,color_code(batch_2))
                                    xlabel('Vm (mV)')
                                    ylabel('d2V/dt2')
                                    hold off
                                    figure(300)
                                    hold on
                                    plot(v_phase,dvdt,color_code(batch_2))
                                    xlabel('Vm (mV)')
                                    ylabel('dV/dt')
                                    hold off
                                elseif plot_switch==1
                                elseif plot_switch==2
                                else   plot_switch==3 %(No Plots)
                                end
                                if save_switch==1
                                    if trial ==1
                                        type='Sustained';
                                    else
                                        type='Transient';
                                    end
                                    cd(['C:\Users\Ariel\Desktop\House\Modeling\Results\' datestr(now,'yy') '_' datestr(now,'mm') '\' datestr(now,'yy') '_' datestr(now,'mm') '_' datestr(now,'dd') '\'])
                                    save(['Threshold_' type '_Ihold_' num2str(floor((cc_mag*10))) '_' num2str(round(10*((cc_mag*10)-floor((cc_mag*10))-1/10*(10*(cc_mag*10)-floor(10*(cc_mag*10)))))) num2str(round(10*(10*(cc_mag*10)-floor(10*(cc_mag*10))-1/10*(100*(cc_mag*10)-floor(100*(cc_mag*10)))))) '_Iexcite_' num2str(floor((cc_mag_2*10))) '_' num2str(round(10*((cc_mag_2*10)-floor((cc_mag_2*10))-1/10*(10*(cc_mag_2*10)-floor(10*(cc_mag_2*10)))))) num2str(round(10*(10*(cc_mag_2*10)-floor(10*(cc_mag_2*10))-1/10*(100*(cc_mag_2*10)-floor(100*(cc_mag_2*10)))))) '_timestamp_' datestr(now,'HH_MM_ss_AM')], 'V', 'I', 'v_phase', 'gb_k_rm', 'gb_na_rm', 'gb_htk_rm', 'gb_ltk_rm', 'gb_h_rm', 'gl','dur','v_hold','v_max')
                                    figure(batch_2+100)
                                    saveas(gcf,['Threshold_' type '_Ihold_' num2str(floor((cc_mag*10))) '_' num2str(round(10*((cc_mag*10)-floor((cc_mag*10))-1/10*(10*(cc_mag*10)-floor(10*(cc_mag*10)))))) num2str(round(10*(10*(cc_mag*10)-floor(10*(cc_mag*10))-1/10*(100*(cc_mag*10)-floor(100*(cc_mag*10)))))) '_Iexcite_' num2str(floor((cc_mag_2*10))) '_' num2str(round(10*((cc_mag_2*10)-floor((cc_mag_2*10))-1/10*(10*(cc_mag_2*10)-floor(10*(cc_mag_2*10)))))) num2str(round(10*(10*(cc_mag_2*10)-floor(10*(cc_mag_2*10))-1/10*(100*(cc_mag_2*10)-floor(100*(cc_mag_2*10)))))) '_timestamp_' datestr(now,'HH_MM_ss_AM')])
                                else
                                end
                            end
                        else
                            I(1:dur/dt)=cc_mag_array(batch);
                            V=single_compart_second_rm(V(1),I,c,cs,trial,epsc_amp_94);
                            if plot_switch==0
                                figure(batch)
                                plot(dt:dt:dur,V)
                                xlabel('time (ms)')
                                ylabel('Vm (mV)')
                            elseif plot_switch==1
                            elseif plot_switch==2
                            else        % Plot_switch==3 (No Plot)
                            end
                            ylim([-80 80])
                            ss(batch)=V(400/dt-1);
                        end                        
                        
                    else
                        [V,Ioutput]=single_compart_second_rm(V(1),I,c,cs,trial,epsc_amp_94);
                        if plot_switch==0
                            figure(batch)
                            plot(dt:dt:dur,V)
                        elseif plot_switch==1
                        elseif plot_switch==2
                        else   plot_switch==3; %(No Plot)
                        end
                        ylim([-80 80])
                    end
                        
%                     V=single_compart_second_rm(V(1),t,I,c,cs,trial,epsc_amp_94);
%                     ss(batch)=V(300/dt);                    
                    
                    if plot_switch==0
                        figure(400)
%                         figure(batch)
                        plot(dt:dt:dur,V)
                        ylabel('voltage (mV)')
                        xlabel('time (ms)')
%                         title(['Response to Current Clamp Step at' num2str(cc_mag*10) 'pA'])
                        title(['gna=' num2str(gb_na_rm) ' mS/cm2; ghtk=' num2str(gb_htk_rm) ' mS/cm2; gltk=' num2str(gb_ltk_rm) ' mS/cm2'])
                        ylim([-80 80])
                        
%                         % SAVE CODE
%                         cd(['C:\Users\Ariel\Desktop\House\Modeling\Results\' datestr(now,'yy') '_' datestr(now,'mm') '\' datestr(now,'yy') '_' datestr(now,'mm') '_' datestr(now,'dd') '\'])
%                         if trial ==1
%                             type='Sustained';
%                         else
%                             type='Transient';
%                         end
%                         title([type ' Neuron Response to Current Clamp Step at ' num2str(cc_mag*10) ' pA'])
%                         title(['gna=' num2str(gb_na_rm) ' mS/cm2; ghtk=' num2str(gb_htk_rm) ' mS/cm2; gltk=' num2str(gb_ltk_rm) ' mS/cm2'])
%                         save(['CCresponse_' type '_ccMag_' num2str(cc_mag*10) '_gNa_' num2str(floor(gb_na_rm)) '_' num2str(floor(10*(gb_na_rm-floor(gb_na_rm)))) '_gKlt_' num2str(floor(gb_ltk_rm)) '_' num2str(floor(10*(gb_ltk_rm-floor(gb_ltk_rm)))) num2str(floor(100*(gb_ltk_rm-floor(gb_ltk_rm)-1/10*(floor(10*(gb_ltk_rm-floor(gb_ltk_rm))))))) '_timestamp_' datestr(now,'HH_MM_ss_AM')], 'VV', 'V', 'plot_syn', 'gb_k_rm', 'gb_na_rm', 'gb_htk_rm', 'gb_ltk_rm', 'gb_h_rm', 'gl','dur','cc_mag') 
%                         saveas(gcf,['CCresponse_' type '_ccMag_' num2str(cc_mag*10) '_gNa_' num2str(floor(gb_na_rm)) '_' num2str(floor(10*(gb_na_rm-floor(gb_na_rm)))) '_gKlt_' num2str(floor(gb_ltk_rm)) '_' num2str(floor(10*(gb_ltk_rm-floor(gb_ltk_rm)))) num2str(floor(100*(gb_ltk_rm-floor(gb_ltk_rm)-1/10*(floor(10*(gb_ltk_rm-floor(gb_ltk_rm))))))) '_tstamp_' datestr(now,'HH_MM_ss_AM')])
%                         % SAVE CODE
                        
                    elseif plot_switch==2
                        subplot(2,ceil(b_length/2),batch)
                        plot(dt:dt:dur,V)
                        ylim([-80 40])
                        ylabel('voltage (mV)')
                        xlabel('time (ms)')
                        title(['Response to Current Clamp Step at ',num2str(10*I(start)),' pA'])
                    elseif plot_switch==1
                        color_code=['r' 'y' 'g' 'c' 'm' 'b' 'k' 'k' 'k' 'k' 'k' 'k' 'k''k' 'k' 'k' 'k' 'k' 'k'];
                        %color_code=['b' 'r' 'g' 'g'];
                        figure(20)
                        plot(dt:dt:dur,V,[color_code(1) '-'],'Linewidth',4)
                        %plot(dt:dt:dur,V,[color_code(batch) '-'],'Linewidth',4)
                        %Add the Rm and Vrest here
                        Vrest(batch) = V(end);
                        Rm(batch) = (V(end)-V(40000))./abs(I(start).*10);
                        ylabel('voltage (mV)')
                        xlabel('time (ms)')
                        title('Response to Current Clamp Step')
                        hold on
                    end
                  
                    %% Create a structure with important parameters to easily output.
                    Outputdata.Vsave(batch,:)=V; %RK 5/14/2015
                    Outputdata.gNasave(batch,:)=gb_na_rm; %RK 8/29/2017
                    Outputdata.gKLsave(batch,:)=gb_ltk_rm; %RK 8/29/2017
                    Outputdata.gHsave(batch,:)=gb_h_rm; %RK 8/29/2017
                    Outputdata.cc_mag_save(batch,:)=cc_mag_array; %RK 8/29/2017
                    Outputdata.time = (0:length(V)-1)*dt;
                    Outputdata.Ioutput(batch)=Ioutput;
                    Outputdata.Ihshift(batch,:)=Ih_shift;
                end
                if b_length == 8 && plot_switch==2
                    legend([num2str(legend_value(1)*10),' pA'],[num2str(legend_value(2)*10),' pA'],[num2str(legend_value(3)*10),' pA'],[num2str(legend_value(4)*10),' pA'],[num2str(legend_value(5)*10),' pA'],[num2str(legend_value(6)*10),' pA'],[num2str(legend_value(7)*10),' pA'],[num2str(legend_value(8)*10),' pA'])
                end
%                 if b_length < 18
%                     legend([num2str(legend_value(1)*10),' pA'],[num2str(legend_value(2)*10),' pA'],[num2str(legend_value(3)*10),' pA'],[num2str(legend_value(4)*10),' pA'],[num2str(legend_value(5)*10),' pA'],[num2str(legend_value(6)*10),' pA'],[num2str(legend_value(7)*10),' pA'],[num2str(legend_value(8)*10),' pA'],[num2str(legend_value(9)*10),' pA'],[num2str(legend_value(10)*10),' pA'],[num2str(legend_value(11)*10),' pA'],[num2str(legend_value(12)*10),' pA'],[num2str(legend_value(13)*10),' pA'],[num2str(legend_value(14)*10),' pA'],[num2str(legend_value(15)*10),' pA'],[num2str(legend_value(16)*10),' pA'],[num2str(legend_value(17)*10),' pA'],[num2str(legend_value(18)*10),' pA'])
%                 end
                hold off
% % % % %             elseif b_switch == 2
% % % % %                 figure(230)
% % % % %                 b_length=7;
% % % % %                 V=ones(b_length,dur/dt)*V(1);
% % % % %                 legend_value=ones(1,b_length);
% % % % %                 for batch=1:b_length
% % % % %                     gb_na_rm=1*batch-1;
% % % % %                     legend_value(batch)=gb_na_rm;
% % % % %                     V(batch,:)=single_compart_second_rm(V(1),t,I,c,cs,trial,epsc_amp_94);
% % % % %                     color_code=['r' 'y' 'g' 'c' 'm' 'b' 'k' 'k' 'k' 'k'];
% % % % %                     if b_length < 10
% % % % %                         plot(dt:dt:dur,V(batch,:),color_code(batch))
% % % % %                         ylabel('voltage (mV)')
% % % % %                         xlabel('time (ms)')
% % % % %                         if trial == 1
% % % % %                             title(['Sustained Neuron Response to Current Clamp Step at ',num2str(10*I(start)),' pA'])
% % % % %                         elseif trial == 2
% % % % %                             title(['Transient Neuron Response to Current Clamp Step at ',num2str(10*I(start)),' pA'])
% % % % %                         else
% % % % %                             title(['Custom Neuron Response to Current Clamp Step at ',num2str(10*I(start)),' pA'])
% % % % %                         end
% % % % %                         hold on
% % % % %                     else
% % % % %                         subplot(2,ceil(b_length/2),batch)
% % % % %                         plot(dt:dt:dur,V(batch,:))
% % % % %                         ylim([-80 40])
% % % % %                         ylabel('voltage (mV)')
% % % % %                         xlabel('time (ms)')
% % % % %                         title(['Response to Current Clamp Step at ',num2str(10*I(start)),' pA'])
% % % % %                     end
% % % % %                 end
% % % % %                 if b_length == 9
% % % % %                     legend([num2str(legend_value(1)),'mS/cm2'],[num2str(legend_value(2)),' mS/cm2'],[num2str(legend_value(3)),' mS/cm2'],[num2str(legend_value(4)),' mS/cm2'],[num2str(legend_value(5)),' mS/cm2'],[num2str(legend_value(6)),' mS/cm2'],[num2str(legend_value(7)),' mS/cm2'],[num2str(legend_value(8)),' mS/cm2'],[num2str(legend_value(9)),' mS/cm2'])
% % % % %                 end
% % % % %                 hold off
            else        % b_switch==0
                V=single_compart_second_rm(V(1),I,c,cs,trial,epsc_amp_94);
                % Plotting the responses
                figure(trial+cs)
                subplot(2,1,1)
                plot(dt:dt:dur,V)
                ylabel('voltage (mV)')
                xlabel('time (ms)')
                if trial==1 && cs==1
                    title('Sustained Neuron response to EPSP Spikes')
                    subplot(2,1,2)
                    plot(dt:dt:dur,epsc_amp_94(1:dur/dt)*10)
                    title('Synaptic Activity')
                    ylabel('Current (pA) at Voltage Clamp = 94 mV')
                    xlabel('time(ms)')
                elseif trial ==1 && cs==2
                    title(['Sustained Neuron response to Current Clamp step at ',num2str(10*I(start)),' pA'])
                    subplot(2,1,2)
                    plot(dt:dt:dur,I(1:dur/dt)*10)
                    title('Current Clamp')
                    ylabel('current (pA)')
                    xlabel('time(ms)')
                elseif trial ==2 && cs==1
                    title('Transient Neuron response to EPSP Spikes')
                    subplot(2,1,2)
                    plot(dt:dt:dur,epsc_amp_94(1:dur/dt)*10)
                    title('Synaptic Activity')
                    ylabel('Current (pA) at Voltage Clamp = 94 mV')
                    xlabel('time(ms)')
                elseif trial ==2 && cs==2
                    title(['Transient Neuron response to Current Clamp at ',num2str(10*I(start)),' pA'])
                    subplot(2,1,2)
                    plot(dt:dt:dur,I(1:dur/dt)*10)
                    title('Current Clamp')
                    ylabel('current (pA)')
                    xlabel('time(ms)')
                elseif trial ==3 && cs==2
                    title(['Type Custom Cell response to Current Clamp at ',num2str(10*I(start)),' pA'])
                    subplot(2,1,2)
                    plot(dt:dt:dur,I(1:dur/dt)*10)
                    title('Current Clamp')
                    ylabel('current (pA)')
                    xlabel('time(ms)')
                end
            end
        end
   % end    
%  end
toc

%% Save Code for CV
% cd('E:\Grad School\MATLAB\Sustained-B_0.91_IH conditions\Spikes Data\Spike Trains'); % Change directory for where you wish to save
% filename=input('EnterFileName:','s');
% save(filename,'Outputdata'); %Output data structure is defined in lines near 463.

