% ****************************************
% FIGURE 3C: morph model passive fits
% ****************************************

% run: nrngui neuron-CB-pas-electrode-embed-fit-pas-VClamp.ses
% it uses Neuron's multiple-run parameter fitter. The final fitted output
% is loaded below:

% saved in Neuron rec 10202034, #4 
fit_rec34_t = ...
    trace([ 'nrn-fit-cap-02_dt_0.025000ms_dy_1e-9nA.bin'], ...
          25e-6, 1e-12, 'Fit', ...
          struct('file_type', 'neuron', 'unit_y', 'A', 'scale_y', 1e3));
plotFigure(superposePlots(plotData(...
  [trace([ '../data/10202002.abf' ]), fit_rec34_t], '', ...
  struct('fixedSize', [3 2.5], 'noTitle', 1, 'ColorOrder', [0 0 0; 1 0 0], ...
         'legendLocation', 'southeast', 'axisLimits', [5.5 9.3 NaN 10]))))

% ****************************************
% FIGURE 4B: voltage traces in response to 50pA current injection
% ****************************************

% use neuron-CB-pas-electrode-embed.ses. Load iclamp-50pA.ses and
% v-graph.ses to input stimulus and draw graph. Save the graph
% directly as EPS file. 

% ****************************************
% FIGURE 4C: electrotonic length
% ****************************************

% load neuron-CB-pas-electrode-embed-IClamp.ses and use
% Tools->Impedance->Shape to make plots for 0 Hz vs 100 Hz. 

% ****************************************
% FIGURE 4E/F/H: VC electrode, exp-decaying CC at various locations
% ****************************************

% Uncomment the section injecting exponential current stimulus in
% neuron-CB-pas-electrode-embed-IClamp.ses and load it in
% nrngui. Adjust the location of the input from the
% PointProcessManager window that will open.

% show voltage graphs for panel F

% change exp-decay time constant for panel H using expDecayIClamp
% function.

% ****************************************
% FIGURE 5A: VC of Kdr channels in morphology
% ****************************************

% started with 
% nrngui neuron-CB-act-electrode-embed-IClamp.ses

% changed PointProcess to VClamp and setup a 50 ms VC at 0mV
% set the Tstop at 45 ms

% load shape-plot.ses. Change to "Shape" plot and make sure set the
% scale to "-40 0".

% change location of Kdr channels using CellBuilder Biophysics tab.

% ****************************************
% FIGURE 5B: aCC/RP2 input resistance when Cd present
% ****************************************

% load recorded data 
load('../data/db_2011-02-04_cell1_aCC.mat')
load('../data/db_2011-02-04_cell2_aCC.mat')

% add bias to cip values (do this only one time after loading!)
feb04_cell2_db(:, 'cip_level_pA') = ...
    feb04_cell2_db(:, 'cip_level_pA') + 12;
% TODO: don't do additional correction below! How about cell1?

% default to Cd2+ traces only
cell1_db = feb04_cell1_db(feb04_cell1_db(:, 'cadmium') == 1, :);
cell2_db = feb04_cell2_db(feb04_cell2_db(:, 'cadmium') == 1, :);

% f-V
plotFigure(plot_superpose({...
    set(setProp(plot_cip_stats(cell1_db, '2010/02/04 #1 RP2', ...
                               'PulsePotAvg', 'PulseSpikeRateISI'), ...
                'fixedSize', [3 3], 'axisLimits', [-60 0 NaN NaN], ...
                'legendLocation', 'NorthOutside'), ...
        'axis_labels', {'membrane voltage [mV]', 'Cd^{2+} firing rate [Hz]'}), ...
    plot_cip_stats(cell2_db, '2010/02/04 #2 aCC', ...
                   'PulsePotAvg', 'PulseSpikeRateISI')}));
print -depsc2 current-clamp-firing-marley-all-Cd-cells-f-V.eps

% I-V
plot_iv = plot_superpose({...
    set(setProp(plot_cip_stats(cell1_db, '2010/02/04 #1 RP2', ...
                               'PulsePotAvg', 'cip_level_pA'), ...
                'fixedSize', [3 3], 'axisLimits', [-60 0 NaN NaN]), ...
        'axis_labels', {'membrane voltage [mV]', 'current [pA]'}), ...
    plot_cip_stats(cell2_db, '2010/02/04 #2 aCC', ...
                 'PulsePotAvg', 'cip_level_pA')});
%               'legendLocation', 'NorthOutside'
plotFigure(plot_iv);
print -depsc2 current-clamp-firing-marley-all-Cd-cells-I-V.eps

% do a I-V derivative plot
% TODO: put these in post-columns
cell1_db = ...
    addColumns(cell1_db, ...
               {'I_V_rate_pA_per_mV'}, ...
               [0; diff(cell1_db(:, 'cip_level_pA').data) ./ diff(cell1_db(:, 'PulsePotAvg').data)]);
g1_db = groupBy(cell1_db(:, {'TraceNum', 'PulsePotAvg', 'cip_level_pA'}), ...
                'TraceNum');
g1_diff_db = diff(g1_db);

g1_data = [g1_db(1:(end-1), 'PulsePotAvg', :).data, ...
                    g1_diff_db(:, 'cip_level_pA', :).data ./ ...
                    g1_diff_db(:, 'PulsePotAvg', :).data];
g1_rate_db = ...
    tests_3D_db(reshape(permute(g1_data, [2 1 3]), ...
                        [size(g1_data, 2), size(g1_data, 1) * size(g1_data, 3)])', ...
                {'PulsePotAvg', 'I_V_rate_pA_per_mV'}, {}, {}, ...
                ['I/V rates, ' cell1_db.id]);

cell2_db = ...
    addColumns(cell2_db, ...
               {'I_V_rate_pA_per_mV'}, ...
               [0; diff(cell2_db(:, 'cip_level_pA').data) ./ diff(cell2_db(:, 'PulsePotAvg').data)]);
g2_db = groupBy(cell2_db(:, {'TraceNum', 'PulsePotAvg', 'cip_level_pA'}), ...
                'TraceNum');
g2_diff_db = diff(g2_db);

g2_data = [g2_db(1:(end-1), 'PulsePotAvg', :).data, ...
                    g2_diff_db(:, 'cip_level_pA', :).data ./ ...
                    g2_diff_db(:, 'PulsePotAvg', :).data];
g2_rate_db = ...
    tests_3D_db(reshape(permute(g2_data, [2 1 3]), ...
                        [size(g2_data, 2), size(g2_data, 1) * size(g2_data, 3)])', ...
                {'PulsePotAvg', 'I_V_rate_pA_per_mV'}, {}, {}, ...
                ['I/V rates, ' cell2_db.id]);

% pA/mV = nS
plot_diff_iv = ...
    plot_superpose({...
        set(setProp(plotScatter(g1_rate_db, ...
                                'PulsePotAvg', 'I_V_rate_pA_per_mV', ...
                                '', '2010/02/04 #1 RP2', ...
                                struct('LineStyle', '+')), ...
                    'noTitle', 1, ...
                    'plotProps', struct('LineWidth', 2), ...
                    'fixedSize', [3 3], 'axisLimits', [-80 0 NaN NaN]), ...
            'axis_labels', {'membrane voltage [mV]', 'input conductance [nS]'}), ...
        plotScatter(g2_rate_db, ...
                    'PulsePotAvg', 'I_V_rate_pA_per_mV', '', '2010/02/04 #2 aCC', ...
                    struct('LineStyle', 'o'))});
%                    'legendLocation', 'NorthOutside'
plotFigure(plot_diff_iv)
print -depsc2 current-clamp-firing-marley-all-Cd-cells-I-V-deriv.eps

% make a stack of IV and nS
plotFigure(plot_stack({plot_iv, plot_diff_iv}, [-60 0 NaN NaN], 'y', ...
                      '', struct('fixedSize', [3 5], 'xLabelsPos', ...
                                 'bottom', 'xTicksPos', 'bottom')));
print -depsc2 compare-input-res-K-chans-aCC-RP2.eps

% ****************************************
% FIGURE 6F-I: Changing locations of Na-K channels in morphology
% ****************************************

% different cases are simulated using
% exp-axon-tail2-chans-axon-last.ses: Last reconstructed section of
% axon. 
% exp-axon-tail2-chans-in-all.ses: All of reconstructed morphology.
% exp-axon-tail2-chans-botdend.ses: Bottom (ipsilateral) dendrite.
% exp-axon-tail2-chans-axon.ses: All of reconstructed axon.
% exp-axon-tail2-chans-ext-axon-70um.ses: Added 2-piece axon and channels
% at distance of 70um.

% use the nrniv window (created with Tools->Misc->Family) to start a
% current injection series with 11 steps. 

% Save the simulated data using the saveLines function defined in
% fitfuncs.hoc, such as in: 
% saveLines(2, 1, 11, "data-axon-tail2-chans-botdend")

% Note that you will need to first run the model without stimulation for
% several hundred ms to reach its steady-state and then run finitialize()
% and then saveState() to save it into a file. This file will be loaded at
% each time init() is called.

% plotted in matlab:

% chans in all case
% ********************

chans_in_all_t = ...
    trace([ 'data-axon-tail2-chans-in-all_11_lines_dt_0.025000ms.bin'], ...
          25e-6, 1e-3, 'Neuron axon tail2 - chans in all', ...
          struct('file_type', 'neuron', 'unit_y', 'V', 'scale_y', 1, ...
                 'channel', 1)); % , 'numTraces', 11 (better use trace2cc below)
plot(chans_in_all_t)

% make a CC obj
chans_in_all_cc = ...
    trace2cc(chans_in_all_t, ...
             [50 550], 0:0.005:0.05, ...
             struct('threshold', 1, 'downThreshold', -1, ...
                    'minInit2MaxAmp', 1, 'minMin2MaxAmp', 1, ...
                    'Ihold', 0, 'paramsStruct', struct));
plot(chans_in_all_cc)

% make a tiny figure w current injections 0, 10, 30, 50pA
plotFigure(plot_abstract(setLevels(chans_in_all_cc, [3 7 11]), ...
                         ' - 10, 30, 50 pA', ...
                         struct('fixedSize', [2 2], ...
                                'noTitle', 1, ...
                                'relativeSizes', [1 4], ...
                                'curUnit', 'pA', ...
                                'xTicksPos', 'bottom', ...
                                'axisLimits', [0 600 NaN NaN])))
print('-depsc2', [ morph_analdir 'cc_plots/plot-axon-tail2-chans-in-all_3steps.eps']);

% extract characteristics
[a_db profs ] = params_tests_db(chans_in_all_cc)
displayRows(a_db(:, {'cip_level_pA', 'PulseSpikeRateISI'}))

% an f-I plot
plotFigure(set(plotScatter(a_db, ...
                           'cip_level_pA', 'PulseSpikeRateISI', 'test', 'model', ...
                           struct('LineStyle', '-', 'fixedSize', [2.5 2], 'noTitle', 1)), ...
                   'axis_labels', {'current [pA]', 'firing rate [Hz]'}))


% chans in bot dend
% ********************

chans_botdend_cc = ...
    trace2cc(trace([ 'data-axon-tail2-chans-botdend_11_lines_dt_0.025000ms.bin'], ...
          25e-6, 1e-3, 'Neuron axon tail2 - chans in bot dend', ...
          struct('file_type', 'neuron', 'unit_y', 'V', 'scale_y', 1, ...
                 'channel', 1)), ...
             [50 550], 0:0.005:0.05, ...
             struct('threshold', 1, 'downThreshold', -1, ...
                    'minInit2MaxAmp', 1, 'minMin2MaxAmp', 1, ...
                    'Ihold', 0, 'paramsStruct', struct));
plot(chans_botdend_cc)

% make a tiny figure w current injections 0, 10, 30, 50pA
plotFigure(plot_abstract(setLevels(chans_botdend_cc, [3 7 11]), ...
                         ' - 10, 30, 50 pA', ...
                         struct('fixedSize', [2 2], ...
                                'noTitle', 1, ...
                                'relativeSizes', [1 4], ...
                                'curUnit', 'pA', ...
                                'xTicksPos', 'bottom', ...
                                'axisLimits', [0 600 NaN NaN])))
% select v axis
a = axis;
a(3:4) = [-70 50];
axis(a);
print('-depsc2', [ morph_analdir 'cc_plots/plot-axon-tail2-chans-botdend_3steps.eps']);

% chans in axon
% ********************

chans_axon_cc = ...
    trace2cc(trace([ 'data-axon-tail2-chans-axon_11_lines_dt_0.025000ms.bin'], ...
          25e-6, 1e-3, 'Neuron axon tail2 - chans in axon', ...
          struct('file_type', 'neuron', 'unit_y', 'V', 'scale_y', 1, ...
                 'channel', 1)), ...
             [50 550], 0:0.005:0.05, ...
             struct('threshold', 1, 'downThreshold', -1, ...
                    'minInit2MaxAmp', 1, 'minMin2MaxAmp', 1, ...
                    'Ihold', 0, 'paramsStruct', struct));
plot(chans_axon_cc)

% make a tiny figure w current injections 0, 10, 30, 50pA
plotFigure(plot_abstract(setLevels(chans_axon_cc, [3 7 11]), ...
                         ' - 10, 30, 50 pA', ...
                         struct('fixedSize', [2 2], ...
                                'noTitle', 1, ...
                                'relativeSizes', [1 4], ...
                                'curUnit', 'pA', ...
                                'xTicksPos', 'bottom', ...
                                'axisLimits', [0 600 NaN NaN])))
% select v axis
set(gca, 'YLim', [-70 50]);
print('-depsc2', [ morph_analdir 'cc_plots/plot-axon-tail2-chans-axon_3steps.eps']);

% chans in axon last segment
% ******************************

chans_axon_last_cc = ...
    trace2cc(trace([ 'data-axon-tail2-chans-axon-last_11_lines_dt_0.025000ms.bin'], ...
          25e-6, 1e-3, 'Neuron axon tail2 - chans in axon last segment', ...
          struct('file_type', 'neuron', 'unit_y', 'V', 'scale_y', 1, ...
                 'channel', 1)), ...
             [50 550], 0:0.005:0.05, ...
             struct('threshold', 1, 'downThreshold', -1, ...
                    'minInit2MaxAmp', 1, 'minMin2MaxAmp', 1, ...
                    'Ihold', 0, 'paramsStruct', struct));
plot(chans_axon_last_cc)

% make a tiny figure w current injections 0, 10, 30, 50pA
plotFigure(plot_abstract(setLevels(chans_axon_last_cc, [3 7 11]), ...
                         ' - 10, 30, 50 pA', ...
                         struct('fixedSize', [2 2], ...
                                'noTitle', 1, ...
                                'relativeSizes', [1 4], ...
                                'curUnit', 'pA', ...
                                'xTicksPos', 'bottom', ...
                                'axisLimits', [0 600 NaN NaN])))
% select v axis
set(gca, 'YLim', [-70 50]);
print('-depsc2', [ morph_analdir 'cc_plots/plot-axon-tail2-chans-axon-last_3steps.eps']);

% chans in ext axon, 70 um from last axon segment
% **************************************************

chans_ext_axon_70um_cc = ...
    trace2cc(trace([ 'data-axon-tail2-chans-ext-axon-70um_11_lines_dt_0.025000ms.bin'], ...
          25e-6, 1e-3, 'Neuron axon tail2 - chans in ext axon 70 um from last axon segment', ...
          struct('file_type', 'neuron', 'unit_y', 'V', 'scale_y', 1, ...
                 'channel', 1)), ...
             [50 550], 0:0.005:0.05, ...
             struct('threshold', 1, 'downThreshold', -1, ...
                    'minInit2MaxAmp', 1, 'minMin2MaxAmp', 1, ...
                    'Ihold', 0, 'paramsStruct', struct));
plot(chans_ext_axon_70um_cc)

% make a tiny figure w current injections 0, 10, 30, 50pA
plotFigure(plot_abstract(setLevels(chans_ext_axon_70um_cc, [3 7 11]), ...
                         ' - 10, 30, 50 pA', ...
                         struct('fixedSize', [2 2], ...
                                'noTitle', 1, ...
                                'relativeSizes', [1 4], ...
                                'curUnit', 'pA', ...
                                'xTicksPos', 'bottom', ...
                                'axisLimits', [0 600 NaN NaN])))
% select v axis
set(gca, 'YLim', [-70 50]);
print('-depsc2', [ morph_analdir 'cc_plots/plot-axon-tail2-chans-ext-axon-70um_3steps.eps']);


% ****************************************
% FIGURE 6J: aCC current reponse
% ****************************************

% same as in fig 2

% ****************************************
% FIGURE 6K,L: voltage offset and spike amp comparison
% ****************************************

% another plot showing spike height vs pA
% ********************

% extract characteristics
[chans_in_all_db] = params_tests_db(chans_in_all_cc);
[chans_botdend_db] = params_tests_db(chans_botdend_cc);
[chans_axon_db] = params_tests_db(chans_axon_cc);
[chans_axon_last_db] = params_tests_db(chans_axon_last_cc);
[chans_ext_axon_70um_db] = params_tests_db(chans_ext_axon_70um_cc);

% func for making line plot
plot_spikeamp = @(a_db, a_label) ...
    set(plotScatter(a_db, ...
                    'cip_level_pA', 'PulseSpikeAmplitudeMean', 'Spike Amp', a_label, ...
                    struct('LineStyle', '-', 'noTitle', 1)), ...
        'axis_labels', {'current [pA]', 'mean spike amplitude [mV]'});

% compare spike amps together
plotFigure(plot_superpose({...
    setProp(plot_spikeamp(chans_in_all_db, 'all'), ...
            'fixedSize', [4 2], 'legendLocation', 'EastOutside', ...
            'axisLimits', [0 50 0 80], ...
            'axisProps', struct('Box', 'off'), ...
            'ColorOrder', [0 0 1; 0 0.5 0; 1 0 0; ...
                    [232 170 228]./255; [191 0 191]./255; [0 0 0]]), ...
    plot_spikeamp(chans_botdend_db, 'botdend'), ...
    plot_spikeamp(chans_axon_db, 'axon'), ...
    plot_spikeamp(chans_axon_last_db, 'axonlast'), ...
    plot_spikeamp(chans_ext_axon_70um_db, 'ext70um') ...
    plot_cip_stats(cell2_db, 'aCC w/ Cd^{2+}', ...
                   'cip_level_pA', 'PulseSpikeAmplitudeMean')
                          }))
print('-depsc2', [ 'plot-chan-dists-compare-spike-amp.eps']);

% do also for offset
plot_offset = @(a_db, a_label) ...
    set(plotScatter(a_db, ...
                    'cip_level_pA', 'PulsePotAvg', 'Voltage Offset', a_label, ...
                    struct('LineStyle', '-', 'noTitle', 1)), ...
        'axis_labels', {'current [pA]', 'voltage offset [mV]'});

plotFigure(plot_superpose({...
    setProp(plot_offset(chans_in_all_db, 'all'), ...
            'axisLimits', [0 50 -65 0], ...
            'fixedSize', [2.5 2], 'noLegends', 1, ...
            'axisProps', struct('Box', 'off'), ...
            'ColorOrder', [0 0 1; 0 0.5 0; 1 0 0; ...
                    [232 170 228]./255; [191 0 191]./255; [0 0 0]]), ...
    plot_offset(chans_botdend_db, 'botdend'), ...
    plot_offset(chans_axon_db, 'axon'), ...
    plot_offset(chans_axon_last_db, 'axonlast'), ...
    plot_offset(chans_ext_axon_70um_db, 'ext70um') ...
    plot_cip_stats(cell2_db, 'aCC w/ Cd^{2+}', ...
                   'cip_level_pA', 'PulsePotAvg') ...
                   }))
print('-depsc2', [ 'plot-chan-dists-compare-offset.eps']);

% ************************************************************
% Figure 7: NaT and NaP components affected by morphology
% ************************************************************

% load exp-axon-tail2-chans-ext-axon-70um-onlyNa.ses

% ********************
% Panel 7A: Voltage distribution across morphology
% ********************

% save voltages with:
% saveLines(1, 1, 14, "data-axon-tail2-axon-70um-vc-Na-Kdr-long-back-85mV-justV")

% $$$ oc>readGraph(1)
% $$$ #0 label=electrode 
% $$$ #2 label=soma 
% $$$ #4 label=axon (prox) 
% $$$ #6 label=axon (mid) 
% $$$ #8 label=axon (dist) 
% $$$ #10 label=top dend 
% $$$ #12 label=bot dend 
% $$$ #14 label=axon ext (dist) 

% $$$ 

vc_t = ...
    trace([ morph_analdir 'data-axon-tail2-axon-70um-vc-Na-Kdr-long-back-85mV-justV_14_lines_dt_0.025000ms.bin'], ...
          25e-6, 1e-3, 'VC volt-dist passive w/ ext axon', ...
          struct('file_type', 'neuron', 'unit_y', 'V', 'scale_y', 1, ...
                 'channel', 1))

% reshape it
vc_t.data = reshape(vc_t.data, size(vc_t.data, 1) / 14, 14);

% select few representative v traces
vc_t.data = vc_t.data(:, [1 2 7 13]);

plotFigure(set(plot_abstract(vc_t, '', ...
                             struct('fixedSize', [3.5 2], ...
                                    'legendLocation', 'EastOutside', ...
                                    'noTitle', 1, ...
                                    'axisLimits', [50 110 -87 50], ...
                                    'axisProps', struct('Box', 'off'))), ...
               'legend', ...
               {'electrode', 'soma', 'dist. axon', 'ext. axon'}))

%print('-depsc2', [ morph_analdir 'cc_plots/plot-Na-vc-Na-Kdr-volt-dists.eps']);
print('-depsc2', [ morph_analdir 'cc_plots/plot-Na-vc-noKdr-volt-dists.eps']);

% ********************
% Panel 7A: Na currents throughout the morphology
% ********************

% to get currents, save voltages with passive and active properties, and
% then subtract one from the other.

% run this to save both passive and active traces:
% oc> runSave("exp-name")
% where "exp-name" is a keyword you want to append to file names.

load_pas_vt = @(name)...
          trace([ name ], ...
          25e-6, 1e-3, 'VC passive w/ ext axon', ... 
         struct('file_type', 'neuron', 'unit_y', 'V', 'scale_y', 1, ...
                 'channel', 1));

load_act_vt = @(name)...
          trace([ name ], ...
          25e-6, 1e-3, 'VC active w/ ext axon', ... 
         struct('file_type', 'neuron', 'unit_y', 'V', 'scale_y', 1, ...
                 'channel', 1));

% axon13 to axon12 axial resistance (attached to soma)
ri = 25.045413 + 6.2445732; % MOhms

% first pas
vaxon12_pas_t = ...
    load_pas_vt('data-ext-axon-70um-onlyNa_axon[12](0.5)_pas_dt_0.025000ms.bin');
vaxon13_pas_t = ...
    load_pas_vt('data-ext-axon-70um-onlyNa_axon[13](0.5)_pas_dt_0.025000ms.bin');
iaxon12_pas_t = (vaxon13_pas_t - vaxon12_pas_t) ./ ri; % mV/MO=nA
iaxon12_pas_t.props.unit_y='A';
%plot(iaxon12_pas_t)

% then act
vaxon12_act_t = ...
    load_act_vt('data-ext-axon-70um-onlyNa_axon[12](0.5)_act_dt_0.025000ms.bin');
vaxon13_act_t = ...
    load_act_vt('data-ext-axon-70um-onlyNa_axon[13](0.5)_act_dt_0.025000ms.bin');
iaxon12_act_t = (vaxon13_act_t - vaxon12_act_t) ./ ri; % mV/MO=nA
iaxon12_act_t.props.unit_y='A';
%plot(iaxon12_act_t)

% take the difference in current
iaxon12_t = iaxon12_act_t - iaxon12_pas_t;
iaxon12_t.props.unit_y='A';
iaxon12_t.dy=1e-9;
plot(iaxon12_t)

% axon6 to axon7 (middle of ipsilateral dendritic tree)
ri = 1.8150921 + 0.83369277; % MOhms

% first pas
vaxon6_pas_t = ...
    load_pas_vt('data-ext-axon-70um-onlyNa_axon[6](0.5)_pas_dt_0.025000ms.bin');
vaxon7_pas_t = ...
    load_pas_vt('data-ext-axon-70um-onlyNa_axon[7](0.5)_pas_dt_0.025000ms.bin');
iaxon6_pas_t = (vaxon7_pas_t - vaxon6_pas_t) ./ ri; % mV/MO=nA
iaxon6_pas_t.props.unit_y='A';
plot(iaxon6_pas_t)

% then act
vaxon6_act_t = ...
    load_act_vt('data-ext-axon-70um-onlyNa_axon[6](0.5)_act_dt_0.025000ms.bin');
vaxon7_act_t = ...
    load_act_vt('data-ext-axon-70um-onlyNa_axon[7](0.5)_act_dt_0.025000ms.bin');
iaxon6_act_t = (vaxon7_act_t - vaxon6_act_t) ./ ri; % mV/MO=nA
iaxon6_act_t.props.unit_y='A';
plot(iaxon6_act_t)

% take the difference in current
iaxon6_t = iaxon6_act_t - iaxon6_pas_t;
iaxon6_t.props.unit_y='A';
iaxon6_t.dy=1e-9;
plot(iaxon6_t)

% start of axon_ext[1] from 0.02 to 0.06 (SIZ location)
ri = 98.22134; % MOhms
vaxon_ext1a_pas_t = ...
    load_pas_vt('data-ext-axon-70um-onlyNa_axon_ext[1](0.02)_pas_dt_0.025000ms.bin');
vaxon_ext1b_pas_t = ...
    load_pas_vt('data-ext-axon-70um-onlyNa_axon_ext[1](0.06)_pas_dt_0.025000ms.bin');
iaxon_ext1_pas_t = (vaxon_ext1a_pas_t - vaxon_ext1b_pas_t) ./ ri; % mV/MO=nA
iaxon_ext1_pas_t.props.unit_y='A';
plot(iaxon_ext1_pas_t)

% act
vaxon_ext1a_act_t = ...
    load_act_vt('data-ext-axon-70um-onlyNa_axon_ext[1](0.02)_act_dt_0.025000ms.bin');
vaxon_ext1b_act_t = ...
    load_act_vt('data-ext-axon-70um-onlyNa_axon_ext[1](0.06)_act_dt_0.025000ms.bin');
iaxon_ext1_act_t = (vaxon_ext1a_act_t - vaxon_ext1b_act_t) ./ ri; % mV/MO=nA
iaxon_ext1_act_t.props.unit_y='A';
plot(iaxon_ext1_act_t)

% sub
iaxon_ext1_t = iaxon_ext1_act_t - iaxon_ext1_pas_t;
iaxon_ext1_t.props.unit_y='A';
iaxon_ext1_t.dy=1e-9;
plot(iaxon_ext1_t)

% start of axon_ext[0] from 0.3 to 0.5 (extended axon - 70 um)
ri = 52.640496; % MOhms
vaxon_ext0a_pas_t = ...
    load_pas_vt('data-ext-axon-70um-onlyNa_axon_ext[0](0.3)_pas_dt_0.025000ms.bin');
vaxon_ext0b_pas_t = ...
    load_pas_vt('data-ext-axon-70um-onlyNa_axon_ext[0](0.5)_pas_dt_0.025000ms.bin');
iaxon_ext0_pas_t = (vaxon_ext0a_pas_t - vaxon_ext0b_pas_t) ./ ri; % mV/MO=nA
iaxon_ext0_pas_t.props.unit_y='A';
iaxon_ext0_pas_t.dy=1e-9;
plot(iaxon_ext0_pas_t)

% act
vaxon_ext0a_act_t = ...
    load_act_vt('data-ext-axon-70um-onlyNa_axon_ext[0](0.3)_act_dt_0.025000ms.bin');
vaxon_ext0b_act_t = ...
    load_act_vt('data-ext-axon-70um-onlyNa_axon_ext[0](0.5)_act_dt_0.025000ms.bin');
iaxon_ext0_act_t = (vaxon_ext0a_act_t - vaxon_ext0b_act_t) ./ ri; % mV/MO=nA
iaxon_ext0_act_t.props.unit_y='A';
iaxon_ext0_act_t.dy=1e-9;
plot(iaxon_ext0_act_t)

% sub
iaxon_ext0_t = iaxon_ext0_act_t - iaxon_ext0_pas_t;
plot(iaxon_ext0_t)

% manicured plot version
plotFigure(plot_superpose({...
  plot_abstract(iaxon12_t, 'axon[12]', ...
                struct('quiet', 1, 'fixedSize', [2.5 2], ...
                       'noTitle', 1, ...
                       'axisLimits', [54 70 -.25 0.01], ...
                       'plotProps', struct('LineWidth', 2), ...
                       'axisProps', struct('Box', 'off'))), ...
  plot_abstract(iaxon6_t, 'axon[6]', struct('quiet', 1)), ...
  plot_abstract(iaxon_ext0_t, 'axon_ext[0]', struct('quiet', 1))}))

print('-depsc2', [ morph_analdir 'vc_plots/plot-vc-ext-axon-70um-onlyNa-85mV_to_-25mV_compare_3compartments.eps']);

% ********************
% Panel 7B: Na current changes with increases in NaT vs NaP conductances
% ********************

% In Neuron, apply NaP/T changes for given current injection

% this is the passive case so that we can use it for capacitance artifact
% subtraction. It's a VC step from -85 mV to -25 mV like in Marley's L3 recording.
na_pas_t = ...
    trace([ 'data-axon-tail2-axon-70um-vc-noKdr-long-back-85mV-passive_4_lines_dt_0.025000ms.bin'], ...
          25e-6, 1e-9, 'VC passive w/ ext axon', ... 
         struct('file_type', 'neuron', 'unit_y', 'A', 'scale_y', 1, ...
                 'channel', 1))

% regular Na channels (w/o Kdr it stays depolarized unless it's properly initialized)
na_t = ...
    trace([ 'data-axon-tail2-axon-70um-vc-noKdr-long-back-85mV-Na_4_lines_dt_0.025000ms.bin'], ...
          25e-6, 1e-9, 'VC to get Na currents - chans 70 um into ext axon', ...
          struct('file_type', 'neuron', 'unit_y', 'A', 'scale_y', 1, ...
                 'channel', 1))

% subtract passive capacitance artifacts
na_sub_t = na_t - na_pas_t;

% 5xNaT conductance from regular
nat5_t = ...
    trace([ 'data-axon-tail2-axon-70um-vc-noKdr-long-back-85mV-Na-5xNaT_4_lines_dt_0.025000ms.bin'], ...
          25e-6, 1e-9, 'Na-VC chans 70 um into ext axon, 2xNaT', ...
          struct('file_type', 'neuron', 'unit_y', 'A', 'scale_y', 1, ...
                 'channel', 1))

% 5xNaP conductance from regular
nap5_t = ...
    trace([ 'data-axon-tail2-axon-70um-vc-noKdr-long-back-85mV-Na-5xNaP_4_lines_dt_0.025000ms.bin'], ...
          25e-6, 1e-9, 'Na-VC chans 70 um into ext axon, 2xNaP', ...
          struct('file_type', 'neuron', 'unit_y', 'A', 'scale_y', 1, ...
                 'channel', 1))

% moved to 50um axon_ext
na_50um_t = ...
    trace([ 'data-axon-tail2-axon-50um-vc-noKdr-long-back-85mV-Na_4_lines_dt_0.025000ms.bin'], ...
          25e-6, 1e-9, 'VC to get Na currents - chans 50 um into ext axon', ...
          struct('file_type', 'neuron', 'unit_y', 'A', 'scale_y', 1, ...
                 'channel', 1))

% set axis size the same so can use same scalebars
% calc_axis - produces axis() command output
calc_axis = @(x,y,w,h) [x, x+w, y, y+h];

% compare all in one panel
plotFigure(plot_superpose({...
  plot_abstract(na_t - na_pas_t, 'orig.', ...
                struct('fixedSize', [2.5 2], 'noTitle', 1, 'quiet', 1, ...
                       'plotProps', struct('LineWidth', 2), ...
                       'axisProps', struct('Box', 'off'))), ...
  plot_abstract(nat5_t - na_pas_t, '5xNaT', struct('quiet', 1)), ...
  plot_abstract(nap5_t - na_pas_t, '5xNaP', struct('quiet', 1)), ...
  plot_abstract(na_50um_t - na_pas_t, '50\mu{}m', struct('quiet', 1))}))

axis(calc_axis(53, -0.15, 25, 0.16))

%print('-depsc2', [ morph_analdir 'cc_plots/plot-vc-long-back-85mV-NaT-NaP-double.eps']);
%print('-depsc2', [ morph_analdir 'cc_plots/plot-vc-long-back-85mV-NaT-NaP-5x.eps']);
print('-depsc2', [ morph_analdir 'cc_plots/plot-vc-long-back-85mV-NaT-NaP-5x+50um.eps']);

% ********************
% Panel 7C/D
% ********************

% move axis to show NaT
axis(calc_axis(275, 0, 5, 0.4))
print('-depsc2', [ morph_analdir 'cc_plots/plot-vc-long-back-85mV-NaT-NaP-5x-show-NaT.eps']);

% move axis to show NaP
axis(calc_axis(380, 0, 25, 0.0007))
print('-depsc2', [ morph_analdir 'cc_plots/plot-vc-long-back-85mV-NaT-NaP-5x-show-NaP.eps']);

% ********************
% Panel 7E: P/T ratio changing with other parameters
% ********************

pt_ratios_db = ...
    tests_db('NaP_NaT_data.csv', {}, {}, ...
             ['Collected data for NaP/NaT '], ...
             struct('csvReadColNames', 1));
% export for paper
string2File(displayRowsTeX(transpose(pt_ratios_db)), 'table-PT-ratios-escape-VC.tex')

border = [0 .03 .1 0];
axis_limits = [0 1e-4 NaN NaN];
color_obsV = [0 0.49 0];
color_PT = [0.87 0.49 0];
plotFigure(plot_inset({...
    plotScatter(pt_ratios_db(1:4, :), 'galeak', 'PTsoma', '', '', ...
                       struct('quiet', 1, 'LineStyle', '-', ...
                              'border', border, ...
                              'axisLimits', axis_limits, ...
                              'plotProps', ...
                              struct('Color', color_PT, ...
                                     'LineWidth', 2), ...
                              'axisProps', ...
                              struct('Box', 'off', ...
                                     'YColor', color_PT))), ...
    plotScatter(pt_ratios_db(1:4, :), 'galeak', 'obsV', '', '', ...
                       struct('quiet', 1, 'LineStyle', '-', ...
                              'border', border, ...
                              'axisLimits', axis_limits, ...
                              'plotProps', ...
                              struct('Color', color_obsV, ...
                                     'LineWidth', 2), ...
                              'axisProps', ...
                              struct('Color', 'none', ...
                                     'YColor', color_obsV, ...
                                     'YAxisLocation', 'right', ...
                                     'Box', 'off')))}, ...
                      [0 0 1 1; 0 0 1 1], '', ...
                      struct('fixedSize', [2.5 2], ...
                             'border', [0 .03 .03 0])))
print -depsc2 plot-scatter-galeak-vs-PTsoma-obsV.eps

% ****************************************
% Figure 8A: compare synaptic mini currents injected at different sites.
% ****************************************

% use exp-axon-tail2-chans-ext-axon-70um-mimic-synapses.ses and
% exp-axon-tail2-chans-ext-axon-70um-mimic-synapses-v-change.ses, which
% is the no-VC version.

% choose axis limits to deo all different synapse locations
axis_limits = [48 72 -12 1];

syn_ext_axon_70um_dend685_t = ...
    trace([ 'data-syn-dend-685_2_lines_dt_0.025000ms.bin'], ...
          25e-6, 1e-12, 'dend685 synapse - Neuron axon tail2 chans 70 um into ext axon', ...
          struct('file_type', 'neuron', 'unit_y', 'A', 'scale_y', 1e3, ...
                 'channel', 1, 'numTraces', 2));
plotFigure(plot_abstract(syn_ext_axon_70um_dend685_t, '', ...
                         struct('fixedSize', [2.5 2], ...
                                'noTitle', 1, ...
                                'axisLimits', axis_limits)))


syn_ext_axon_70um_dend520_t = ...
    trace([ 'data-syn-dend-520_2_lines_dt_0.025000ms.bin'], ...
          25e-6, 1e-12, 'dend520 synapse - Neuron axon tail2 chans 70 um into ext axon', ...
          struct('file_type', 'neuron', 'unit_y', 'A', 'scale_y', 1e3, ...
                 'channel', 1, 'numTraces', 2));
plotFigure(plot_abstract(syn_ext_axon_70um_dend520_t, '', ...
                         struct('fixedSize', [2.5 2], ...
                                'noTitle', 1, ...
                                'axisLimits', axis_limits)))

syn_ext_axon_70um_dend357_t = ...
    trace([ 'data-syn-dend-357_2_lines_dt_0.025000ms.bin'], ...
          25e-6, 1e-12, 'dend357 synapse - Neuron axon tail2 chans 70 um into ext axon', ...
          struct('file_type', 'neuron', 'unit_y', 'A', 'scale_y', 1e3, ...
                 'channel', 1, 'numTraces', 2));
plotFigure(plot_abstract(syn_ext_axon_70um_dend357_t, '', ...
                         struct('fixedSize', [2.5 2], ...
                                'noTitle', 1, ...
                                'axisLimits', axis_limits)))

syn_ext_axon_70um_dend513_t = ...
    trace([ 'data-syn-dend-513_2_lines_dt_0.025000ms.bin'], ...
          25e-6, 1e-12, 'dend513 synapse - Neuron axon tail2 chans 70 um into ext axon', ...
          struct('file_type', 'neuron', 'unit_y', 'A', 'scale_y', 1e3, ...
                 'channel', 1, 'numTraces', 2));
plotFigure(plot_abstract(syn_ext_axon_70um_dend513_t, '', ...
                         struct('fixedSize', [2.5 2], ...
                                'noTitle', 1, ...
                                'axisLimits', axis_limits)))

% plot all of them in stack
plotFigure(plot_stack({plot_abstract(syn_ext_axon_70um_dend685_t, '', ...
                         struct('noTitle', 1)), ...
                    plot_abstract(syn_ext_axon_70um_dend520_t, '', ...
                         struct('noTitle', 1)), ...
                    plot_abstract(syn_ext_axon_70um_dend357_t, '', ...
                         struct('noTitle', 1)), ...
                    plot_abstract(syn_ext_axon_70um_dend513_t, '', ...
                         struct('noTitle', 1))}, axis_limits, 'x', '', ...
                      struct('yTicksPos', 'left', 'yLabelsPos', 'left', ...
                             'fixedSize', [6 2])))
% export this and modify EPS
print('-depsc2', [ morph_analdir 'vc_plots/plot-4-synapses-6pA-epsc-' ...
                   'axon-tail2-chans-ext-axon-70um-iclamp.eps']);

% ****************************************
% Figure 8B: single AP generation
% ****************************************

% load V graph that shows single AP generation
v_syn_ext_axon_70um_dend685_AP_t = ...
    trace([ 'data-v-syn-dend-685-AP_3_lines_dt_0.025000ms.bin'], ...
          25e-6, 1e-3, 'dend685 synapse AP - Neuron axon tail2 chans 70 um into ext axon', ...
          struct('file_type', 'neuron', 'unit_y', 'V', 'scale_y', 1, ...
                 'channel', 1, 'numTraces', 3));
plotFigure(plot_abstract(v_syn_ext_axon_70um_dend685_AP_t, '', ...
                         struct('fixedSize', [2.5 2], ...
                                'noTitle', 1)))
legend('elec', 'syn', 'ext\_axon')
print('-depsc2', [ morph_analdir 'cc_plots/plot-4-synapses-AP-elec-syn-axon-' ...
                   'axon-tail2-chans-ext-axon-70um-iclamp.eps']);

% ****************************************
% Figure 8C: curve of presynaptic stim freq vs output current
% ****************************************

syn1_freq_db = ...
    tests_db([[10, 50, 100, 200, 333, 1000]; ...
              [6,   6,   7,  11,  17,   37]]', ...
             {'presyn_freq_Hz', 'soma_cur_pA'}, {}, ...
             'synapse 685 stim w/ 10 spikes');
plotFigure(plotScatter(syn1_freq_db, 'presyn_freq_Hz', 'soma_cur_pA', ...
                       '', '', ...
                       struct('LineStyle', '-o', 'fixedSize', [2.5 2], ...
                              'noTitle', 1, 'tightLimits', 1)))
print('-depsc2', ...
      [ morph_analdir 'synapse_plots/plot-curve-dend685-num-EPSCs-10-presyn-freq-vs-soma-cur-vc-60mV.eps']);

% ****************************************
% Figure 8D: voltage/current graphs that show saturation at dend513 synapse
% ****************************************

% w/o VC
v_syn_dend513_180_10_saturate_noVC_t = ...
    trace([ 'data-v-syn-dend-513-180-EPSCs-10x-1ms-saturating-noVC_5_lines_dt_0.025000ms.bin'], ...
          25e-6, 1e-3, 'dend513 synapse 180x 10 EPSCs at 1ms interval V, no VC', ...
          struct('file_type', 'neuron', 'unit_y', 'V', 'scale_y', 1, ...
                 'channel', 1, 'numTraces', 5));
plotFigure(plot_abstract(v_syn_dend513_180_10_saturate_noVC_t, '', ...
                         struct('fixedSize', [2.5 2], ...
                                'noTitle', 1)))
legend('soma', 'syn', 'axon')
print('-depsc2', ...
      [ morph_analdir 'synapse_plots/plot-v-syn-axon-dend513-180x-EPSCs-10-at-1ms-each-noVC.eps']);


% graph current
syn_dend513_180_10_saturate_t = ...
    trace([ 'data-i-vclamp-syn-dend-513-180-EPSCs-10x-1ms-saturating_2_lines_dt_0.025000ms.bin'], ...
          25e-6, 1e-12, 'dend513 synapse 180x 10 EPSCs at 1ms interval', ...
          struct('file_type', 'neuron', 'unit_y', 'A', 'scale_y', 1e3, ...
                 'channel', 1, 'numTraces', 2));
plotFigure(plot_abstract(syn_dend513_180_10_saturate_t, '', ...
                         struct('fixedSize', [2.5 2], ...
                                'noTitle', 1)))
legend('vc', 'syn')
% truncated peak at -1981 pA
print('-depsc2', [ morph_analdir 'synapse_plots/plot-syn-vc-dend513-180x-EPSCs-10-at-1ms-each-vc-60mV-100ms.eps']);

% ****************************************
% Figure 9A: 4x synapse case
% ****************************************

% use
% exp-axon-tail2-chans-ext-axon-70um-mimic-synapses-sustained-currents.ses

% voltage
v_4dends_50x_10_EPSCs_interval_10ms_noVC_t = ...
    trace([ 'data-v-syn-4dends-50-EPSCs-10x-10ms-noVC_6_lines_dt_0.025000ms.bin'], ...
          25e-6, 1e-3, '4 dends synapse 50x 10 EPSCs at 10ms interval V, no VC', ...
          struct('file_type', 'neuron', 'unit_y', 'V', 'scale_y', 1, ...
                 'channel', 1, 'numTraces', 6));
plotFigure(plot_abstract(v_4dends_50x_10_EPSCs_interval_10ms_noVC_t, '', ...
                         struct('fixedSize', [2.5 2], ...
                                'noTitle', 1)))
legend('axon', 'soma', 'syn')
print('-depsc2', ...
      [ morph_analdir 'synapse_plots/plot-v-syn-axon-4dend-50x-EPSCs-10-at-10ms-each-noVC.eps']);

% current
i_4dends_50x_10_EPSCs_interval_10ms_noVC_t = ...
    trace([ 'data-i-syn-4dends-50-EPSCs-10x-10ms-VC-60mV_5_lines_dt_0.025000ms.bin'], ...
          25e-6, 1e-12, '4 dends synapse 50x 10 EPSCs at 10ms interval', ...
          struct('file_type', 'neuron', 'unit_y', 'A', 'scale_y', 1e3, ...
                 'channel', 1, 'numTraces', 5));
plotFigure(plot_abstract(i_4dends_50x_10_EPSCs_interval_10ms_noVC_t, '', ...
                         struct('fixedSize', [2.5 2], ...
                                'noTitle', 1)))
legend('I_{VC}', 'syn1', 'syn2', 'syn3', 'syn4')
print('-depsc2', [ morph_analdir 'synapse_plots/plot-i-vc-4dends-50x-EPSCs-10-at-10ms-each-vc-60mV.eps']);

% ****************************************
% Figure 9B: 10x synapse case with exactly 10 spikes @ 100 Hz
% ****************************************

% use exp-axon-tail2-chans-ext-axon-70um-10x-mimic-sustained.ses, but no
% randomization.

% voltage
v_10dends_20x_EPSCs_10_interval_10ms_noVC_t = ...
    trace([ 'data-v-syn-10dends-20-EPSCs-10x-10ms-noVC_6_lines_dt_0.025000ms.bin'], ...
          25e-6, 1e-3, '10 dends synapse 20x 10 EPSCs at 10ms interval V, no VC', ...
          struct('file_type', 'neuron', 'unit_y', 'V', 'scale_y', 1, ...
                 'channel', 1, 'numTraces', 6));
plotFigure(plot_abstract(v_10dends_20x_EPSCs_10_interval_10ms_noVC_t, '', ...
                         struct('fixedSize', [2.5 2], ...
                                'noTitle', 1)))
legend('axon', 'bottom dend', 'top dend', 'soma')
print('-depsc2', ...
      [ morph_analdir 'synapse_plots/plot-v-syn-axon-10dend-20x-EPSCs-10-at-10ms-each-noVC.eps']);

% current
i_10syns_20x_EPSCs_10x_interval_10ms_VC_t = ...
    trace([ 'data-i-syn-10syns-20-EPSCs-10x-10ms-VC-60mV_6_lines_dt_0.025000ms.bin'], ...
          25e-6, 1e-12, '10 synapses 20x EPSCs 10 at 10ms interval', ...
          struct('file_type', 'neuron', 'unit_y', 'A', 'scale_y', 1e3, ...
                 'channel', 1, 'numTraces', 6));
plotFigure(plot_abstract(i_10syns_20x_EPSCs_10x_interval_10ms_VC_t, '', ...
                         struct('fixedSize', [2.5 2], ...
                                'noTitle', 1)))
legend('I_{VC}', 'syn1', 'syn3', 'syn5', 'syn7', 'syn9')
print('-depsc2', [ morph_analdir 'synapse_plots/plot-i-vc-10dends-20x-EPSCs-10-at-10ms-each-vc-60mV.eps']);

% ****************************************
% Figure 9D: EPSCx vs I curves
% ****************************************

% make curve of 1 synapse 180x @ 1 kHz (saturation)
syn1_1kHz_db = ...
    tests_db([[10, 20, 50, 120, 180]; [200, 330, 380, 405, 418]]', ...
             {'syn_mag_times', 'syn_cur_pA'}, {}, ...
             'synapse 513 saturation w/ 10 spikes @ 1kHz');
plotFigure(plotScatter(syn1_1kHz_db, 'syn_mag_times', 'syn_cur_pA', ...
                       '', '', ...
                       struct('LineStyle', '-o', 'fixedSize', [2.5 2], ...
                              'noTitle', 1, 'tightLimits', 1)))
print('-depsc2', ...
      [ morph_analdir 'synapse_plots/plot-curve-dend513-num-EPSCs-10-at-1ms-each-vs-peak-syn-cur-vc-60mV.eps']);

% curve for 4 synapses 50x @ 100 Hz
syn4_100Hz_db = ...
    tests_db([4*[10, 20, 30, 40, 50]; [166, 275, 347, 408, 452]]', ...
             {'num_EPSCs', 'soma_cur_pA'}, {}, ...
             '4 synapses w/ shifted 10 spikes @ 100Hz');

% data for 10 spikes at 10 synapses @ 100 Hz
syn10_100Hz_db = ...
    tests_db([10*[1 5 10 15 20 25 30 35]; [46 194 324 395 462 511 552 583]]', ...
             {'num_EPSCs', 'soma_cur_pA'}, {}, ...
             '10 synapses w/ shifted 10 spikes @ 100Hz');

plotFigure(superposePlots([plotScatter(syn1_1kHz_db, 'syn_mag_times', 'syn_cur_pA', ...
                       '', '1 syn @ 1kHz', ...
                       struct('LineStyle', '-o', 'fixedSize', [2.5 2], ...
                              'noTitle', 1, 'axisLimits', [0 350 0 650])), ...
                    plotScatter(syn4_100Hz_db, 'num_EPSCs', 'soma_cur_pA', ...
                       '', '4 syn @ 100Hz', ...
                       struct('LineStyle', '-+')), ...
                    plotScatter(syn10_100Hz_db, 'num_EPSCs', 'soma_cur_pA', ...
                       '', '10 syn @ 100Hz', ...
                       struct('LineStyle', '-d'))]))

print('-depsc2', ...
      [ morph_analdir 'synapse_plots/plot-curve-num-EPSCs-vs-peak-soma-cur-all-conds-vc-60mV.eps']);

% ========================================
% Figure 9E: input frequency with 10 synapses
% ========================================

syn10_freq_db = ...
    tests_db([[10 50 100 200]; ...
              [106 328 462 634]]', ...
             {'presyn_freq_Hz', 'soma_cur_pA'}, {}, ...
             '10 synapses 20x EPSC stim w/ 10 spikes');
plotFigure(plotScatter(syn10_freq_db, 'presyn_freq_Hz', 'soma_cur_pA', ...
                       '', '', ...
                       struct('LineStyle', '-o', 'fixedSize', [2.5 2], ...
                              'noTitle', 1, 'tightLimits', 1)))
axis([0 200 0 650])
legend('10 syn, 20x EPSC')
print('-depsc2', ...
      [ morph_analdir 'synapse_plots/plot-curve-10dend-num-EPSCs-20-presyn-freq-vs-soma-cur-vc-60mV.eps']);

% ========================================
% Figure 9F: SRCs
% ========================================

% for model; use exp-axon-tail2-chans-ext-axon-70um-10x-mimic-sustained.ses,
% with randomization.

% recorded SRC

synapse_SRCL3_datadir = ...
    [ '../data/' ];

src1_t = trace([synapse_SRCL3_datadir '12829006.abf'])
plot(src1_t, '', struct('timeScale', 's'))
% @ 9.8s: mag 500 pA, 600 ms
plot(src1_t, '', ...
     struct('timeScale', 's', ...
            'axisLimits', ...
            [9.7 10.7 -540 40], ...
            'fixedSize', [2.5 2], 'noTitle', 1))
print('-depsc2', [synapse_analdir 'src-12829006-t9.800-500pA-600ms.eps'] )